劉小澤寫于19.6.14-15-第二單元第二蒜危、三講：獲取Github代碼包以及準備工作
筆記目的：根據(jù)生信技能樹的單細胞轉(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
首先進入RNA-seq目錄辐赞，從step0-step9是對常規(guī)轉(zhuǎn)錄組的一個回顧
準備工作之R包
從step0開始响委，代碼注釋蠻詳細的晃酒，我會挑選重要的部分寫到這里贝次，其他可以自行看代碼學習彰导，下面就是主要利用Rstudio進行操作了
一個好習慣位谋，做新項目時記得清空之前的變量，使用
rm(list = ls())-
有的R包比較大笋轨，經(jīng)常需要加載其他的動態(tài)庫dynamically loaded libraries (DLLs)爵政，例如：
> length(loadedNamespaces()) [1] 34 > library(Seurat) #加載一個seurat包會出現(xiàn)接近60個依賴的動態(tài)庫 > length(loadedNamespaces()) [1] 96如果不設置，就會因為加載數(shù)量超限制而報錯(https://developer.r-project.org/Blog/public/2018/03/23/maximum-number-of-dlls/)
在R3.3版本中洁灵，只能有100個固定的動態(tài)庫限制徽千，到了3.4版本以后汤锨，就能夠使用
Sys.setenv(R_MAX_NUM_DLLS=xxx)進行設置，而這個數(shù)字根據(jù)個人情況設定 在新建數(shù)據(jù)框時會自動將字符串的列當做是因子型向量荠诬，但是我們常常還需要對字符進行修改柑贞，因此需要先將這個設置取消：
options(stringsAsFactors = F)因為Bioconductor下載方法的變動聂抢，要學會使用
BiocManager::install這個命令琳疏，例如：BiocManager::install(c( 'scran'),ask = F,update = F)空盼，新加的兩個選項表示：不要問我要不要下載，直接下台汇！還有不要問我要不要更新篱瞎，不更新！【除非不升級就報錯】-
下載包存在網(wǎng)絡的限制俐筋，畢竟R語言是國外開發(fā)牵素，因此可以通過
options()$repos看看常規(guī)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/")# 總結(jié)一下请琳，可以先用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包重復文章的基因數(shù)量圖腰奋、聚類圖单起、基因在聚類圖中的熱圖、每個基因表達量在不同cluster的小提琴圖
準備工作之表達矩陣
看到文章中有兩個表達矩陣嘀倒，其中第一個是原始表達矩陣(均為整數(shù)），第二個是rpkm是表達量歸一化后的值(包含了小數(shù))局冰，因此也能說明為何第二個文件比第一個要大测蘑。
RPKM這個指標可以這樣理解：R表示reads，K表示基因長度康二，M表示文庫大小碳胳，它實際上做的事情也就是去掉基因長度和測序文庫的差異對reads比對數(shù)量的影響
好，先說說為什么要去掉文庫大小差異：以這篇文章中的圖片為例：https://sci-hub.tw/https://doi.org/10.1186/s13059-018-1466-5
比如有兩個樣本沫勿，要比較三個基因ABC的表達量挨约，圖中越高表示比對到這個基因的reads數(shù)越多，因此在同一個樣本中可以看到C>B>A产雹，但是不同的兩個樣本呢诫惭？
測序量不同，比大小是不公平的
舉個夸張的例子：上面??的樣本（簡稱"樣本1"）中一共比對了100萬條reads蔓挖，其中C基因比對到了100條共缕；下面??的樣本（簡稱"樣本2"）中一共比對了100條reads惕橙，其中C基因比對了10條。雖然最終的數(shù)據(jù)顯示：樣本1中C基因比樣本2的C基因比對reads數(shù)多了90條汹来，但是考慮到實際樣本情況就是利耍，樣本2中C基因可是占據(jù)了總比對量的十分之一睬愤，而樣本1呢待榔？很小很小…屁使。因此去掉M(也就是每個樣本的測序文庫大小，以Million為單位)的影響匆瓜，才是比較客觀的咽扇。
同樣的，有的基因長陕壹，有的基因短，開發(fā)RPKM的人就想：基因長的比對到的reads也會更多树埠，因此也去掉了這個差異(除以K)
但是糠馆！這個概念目前在統(tǒng)計上是錯誤的，因此并不建議使用這個指標
操作表達矩陣
讀取
# 保留頭信息怎憋，并設置分隔符為制表符tab
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',header = T ,sep = '\t')
# 讀進來以后又碌，簡單查看一下
a[1:6,1:4]
過濾
可以看到很多基因?qū)谋磉_量都是0
下面會用到循環(huán)九昧，但是為了方便理解，先拿其中一行為例：
x=a[1,] #比如將第一行提取出來賦值給x
# 將x中的值與1作比較(利用了R語言的循環(huán)補齊毕匀，也就是說铸鹰，它會將768個值一個一個去和1做比較，然后返回邏輯值TRUE或者FALSE)
x>1
# 然后利用table()函數(shù)檢查x中有多少是TRUE皂岔，多少是FALSE
table(x>1)
# FALSE TRUE
# 766 2
# 可以看到第一行這個基因在768個細胞中只有兩個細胞有表達蹋笼，我們認為：這兩個細胞也不好分組，cluster聚類也沒有什么意義躁垛，因此可以去掉
# 但是這個細胞量設置成多少合適呢剖毯？總不能不能一股腦全設成2吧
floor(ncol(a)/50) # 用總列數(shù)除以50然后向下取整，結(jié)果就是15
# 也就是說教馆，只要一行中至少要在15個樣本中有表達量
# 上面知道了 x>1 返回邏輯值0和1逊谋，0為FALSE，1為TRUE⊥疗蹋現(xiàn)在我們要找一行中總共有多少TRUE胶滋，就用sum計算一下(因為會忽略掉0的影響)
sum(x>1) > floor(ncol(a)/50)
# 當然第一行會返回FALSE，也就表明我們要去掉這一行內(nèi)容
a[sum(x>1) > floor(ncol(a)/50),]# 就把不符合要求的第一行去掉了
上面悲敷，我們對一行的篩選與過濾有了認識究恤，那么一個表達矩陣有2萬多行，怎樣實現(xiàn)循環(huán)操作呢镀迂？
# 專業(yè)的事情交給專業(yè)的工具去處理=》apply
# 要使用apply函數(shù)先要明白三個問題：對誰進行操作丁溅？對行還是列進行操作？操作什么探遵？
apply(a, 1, function(x) {sum(x>1) > floor(ncol(a)/50)})
# 1：對a這個矩陣進行操作
# 2：對行(也就是1表示)進行操作[補充：如果對列操作窟赏，用2表示]
# 3：操作什么？復雜的操作先寫上 function(x){}箱季，這是一個標準格式涯穷，然后大括號中是要進行操作的函數(shù)，于是我們就可以將我們之前寫的那一行粘到這里藏雏，最后仍然是邏輯值
最后拷况，有多少行就會返回多少個apply判斷的邏輯值，顯示FALSE的就是要過濾掉的掘殴，于是再用行篩選完成整個操作赚瘦，并賦值給一個新變量：
dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
dim(dat)
# 12198 768 最終就保留了12198個基因
其實原文保留的更少，原文只有10835個基因
