理論基礎：線性模型,， 設計矩陣和比較矩陣
這部分內容最先在 RNA-Seq Data Analysis 的8.5.3節(jié)看到，剛開始一點都不理解,，但是學完生物統(tǒng)計之后,，我認為這是理解所有差異基因表達分析R包的關鍵。
基本上,，統(tǒng)計課都會介紹如何使用t檢驗用來比較兩個樣本之間的差異，然后在樣本比較多的時候使用方差分析確定樣本間是否有差異,。當然前是樣本來自于正態(tài)分布的群體,，或者隨機獨立大量抽樣。
對于基因芯片的差異表達分析而言,，由于普遍認為其數據是服從正態(tài)分布,，因此差異表達分析無非就是用t檢驗和或者方差分析應用到每一個基因上。高通量一次性找的基因多,，于是就需要對多重試驗進行矯正,，控制假陽性,。目前在基因芯片的分析用的最多的就是limma,。
但是,，高通量測序(HTS)的read count普遍認為是服從泊松分布（當然有其他不同意見）,，不可能直接用正態(tài)分布的t檢驗和方差分析。 當然我們可以簡單粗暴的使用對于的非參數檢驗的方法,，但是統(tǒng)計力不夠,，結果的p值矯正之估計一個差異基因都找不到,。老板花了一大筆錢，結果卻說沒有差異基因,，是個負結果,，于是好幾千經費打了水漂，他肯定是不樂意的,。因此,，還是得要用參數檢驗的方法,，于是就要說到方差分析和線性模型之間的關系了,。
線性回歸和方差分析是同一時期發(fā)展出的兩套方法,。在我本科階段的田間統(tǒng)計學課程中就介紹用方差分析（ANOVA）分析不同肥料處理后的產量差異,，實驗設計如下
| 肥料 |
重復1 |
重復2 |
重復3 |
重復4 |
| A1 |
... |
... |
... |
... |
| A2 |
... |
... |
... |
... |
| A3 |
... |
... ... |
... |
... |
這是最簡單的單因素方差分析,，每一個結果都可以看成 yij = ai + u + eij,， 其中u是總體均值，ai是每一個處理的差異,，eij是隨機誤差,。
注：方差分析(Analysis of Variance, ANAOVA)名字聽起來好像是檢驗方差,，但其實是為了判斷樣本之間的差異是否真實存在，為此需要證明不同處理內的方差顯著性大于不同處理間的方差,。
線性回歸 一般是用于量化的預測變量來預測量化的響應變量,。比如說體重與身高的關系建模：
當然線性回歸也可用處理名義型或有序型因子（也就是離散變量）作為預測變量，如果要畫圖的話,，就是下面這個情況,。
如果我們需要通過一個實驗找到不同處理后對照組和控制組的基因變化，那么基因表達可以簡單寫成,， y = a + b · treament + e,。 和之前的 yij = ai + u + eij 相比，你會發(fā)現公式是如此的一致,。 這是因為線性模型和方差分析都是廣義線性模型(generalizing linear models, GLM)在正態(tài)分布的預測變量的特殊形式,。而GLM本身只要采用合適的連接函數是可以處理對任意類型的變量進行建模的。
目前認為read count之間的差異是符合負二項分布,，也叫gamma-Possion分布,。那么問題來了，如何用GLM或者LM分析兩個處理件的差異呢,？其實可以簡單的用上圖的擬合直線的斜率來解釋,，如果不同處理之間存在差異，那么這個擬合線的斜率必定不為零,，也就是與X軸平行,。但是這是一種便于理解的方式（雖然你也未必能理解），實際更加復雜,，考慮因素更多,。
注1 負二向分布有兩個參數，均值(mean)和離散值（dispersion). 離散值描述方差偏離均值的程度,。泊松分布可以認為是負二向分布的離散值為1,，也就是均值等于方差（mean=variance）的情況。
注2 這部分涉及大量的統(tǒng)計學知識,，不懂就用維基百科一個個查清楚,。
聊完了線性模型和方差分析，下面的設計矩陣（design matrix）就很好理解了,， 其實就是用來告訴不同的差異分析函數應該如何對待變量,。比如說我們要研究的KD和control之間變化，設計矩陣就是
| 樣本 |
處理 |
| sample1 |
control |
| sample2 |
control |
| sample3 |
KD |
| sample4 |
KD |
那么比較矩陣（contrast matrix)就是告訴差異分析函數應該如何對哪個因素進行比較,， 這里就是比較不同處理下表達量的變化。
標準化一二事
其實read count如何標準化的方法有很多,，最常用的是FPKM和RPKM,，雖然它們其實是錯的--FPKM/RPKM是錯的,。
我推薦閱讀 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同標準化方法之間的差異。
有一些方法是要求原始數據,，有一些則要求經過某類標準化后的數據,，記得區(qū)分。
探索性分析一二事
使用DESeq2進行差異基因分析
關于DESeq2分析差異表達基因,，其實在https://www./help/workflows/rnaseqGene/ 里面介紹的非常清楚了,。
我們已經準備好了count matrix，接下來就是把數據導入DESeq2,。DESeq2導入數據的方式有如下4種,，基本覆蓋了主流read count軟件的結果。
注 DESeq2要求的數據是raw count,， 沒必要進行FPKM/TPM/RPFKM/TMM標準化,。
| function |
package |
framework |
output |
DESeq2 input function |
| summarizeOverlaps |
GenomicAlignments |
R/Bioconductor |
SummarizedExperiment |
DESeqDataSet |
| featureCounts |
Rsubread |
R/Bioconductor |
matrix |
DESeqDataSetFromMatrix |
| tximport |
tximport |
R/Bioconductor |
list of matrices |
DESeqDataSetFromTximport |
| htseq-count |
HTSeq |
Python |
files |
DESeqDataSetFromHTSeq |
本來我們是可以用DESeq2為htseq-count專門提供的 DESeqDataSetFromHTSeq ，然而很尷尬數據不夠要自己湊數,，所以只能改用 DESeqDataSetFromMatrix了 :cold_sweat:
導入數據,，構建 DESeq2 所需的 DESeqDataSet 對象
|
|
library(DESeq2) countData <- raw_count_filt[,2:5] condition <- factor(c("control","KD","KD","control")) dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition ) |
注： 這一步到下一步之間可以過濾掉一些low count數據，節(jié)省內存,，提高運行速度
|
|
nrow(dds) dds <- dds[ rowSums(counts(dds)) > 1, ] nrow(dds) |
使用DESeq進行差異表達分析： DESeq包含三步,，estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons),， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest）,，可以分布運行，也可用一步到位,，最后返回 results可用的DESeqDataSet對象,。
|
|
dds <- DESeq(dds) # 出現如下提示信息，說明運行成功 estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing |
用results獲取結果： results的參數非常的多,，這里不好具體展開 :pensive: 但是你們會自己看的吧
我們可用mcols查看每一項結果的具體含義,，比如說 log2FoldChange 表示倍數變化取log2結果，還能畫個火山圖,。一般簡單粗暴的用2到3倍作為閾值,，但是對于低表達的基因，3倍也是噪音,，那些高表達的基因,，1.1倍都是生物學顯著了。更重要的沒有考慮到組內變異,，沒有統(tǒng)計學意義,。 padj 就是用BH對多重試驗進行矯正。
|
|
mcols(res, use.names = TRUE) DataFrame with 6 rows and 2 columns type description <character> <character> baseMean intermediate mean of normalized counts for all samples log2FoldChange results log2 fold change (MLE): condition KD vs control lfcSE results standard error: condition KD vs control stat results Wald statistic: condition KD vs control pvalue results Wald test p-value: condition KD vs control padj results BH adjusted p-values |
用summary看描述性的結果，大致是上調的基因占總體的11%,，下調的是7.1%（KD vs control）
|
|
summary(res) out of 29469 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 3154, 11% LFC < 0 (down) : 2095, 7.1% outliers [1] : 0, 0% low counts [2] : 15111, 51% (mean count < 22) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results |
畫個MA圖,，還能標注p值最小的基因。
An MA plot is an application of a Bland–Altman plot for visual representation of genomic data. The plot visualises the differences between measurements taken in two samples, by transforming the data onto M (log ratio) and A (mean average) scales, then plotting these values. Though originally applied in the context of two channel DNA microarray gene expression data, MA plots are also used to visualise high-throughput sequencing analysis --From wikipeida
M表示log fold change,，衡量基因表達量變化,，上調還是下調。A表示每個基因的count的均值,。根據summary可知,，low count的比率很高，所以大部分基因表達量不高,，也就是集中在0的附近（log2(1)=0,，也就是變化1倍）.提供了模型預測系數的分布總覽。
下圖是沒有經過 statistical moderation平緩log2 fold changes的情況
|
|
plotMA(res, ylim = c(-5,5)) topGene <- rownames(res)[which.min(res$padj)] with(res[topGene, ], { points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2) text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue") }) |
如果經過lfcShrink 收縮log2 fold change,， 結果會好看很多
|
|
res.shrink <- lfcShrink(dds, contrast = c("condition","KD","control"), res=res) plotMA(res.shrink, ylim = c(-5,5)) topGene <- rownames(res)[which.min(res$padj)] with(res[topGene, ], { points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2) text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue") }) |
當然還有火山圖,，不過留給其他方法作圖，我們先把差異表達的基因找出來,。
|
|
res.deseq2 <- subset(res, padj < 0.05) |
一般p value 小于0.05就是顯著了, 顯著性不代表結果正確,，只用于給后續(xù)的富集分析和GSEA提供排序標準和篩選而已。關于P值的吐槽簡直無數,， 請多注意,。
使用edgeR進行差異基因分析
edgeR在函數說明中稱其不但可以分析SAGE， CAGE的RNA-Seq,，Tag-RNA,，或RNA-seq， 也能分析ChIP-Seq和CRISPR得到的read counts數據,。嗯,，我信了:confused:！
edgeR使用DGEList函數讀取count matrix數據,，也就說你需要提供一個現成的matrix數據,，而不是指望它能讀取單獨的文件，然后進行合并(當然機智的我發(fā)現,，其實可以用 tximport 或 DESeqDataSetFromHTSeq 讀取單獨的文件,，然后傳遞給DGEList)
第一步： 構建DGEList對象
|
|
library(edgeR) group <- factor(c("control","KD","KD","control")) genelist <- DGEList(counts=raw_count_filt[,2:5], group = group) |
第二步： 過濾 low counts數據。與DESeq2的預過濾不同,，DESeq2的預過濾只是為了改善后續(xù)運算性能,，在運行過程中依舊會自動處理low count數據，edgeR需要在分析前就要排除那些low count數據,，而且非常嚴格,。從生物學角度，有生物學意義的基因的表達量必須高于某一個閾值。從統(tǒng)計學角度上,， low count的數據不太可能有顯著性差異,，而且在多重試驗矯正階段還會拖后腿。 綜上所訴,，放心大膽的過濾吧。
根據經驗（又是經驗 :dog: ）,， 基因至少在某一些文庫的count超過10 ~ 15 才被認為是表達,。這一步全靠嘗試， 剔除太多就緩緩,，剔除太少就嚴格點,。 我們可以簡單的對每個基因的raw count進行比較，但是建議用CPM（count-per-million)標準化 后再比較,，避免了文庫大小的影響,。
|
|
# 簡單粗暴的方法 keep <- rowSums(genelist$count) > 50 # 利用CPM標準化 keep <- rowSums(cpm(genelist) > 0.5 ) >=2 table(keep) genelist.filted <- genelist[keep, ,keep.lib.sizes=FALSE] |
這里的0.5(即閾值）等于 10/(最小的文庫的 read count數 /1000000)，keep.lib.size=FALSE表示重新計算文庫大小,。
第三步： 根據組成偏好(composition bias)標準化,。edgeR的calcNormFactors函數使用TMM算法對DGEList標準化
|
|
genelist.norm <- calcNormFactors(genelist.filted) |
注 大部分的mRNA-Seq數據分析用TMM標準化就行了，但是也有例外,，比如說single-cell RNA-Seq(Lun, Bach, and Marioni 2016), 還有就是global differential expression,， 基因組一半以上的基因都是差異表達的，請盡力避免,，(D. Wu et al. 2013),， 不然就需要用到內參進行標準化了(Risso et al. 2014).
第四步： 實驗設計矩陣(Design matrix)， 類似于DESeq2中的design參數,。 edgeR的線性模型和差異表達分析需要定義一個實驗設計矩陣,。很直白的就能發(fā)現是1vs0
|
|
design <- model.matrix(~0+group) colnames(design) <- levels(group) design control KD 1 1 0 2 0 1 3 0 1 4 1 0 |
第五步： 估計離散值（Dispersion）。前面已經提到負二項分布（negative binomial,，NB)需要均值和離散值兩個參數,。edgeR對每個基因都估測一個經驗貝葉斯穩(wěn)健離散值（mpirical Bayes moderated dispersion），還有一個公共離散值（common dispersion,，所有基因的經驗貝葉斯穩(wěn)健離散值的均值）以及一個趨勢離散值
|
|
genelist.Disp <- estimateDisp(genelist.norm, design, robust = TRUE) plotBCV(genelist.Disp) |
還可以進一步通過quasi-likelihood (QL)擬合NB模型,，用于解釋生物學和技術性導致的基因特異性變異 (Lund et al. 2012; Lun, Chen, and Smyth 2016).
|
|
fit <- glmQLFit(genelist.Disp, design, robust=TRUE) head(fit$coefficients) |
注1 估計離散值這個步驟其實有許多estimate*Disp函數。當不存在實驗設計矩陣(design matrix)的時候,，estimateDisp 等價于 estimateCommonDisp 和 estimateTagwiseDisp ,。而當給定實驗設計矩陣(design matrix)時， estimateDisp 等價于 estimateGLMCommonDisp, estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp,。 其中tag與gene同義,。
注2 其實這里的第三， 四， 五步對應的就是DESeq2的DESeq包含的2步,，標準化和離散值估測,。
第六步: 差異表達檢驗（1）。這一步主要構建比較矩陣,，類似于DESeq2中的results函數的 contrast 參數,。
|
|
cntr.vs.KD <- makeContrasts(control-KD, levels=design) res <- glmQLFTest(fit, contrast=cntr.vs.KD) ig.edger <- res$table[p.adjust(res$table$PValue, method = "BH") < 0.01, ] |
這里用的是glmQLFTest而不是glmLRT是因為前面用了glmQLTFit進行擬合，所以需要用QL F-test進行檢驗,。如果前面用的是glmFit,，那么對應的就是glmLRT. 作者稱QL F-test更加嚴格。多重試驗矯正用的也是BH方法,。
后續(xù)就是提取顯著性差異的基因用作下游分析,，做一些圖看看
|
|
topTags(res,n=10) is.de <- decideTestsDGE(res) summary(is.de) plotMD(res, status=is.de, values=c(1,-1), col=c("red","blue"), legend="topright") |
第六步：差異表達檢驗（2）。上面找到的顯著性差異的基因,，沒有考慮效應值,，也就是具體變化了多少倍。我們也可用找表達量變化比較大的基因,，對應的函數是 glmTreat,。
|
|
tr <- glmTreat(fit, contrast=B.LvsP, lfc=log2(1.5)) s |
使用limma進行差異分析
經過上面兩個方法的洗禮，基本上套路你也就知道了,，我先簡單小結一下,，然后繼續(xù)介紹limma包的 voom 。
- 導入read count,， 保存為專門的對象用于后續(xù)分析
- 原始數據過濾,，根據標準化read count 或者 raw count 作為篩選標準
- raw read count 標準化
- 通過各種算法（如經驗貝葉斯，EM）預測dispersion離散值
- 廣義線性模型擬合數據
- 差異分析,，也就是統(tǒng)計檢驗部分
Limma原先用于處理基因表達芯片數據,，可是說是這個領域的老大 :sunglasses: 。如果你仔細看edgeR導入界面,，你就會發(fā)現,，edgeR有一部分功能依賴于limma包。Limma采用經驗貝葉斯模型（ Empirical Bayesian model）讓結果更穩(wěn)健,。
在處理RNA-Seq數據時,，raw read count先被轉成log2-counts-per-million (logCPM)，然后對mean-variance關系建模,。建模有兩種方法：
- 精確權重法（precision weights）也就是“voom"
- 經驗貝葉斯先驗趨勢（empirical Bayes prior trend）,，也就是”limma-trend“
數據預處理： Limma使用edgeR的DGEList對象，并且過濾方法都是一致的,，對應edgeR的第一步,第二步,， 第三步
|
|
library(edgeR) library(limma) group <- factor(c("control","KD","KD","control")) genelist <- DGEList(counts=raw_count_filt[,2:5], group = group) ### filter base use CPM keep <- rowSums(cpm(genelist) > 0.5 ) >=2 table(keep) genelist.filted <- genelist[keep, ,keep.lib.sizes=FALSE] ### normalizaition x <- calcNormFactors(x, method = "TMM") |
差異表達分析: 使用”limma-trend“
|
|
design <- model.matrix(~0+group) colnames(design) <- levels(group) logCPM <- cpm(genelist.norm, log=TRUE, prior.count=3) fit <- lmFit(logCPM, design) fit <- eBayes(fit, trend=TRUE) topTable(fit, coef=ncol(design)) |
差異表達分析: 使用”limma-voom“
|
|
### DGE with voom v <- voom(genelist.norm, design, plot=TRUE) #v <- voom(counts, design, plot=TRUE) fit <- lmFit(v, design) fit <- eBayes(fit) all <- topTable(fit, coef=ncol(design), number=10000) sig.limma <- all[all$adj.P.Val < 0.01, ] fit <- treat(fit, lfc=log2(1.2)) topTreat(fit, coef=ncol(design)) |
如果分析基因芯片數據,，必須好好讀懂LIMMA包。
不同軟件包分析結果比較
基本上每一個包,，我都提取了各種的顯著性基因,，比較就需要用韋恩圖了，但是我偏不 :stuck_out_tongue: 我要用UpSetR.
|
|
library(UpSetR) input <- fromList(list(edgeR=rownames(sig.edger), DESeq2=rownames(sig.deseq2), limma=rownames(sig.limma))) |
感覺limma的結果有點奇怪,，有生之年在折騰吧,。
使用GFOLD進行無重復樣本的差異基因分析
|
