㈠ 轉錄組入門（7）：差異表達分析
原先三個樣本的HTSeq-count計數的數據可以在我的GitHub中找到，但是前面已經說過Jimmy失誤讓我們分析的人類就只有3個樣本， 另外一個樣本需要從另一批數據獲取（請注意batch effect)，所以不能保證每一組都有兩個重復。
我一直堅信」你並不孤獨「這幾個字，遇到這種情況的人肯定不止我一個，於是我找到了幾種解決方法
以上方法都會在後續進行介紹，但是我們DESeq2必須得要有重復的問題亟待解決，沒辦法我只能自己瞎編了。雖然是編，我們也要有模有樣，不能直接復制一份，要考慮到高通量測序的read是默認符合泊松分布的。我是這樣編的。
這僅僅是一種填坑的方法而已，更好模擬數據的方法需要參閱更加專業的文獻， 有生之年 我希望能補上這一個部分。
這部分內容最先在 RNA-Seq Data Analysis 的8.5.3節看到，剛開始一點都不理解，但是學完生物統計之後，我認為這是理解所有差異基因表達分析R包的關鍵。
基本上，統計課都會介紹如何使用 t檢驗 用來比較兩個樣本之間的差異，然後在樣本比較多的時候使用 方差分析 確定樣本間是否有差異。當然前是樣本來自於正態分布的群體，或者隨機獨立大量抽樣。
對於基因晶元的差異表達分析而言，由於普遍認為其數據是服從正態分布，因此差異表達分析無非就是用t檢驗和或者方差分析應用到每一個基因上。高通量一次性找的基因多，於是就需要對多重試驗進行矯正，控制假陽性。目前在基因晶元的分析用的最多的就是 limma 。
但是 ，高通量測序(HTS)的read count普遍認為是服從泊松分布（當然有其他不同意見），不可能直接用正態分布的 t檢驗 和 方差分析 。 當然我們可以簡單粗暴的使用對於的 非參數檢驗 的方法，但是統計力不夠，結果的p值矯正之估計一個差異基因都找不到。老闆花了一大筆錢，結果卻說沒有差異基因，是個負結果，於是好幾千經費打了水漂，他肯定是不樂意的。因此，還是得要用參數檢驗的方法，於是就要說到方差分析和線性模型之間的關系了。
線性回歸和方差分析是同一時期發展出的兩套方法。在我本科階段的田間統計學課程中就介紹用 方差分析 （ANOVA）分析不同肥料處理後的產量差異，實驗設計如下
這是最簡單的單因素方差分析，每一個結果都可以看成 yij = ai + u + eij， 其中u是總體均值，ai是每一個處理的差異，eij是隨機誤差。
注 ：方差分析(Analysis of Variance, ANAOVA)名字聽起來好像是檢驗方差，但其實是為了判斷樣本之間的差異是否真實存在，為此需要證明不同處理內的方差顯著性大於不同處理間的方差。
線性回歸 一般是用於量化的預測變數來預測量化的響應變數。比如說體重與身高的關系建模：
當然線性回歸也可用處理名義型或有序型因子（也就是離散變數）作為預測變數，如果要畫圖的話，就是下面這個情況。
如果我們需要通過一個實驗找到不同處理後對照組和控制組的基因變化，那麼基因表達可以簡單寫成， y = a + b · treament + e。 和之前的 yij = ai + u + eij 相比，你會發現公式是如此的一致。 這是因為線性模型和方差分析都是 廣義線性模型 (generalizing linear models, GLM)在正態分布的預測變數的特殊形式。而GLM本身只要採用合適的 連接函數 是可以處理對任意類型的變數進行建模的。
目前認為read count之間的差異是符合負二項分布，也叫gamma-Possion分布。那麼問題來了，如何用GLM或者LM分析兩個處理件的差異呢？其實可以簡單的用上圖的擬合直線的斜率來解釋，如果不同處理之間存在差異，那麼這個擬合線的斜率必定不為零，也就是與X軸平行。但是這是一種便於理解的方式（雖然你也未必能理解），實際更加復雜，考慮因素更多。
注1 負二向分布有兩個參數，均值(mean)和離散值（dispersion). 離散值描述方差偏離均值的程度。泊松分布可以認為是負二向分布的離散值為1，也就是均值等於方差（mean=variance）的情況。
注2 這部分涉及大量的統計學知識，不懂就用維基網路一個個查清楚。
聊完了線性模型和方差分析，下面的設計矩陣（design matrix）就很好理解了， 其實就是用來告訴不同的差異分析函數應該如何對待變數。比如說我們要研究的KD和control之間變化，設計矩陣就是
那麼比較矩陣（contrast matrix)就是告訴差異分析函數應該如何對哪個因素進行比較， 這里就是比較不同處理下表達量的變化。
其實read count如何標准化的方法有很多，最常用的是FPKM和RPKM，雖然它們其實是錯的-- FPKM/RPKM是錯的 。
我推薦閱讀 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同標准化方法之間的差異。
有一些方法是要求原始數據，有一些則要求經過某類標准化後的數據，記得區分。
關於DESeq2分析差異表達基因，其實在 https://www.bioconctor.org/help/workflows/rnaseqGene/ 裡面介紹的非常清楚了。
我們已經准備好了count matrix，接下來就是把數據導入DESeq2。DESeq2導入數據的方式有如下4種，基本覆蓋了主流read count軟體的結果。
注 DESeq2要求的數據是raw count， 沒必要進行FPKM/TPM/RPFKM/TMM標准化。
本來我們是可以用DESeq2為htseq-count專門提供的 DESeqDataSetFromHTSeq ，然而很尷尬數據不夠要自己湊數，所以只能改用 DESeqDataSetFromMatrix 了 :cold_sweat:
導入數據，構建 DESeq2 所需的 DESeqDataSet 對象
注 ： 這一步到下一步之間可以過濾掉一些low count數據，節省內存，提高運行速度
使用 DESeq 進行差異表達分析： DESeq 包含三步，estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons)， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest），可以分布運行，也可用一步到位，最後返回 results 可用的DESeqDataSet對象。
用results獲取結果： results的參數非常的多，這里不好具體展開 :pensive: 但是你們會自己看的吧
我們可用mcols查看每一項結果的具體含義，比如說 log2FoldChange 表示倍數變化取log2結果，還能畫個火山圖。一般簡單粗暴的用2到3倍作為閾值，但是對於低表達的基因，3倍也是噪音，那些高表達的基因，1.1倍都是生物學顯著了。更重要的沒有考慮到組內變異，沒有統計學意義。 padj 就是用BH對多重試驗進行矯正。
用summary看描述性的結果，大致是上調的基因占總體的11%，下調的是7.1%（KD vs control）
畫個MA圖，還能標注p值最小的基因。
下圖是沒有經過 statistical moderation平緩log2 fold changes的情況
如果經過 lfcShrink 收縮log2 fold change， 結果會好看很多
當然還有火山圖，不過留給其他方法作圖，我們先把差異表達的基因找出來。
一般p value 小於0.05就是顯著了, 顯著性不代表結果正確，只用於給後續的富集分析和GSEA提供排序標准和篩選而已。關於P值的吐槽簡直無數， 請多注意。
edgeR在函數說明中稱其不但可以分析SAGE， CAGE的RNA-Seq，Tag-RNA，或RNA-seq， 也能分析ChIP-Seq和CRISPR得到的read counts數據。嗯，我信了:confused:！
edgeR使用 DGEList 函數讀取count matrix數據，也就說你需要提供一個現成的matrix數據，而不是指望它能讀取單獨的文件，然後進行合並(當然機智的我發現，其實可以用 tximport 或 DESeqDataSetFromHTSeq 讀取單獨的文件，然後傳遞給 DGEList )
第一步： 構建DGEList對象
第二步： 過濾 low counts數據。與DESeq2的預過濾不同，DESeq2的預過濾只是為了改善後續運算性能，在運行過程中依舊會自動處理low count數據，edgeR需要在分析前就要排除那些low count數據，而且非常嚴格。從生物學角度，有生物學意義的基因的表達量必須高於某一個閾值。從統計學角度上， low count的數據不太可能有顯著性差異，而且在多重試驗矯正階段還會拖後腿。 綜上所訴，放心大膽的過濾吧。
根據經驗（又是經驗 :dog: ）， 基因至少在某一些文庫的count超過10 ~ 15 才被認為是表達。這一步全靠嘗試， 剔除太多就緩緩，剔除太少就嚴格點。 我們可以簡單的對每個基因的raw count進行比較，但是建議用CPM（count-per-million) 標准化 後再比較，避免了 文庫大小 的影響。
這里的0.5(即閾值）等於 10/(最小的文庫的 read count數 /1000000)，keep.lib.size=FALSE表示重新計算文庫大小。
第三步： 根據組成偏好(composition bias)標准化。edgeR的 calcNormFactors 函數使用 TMM演算法 對DGEList標准化
注 大部分的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的線性模型和差異表達分析需要定義一個實驗設計矩陣。很直白的就能發現是1vs0
第五步： 估計離散值（Dispersion）。前面已經提到負二項分布（negative binomial，NB)需要均值和離散值兩個參數。edgeR對每個基因都估測一個經驗貝葉斯穩健離散值（mpirical Bayes moderated dispersion），還有一個公共離散值（common dispersion，所有基因的經驗貝葉斯穩健離散值的均值）以及一個趨勢離散值
還可以進一步通過quasi-likelihood (QL)擬合NB模型，用於解釋生物學和技術性導致的基因特異性變異 (Lund et al. 2012; Lun, Chen, and Smyth 2016).
注1 估計離散值這個步驟其實有許多 estimate*Disp 函數。當不存在實驗設計矩陣(design matrix)的時候， estimateDisp 等價於 estimateCommonDisp 和 estimateTagwiseDisp 。而當給定實驗設計矩陣(design matrix)時， estimateDisp 等價於 estimateGLMCommonDisp , estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp 。 其中tag與gene同義。
注2 其實這里的第三， 四， 五步對應的就是DESeq2的 DESeq 包含的2步，標准化和離散值估測。
第六步: 差異表達檢驗（1）。這一步主要構建比較矩陣，類似於DESeq2中的 results 函數的 contrast 參數。
這里用的是 glmQLFTest 而不是 glmLRT 是因為前面用了glmQLTFit進行擬合，所以需要用QL F-test進行檢驗。如果前面用的是 glmFit ，那麼對應的就是 glmLRT . 作者稱QL F-test更加嚴格。多重試驗矯正用的也是BH方法。
後續就是提取顯著性差異的基因用作下游分析，做一些圖看看
第六步：差異表達檢驗（2）。上面找到的顯著性差異的基因，沒有考慮效應值，也就是具體變化了多少倍。我們也可用找表達量變化比較大的基因，對應的函數是 glmTreat 。
經過上面兩個方法的洗禮，基本上套路你也就知道了，我先簡單小結一下，然後繼續介紹limma包的 voom 。
Limma原先用於處理基因表達晶元數據，可是說是這個領域的老大 :sunglasses: 。如果你仔細看edgeR導入界面，你就會發現，edgeR有一部分功能依賴於limma包。Limma採用經驗貝葉斯模型（ Empirical Bayesian model）讓結果更穩健。
在處理RNA-Seq數據時，raw read count先被轉成log2-counts-per-million (logCPM)，然後對mean-variance關系建模。建模有兩種方法：
數據預處理 ： Limma使用edgeR的DGEList對象，並且過濾方法都是一致的，對應edgeR的第一步,第二步， 第三步
差異表達分析 : 使用」limma-trend「
差異表達分析 : 使用」limma-voom「
如果分析基因晶元數據，必須好好讀懂LIMMA包。
基本上每一個包，我都提取了各種的顯著性基因，比較就需要用韋恩圖了，但是我偏不 :stuck_out_tongue: 我要用UpSetR.
感覺limma的結果有點奇怪，有生之年在折騰吧。
好吧，這部分我鴿了
[1] Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data
[2] https://www.bioconctor.org/help/workflows/rnaseqGene/
[3] https://www.bioconctor.org/help/workflows/RnaSeqGeneEdgeRQL/
[4] https://www.bioconctor.org/help/workflows/RNAseq123/