獲取Github程式碼包以及準備工作

  • 2020 年 3 月 31 日
  • 筆記

先下載程式碼包

github程式碼在：https://github.com/jmzeng1314/scRNA_smart_seq2/archive/master.zip

首先進入RNA-seq目錄，從step0-step9是對常規轉錄組的一個回顧

準備工作之R包

從step0開始，程式碼注釋蠻詳細的，我會挑選重要的部分寫到這裡，其他可以自行看程式碼學習，下面就是主要利用Rstudio進行操作了

  • 一個好習慣，做新項目時記得清空之前的變數，使用rm(list = ls())
  • 有的R包比較大，經常需要載入其他的動態庫dynamically loaded libraries (DLLs)，例如： > length(loadedNamespaces()) [1] 34 > library(Seurat) #載入一個seurat包會出現接近60個依賴的動態庫 > length(loadedNamespaces()) [1] 96 如果不設置，就會因為載入數量超限制而報錯(https://developer.r-project.org/Blog/public/2018/03/23/maximum-number-of-dlls/) 在R3.3版本中，只能有100個固定的動態庫限制，到了3.4版本以後，就能夠使用Sys.setenv(R_MAX_NUM_DLLS=xxx)進行設置，而這個數字根據個人情況設定
  • 在新建數據框時會自動將字元串的列當做是因子型向量，但是我們常常還需要對字元進行修改，因此需要先將這個設置取消：options(stringsAsFactors = F)
  • 因為Bioconductor下載方法的變動，要學會使用BiocManager::install這個命令，例如：BiocManager::install(c( 'scran'),ask = F,update = F)，新加的兩個選項表示：不要問我要不要下載，直接下！還有不要問我要不要更新，不更新！【除非不升級就報錯】
  • 下載包存在網路的限制，畢竟R語言是國外開發，因此可以通過options()$repos看看常規CRAN安裝R包的使用鏡像(一般情況下是rstudio公司的)，但是這裡我們可以自行設置：比如設置成清華源：options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))，另外Bioconductor也有自己的鏡像設置：修改一下options即可，options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") # 總結一下，可以先用if判斷再進行設置 if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F) if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
  • 如何使用if判斷語句進行包的安裝： if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

最後，就是安裝所有必備的R包(包括CRAN和Bioconductor)

# 快速安裝cran包  cran_pkgs <- c("ggfortify","FactoMineR","factoextra")  for (pkg in cran_pkgs){    if (! require(pkg,character.only=T) ) {      install.packages(pkg,ask = F,update = F)      require(pkg,character.only=T)      }  }  # Bioconductor包  library(BiocManager)  bioc_pkgs <- c("scran","TxDb.Mmusculus.UCSC.mm10.knownGene","org.Mm.eg.db","genefu","org.Hs.eg.db","TxDb.Hsapiens.UCSC.hg38.knownGene")  for (pkg in bioc_pkgs){    if (! require(pkg,character.only=T) ) {      BiocManager::install(pkg,ask = F,update = F)      require(pkg,character.only=T)      }   }

目的：利用R包重複文章的基因數量圖、聚類圖、基因在聚類圖中的熱圖、每個基因表達量在不同cluster的小提琴圖

準備工作之表達矩陣

看到文章中有兩個表達矩陣，其中第一個是原始表達矩陣(均為整數），第二個是rpkm是表達量歸一化後的值(包含了小數)，因此也能說明為何第二個文件比第一個要大。

RPKM這個指標可以這樣理解：R表示reads，K表示基因長度，M表示文庫大小，它實際上做的事情也就是去掉基因長度和測序文庫的差異對reads比對數量的影響

好，先說說為什麼要去掉文庫大小差異：以這篇文章中的圖片為例：https://sci-hub.tw/https://doi.org/10.1186/s13059-018-1466-5

比如有兩個樣本，要比較三個基因ABC的表達量，圖中越高表示比對到這個基因的reads數越多，因此在同一個樣本中可以看到C>B>A，但是不同的兩個樣本呢？

測序量不同，比大小是不公平的

舉個誇張的例子：上面?的樣本（簡稱"樣本1"）中一共比對了100萬條reads，其中C基因比對到了100條；下面?的樣本（簡稱"樣本2"）中一共比對了100條reads，其中C基因比對了10條。雖然最終的數據顯示：樣本1中C基因比樣本2的C基因比對reads數多了90條，但是考慮到實際樣本情況就是，樣本2中C基因可是佔據了總比對量的十分之一，而樣本1呢？很小很小…。因此去掉M(也就是每個樣本的測序文庫大小，以Million為單位)的影響，才是比較客觀的。

同樣的，有的基因長，有的基因短，開發RPKM的人就想：基因長的比對到的reads也會更多，因此也去掉了這個差異(除以K)

但是！這個概念目前在統計上是錯誤的，因此並不建議使用這個指標

操作表達矩陣

讀取 

# 保留頭資訊，並設置分隔符為製表符tab  a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',header = T ,sep = 't')    # 讀進來以後，簡單查看一下  a[1:6,1:4]

過濾

可以看到很多基因對應的表達量都是0

下面會用到循環，但是為了方便理解，先拿其中一行為例：

x=a[1,] #比如將第一行提取出來賦值給x  # 將x中的值與1作比較(利用了R語言的循環補齊，也就是說，它會將768個值一個一個去和1做比較，然後返回邏輯值TRUE或者FALSE)  x>1  # 然後利用table()函數檢查x中有多少是TRUE，多少是FALSE  table(x>1)  # FALSE  TRUE  #   766     2  # 可以看到第一行這個基因在768個細胞中只有兩個細胞有表達，我們認為：這兩個細胞也不好分組，cluster聚類也沒有什麼意義，因此可以去掉  # 但是這個細胞量設置成多少合適呢？總不能不能一股腦全設成2吧  floor(ncol(a)/50) # 用總列數除以50然後向下取整，結果就是15  # 也就是說，只要一行中至少要在15個樣本中有表達量  # 上面知道了 x>1 返回邏輯值0和1，0為FALSE，1為TRUE。現在我們要找一行中總共有多少TRUE，就用sum計算一下(因為會忽略掉0的影響)  sum(x>1) > floor(ncol(a)/50)  # 當然第一行會返回FALSE，也就表明我們要去掉這一行內容  a[sum(x>1) > floor(ncol(a)/50),]# 就把不符合要求的第一行去掉了

上面，我們對一行的篩選與過濾有了認識，那麼一個表達矩陣有2萬多行，怎樣實現循環操作呢？

# 專業的事情交給專業的工具去處理=》apply  # 要使用apply函數先要明白三個問題：對誰進行操作？對行還是列進行操作？操作什麼？  apply(a, 1, function(x) {sum(x>1) > floor(ncol(a)/50)})  # 1：對a這個矩陣進行操作  # 2：對行(也就是1表示)進行操作[補充：如果對列操作，用2表示]  # 3：操作什麼？複雜的操作先寫上 function(x){}，這是一個標準格式，然後大括弧中是要進行操作的函數，於是我們就可以將我們之前寫的那一行粘到這裡，最後仍然是邏輯值

最後，有多少行就會返回多少個apply判斷的邏輯值，顯示FALSE的就是要過濾掉的，於是再用行篩選完成整個操作，並賦值給一個新變數：

dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]  dim(dat)  # 12198   768 最終就保留了12198個基因

其實原文保留的更少，原文只有10835個基因

VirMach 便宜 VPS

QNews

QNews