劉小澤寫于19.6.14-15-第二單元第二、三講：獲取Github代碼包以及準(zhǔn)備工作

筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

先下載代碼包

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

首先進(jìn)入RNA-seq目錄，從step0-step9是對(duì)常規(guī)轉(zhuǎn)錄組的一個(gè)回顧

準(zhǔn)備工作之R包

從step0開始，代碼注釋蠻詳細(xì)的，我會(huì)挑選重要的部分寫到這里，其他可以自行看代碼學(xué)習(xí)，下面就是主要利用Rstudio進(jìn)行操作了

• 一個(gè)好習(xí)慣，做新項(xiàng)目時(shí)記得清空之前的變量，使用rm(list = ls())

• 有的R包比較大，經(jīng)常需要加載其他的動(dòng)態(tài)庫(kù)dynamically loaded libraries (DLLs)，例如：

  > length(loadedNamespaces()) 
  [1] 34
  > library(Seurat) #加載一個(gè)seurat包會(huì)出現(xiàn)接近60個(gè)依賴的動(dòng)態(tài)庫(kù)
  > length(loadedNamespaces())
  [1] 96
  

  如果不設(shè)置，就會(huì)因?yàn)?strong>加載數(shù)量超限制而報(bào)錯(cuò)(https://developer.r-project.org/Blog/public/2018/03/23/maximum-number-of-dlls/)

  在R3.3版本中，只能有100個(gè)固定的動(dòng)態(tài)庫(kù)限制，到了3.4版本以后，就能夠使用Sys.setenv(R_MAX_NUM_DLLS=xxx)進(jìn)行設(shè)置，而這個(gè)數(shù)字根據(jù)個(gè)人情況設(shè)定

• 在新建數(shù)據(jù)框時(shí)會(huì)自動(dòng)將字符串的列當(dāng)做是因子型向量，但是我們常常還需要對(duì)字符進(jìn)行修改，因此需要先將這個(gè)設(shè)置取消：options(stringsAsFactors = F)

• 因?yàn)锽ioconductor下載方法的變動(dòng)，要學(xué)會(huì)使用BiocManager::install這個(gè)命令，例如：BiocManager::install(c( 'scran'),ask = F,update = F)，新加的兩個(gè)選項(xiàng)表示：不要問(wèn)我要不要下載，直接下！還有不要問(wèn)我要不要更新，不更新！【除非不升級(jí)就報(bào)錯(cuò)】

• 下載包存在網(wǎng)絡(luò)的限制，畢竟R語(yǔ)言是國(guó)外開發(fā)，因此可以通過(guò)options()$repos看看常規(guī)CRAN安裝R包的使用鏡像(一般情況下是rstudio公司的)，但是這里我們可以自行設(shè)置：比如設(shè)置成清華源：options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))，另外Bioconductor也有自己的鏡像設(shè)置：修改一下options即可，options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

  # 總結(jié)一下，可以先用if判斷再進(jìn)行設(shè)置
  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判斷語(yǔ)句進(jìn)行包的安裝：

  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包重復(fù)文章的基因數(shù)量圖、聚類圖、基因在聚類圖中的熱圖、每個(gè)基因表達(dá)量在不同cluster的小提琴圖

準(zhǔn)備工作之表達(dá)矩陣

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

RPKM這個(gè)指標(biāo)可以這樣理解：R表示reads，K表示基因長(zhǎng)度，M表示文庫(kù)大小，它實(shí)際上做的事情也就是去掉基因長(zhǎng)度和測(cè)序文庫(kù)的差異對(duì)reads比對(duì)數(shù)量的影響

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

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

測(cè)序量不同，比大小是不公平的

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

同樣的，有的基因長(zhǎng)，有的基因短，開發(fā)RPKM的人就想：基因長(zhǎng)的比對(duì)到的reads也會(huì)更多，因此也去掉了這個(gè)差異(除以K)

但是！這個(gè)概念目前在統(tǒng)計(jì)上是錯(cuò)誤的，因此并不建議使用這個(gè)指標(biāo)

操作表達(dá)矩陣

讀取

# 保留頭信息，并設(shè)置分隔符為制表符tab
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',header = T ,sep = '\t')

# 讀進(jìn)來(lái)以后，簡(jiǎn)單查看一下
a[1:6,1:4]

過(guò)濾

可以看到很多基因?qū)?yīng)的表達(dá)量都是0

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

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

上面，我們對(duì)一行的篩選與過(guò)濾有了認(rèn)識(shí)，那么一個(gè)表達(dá)矩陣有2萬(wàn)多行，怎樣實(shí)現(xiàn)循環(huán)操作呢？

# 專業(yè)的事情交給專業(yè)的工具去處理=》apply
# 要使用apply函數(shù)先要明白三個(gè)問(wèn)題：對(duì)誰(shuí)進(jìn)行操作？對(duì)行還是列進(jìn)行操作？操作什么？
apply(a, 1, function(x) {sum(x>1) > floor(ncol(a)/50)})
# 1：對(duì)a這個(gè)矩陣進(jìn)行操作
# 2：對(duì)行(也就是1表示)進(jìn)行操作[補(bǔ)充：如果對(duì)列操作，用2表示]
# 3：操作什么？復(fù)雜的操作先寫上 function(x){}，這是一個(gè)標(biāo)準(zhǔn)格式，然后大括號(hào)中是要進(jìn)行操作的函數(shù)，于是我們就可以將我們之前寫的那一行粘到這里，最后仍然是邏輯值

最后，有多少行就會(huì)返回多少個(gè)apply判斷的邏輯值，顯示FALSE的就是要過(guò)濾掉的，于是再用行篩選完成整個(gè)操作，并賦值給一個(gè)新變量：

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

其實(shí)原文保留的更少，原文只有10835個(gè)基因