這個步驟推薦在R里面做膊存，載入表達矩陣，然后設置好分組信息禁灼，統(tǒng)一用DEseq2進行差異分析管挟，當然也可以走走edgeR或者limma的voom流程。
基本任務是得到差異分析結果弄捕，進階任務是比較多個差異分析結果的異同點僻孝。

目錄

• 數(shù)據(jù)填坑
• 理論基礎：線性模型， 設計矩陣和比較矩陣
• 標準化一二事
• 探索性分析一二事
• 使用DESeq2進行差異基因分析
• 使用edgeR進行差異基因分析
• 使用limma進行差異基因分析
  • 不同軟件包分析結果比較
• 使用GFOLD進行無重復樣本的差異基因分析
• 不同差異表達分析的比較

數(shù)據(jù)填坑

原先三個樣本的HTSeq-count計數(shù)的數(shù)據(jù)可以在我的GitHub中找到守谓，但是前面已經說過Jimmy失誤讓我們分析的人類就只有3個樣本穿铆， 另外一個樣本需要從另一批數(shù)據(jù)獲取（請注意batch effect)斋荞，所以不能保證每一組都有兩個重復荞雏。

我一直堅信”你并不孤獨“這幾個字，遇到這種情況的人肯定不止我一個平酿，于是我找到了幾種解決方法

• 使用edgeR凤优，指定dispersion值
• 無重復轉錄組數(shù)據(jù)推薦用同濟大學的GFOLD

以上方法都會在后續(xù)進行介紹，但是我們DESeq2必須得要有重復的問題亟待解決染服，沒辦法我只能自己瞎編了别洪。雖然是編，我們也要有模有樣柳刮，不能直接復制一份挖垛，要考慮到高通量測序的read是默認符合泊松分布的。我是這樣編的秉颗。

• 計算KD重復組的均值差痢毒，作為泊松分布的均值
• 使用概率函數(shù)rpois()隨機產生一個數(shù)值，前一步的均值作為lambda蚕甥，
• 對一些read count 低于均值的直接加上對應KD重復組之間的差值

# import data if sample are small
options(stringsAsFactors = FALSE)
control <- read.table("F:/Data/RNA-Seq/matrix/SRR3589956.count",
                       sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
                    sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.count",
                    sep="\t",col.names = c("gene_id","rep2"))
# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]

ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL
## the sample problem
delta_mean <- abs(mean(raw_count_filt$rep1) - mean(raw_count_filt$rep2))

sampleNum <- length(raw_count_filt$control)
sampleMean <- mean(raw_count_filt$control)
control2 <- integer(sampleNum)

for (i in 1:sampleNum){
  if(raw_count_filt$control[i] < sampleMean){
    control2[i] <- raw_count_filt$control[i] + abs(raw_count_filt$rep1[i] - raw_count_filt$rep2[i])
  }
  else{
    control2[i] <- raw_count_filt$control[i] + rpois(1,delta_mean)
  }
}
# add data to raw_count
raw_count_filt$control2 <- control2

這僅僅是一種填坑的方法而已哪替，更好模擬數(shù)據(jù)的方法需要參閱更加專業(yè)的文獻， 有生之年 我希望能補上這一個部分菇怀。

理論基礎：線性模型凭舶， 設計矩陣和比較矩陣

這部分內容最先在 RNA-Seq Data Analysis 的8.5.3節(jié)看到，剛開始一點都不理解爱沟，但是學完生物統(tǒng)計之后帅霜，我認為這是理解所有差異基因表達分析R包的關鍵。

基本上呼伸，統(tǒng)計課都會介紹如何使用t檢驗用來比較兩個樣本之間的差異身冀，然后在樣本比較多的時候使用方差分析確定樣本間是否有差異。當然前是樣本來自于正態(tài)分布的群體，或者隨機獨立大量抽樣搂根。

對于基因芯片的差異表達分析而言珍促，由于普遍認為其數(shù)據(jù)是服從正態(tài)分布，因此差異表達分析無非就是用t檢驗和或者方差分析應用到每一個基因上剩愧。高通量一次性找的基因多猪叙，于是就需要對多重試驗進行矯正，控制假陽性隙咸。目前在基因芯片的分析用的最多的就是limma沐悦。

但是，高通量測序(HTS)的read count普遍認為是服從泊松分布（當然有其他不同意見）五督，不可能直接用正態(tài)分布的t檢驗和方差分析。 當然我們可以簡單粗暴的使用對于的非參數(shù)檢驗的方法瓶殃，但是統(tǒng)計力不夠充包，結果的p值矯正之估計一個差異基因都找不到。老板花了一大筆錢遥椿，結果卻說沒有差異基因基矮，是個負結果，于是好幾千經費打了水漂冠场，他肯定是不樂意的家浇。因此，還是得要用參數(shù)檢驗的方法碴裙，于是就要說到方差分析和線性模型之間的關系了钢悲。

線性回歸和方差分析是同一時期發(fā)展出的兩套方法。在我本科階段的田間統(tǒng)計學課程中就介紹用方差分析（ANOVA）分析不同肥料處理后的產量差異舔株，實驗設計如下

這是最簡單的單因素方差分析莺琳，每一個結果都可以看成 yij = ai + u + eij， 其中u是總體均值载慈，ai是每一個處理的差異惭等，eij是隨機誤差。

image

注：方差分析(Analysis of Variance, ANAOVA)名字聽起來好像是檢驗方差办铡，但其實是為了判斷樣本之間的差異是否真實存在辞做，為此需要證明不同處理內的方差顯著性大于不同處理間的方差。

線性回歸 一般是用于量化的預測變量來預測量化的響應變量寡具。比如說體重與身高的關系建模：

image

當然線性回歸也可用處理名義型或有序型因子（也就是離散變量）作為預測變量秤茅，如果要畫圖的話，就是下面這個情況晒杈。

image

如果我們需要通過一個實驗找到不同處理后對照組和控制組的基因變化嫂伞，那么基因表達可以簡單寫成， y = a + b · treament + e。 和之前的 yij = ai + u + eij 相比帖努，你會發(fā)現(xiàn)公式是如此的一致撰豺。 這是因為線性模型和方差分析都是廣義線性模型(generalizing linear models, GLM)在正態(tài)分布的預測變量的特殊形式。而GLM本身只要采用合適的連接函數(shù)是可以處理對任意類型的變量進行建模的拼余。

目前認為read count之間的差異是符合負二項分布污桦，也叫gamma-Possion分布。那么問題來了匙监，如何用GLM或者LM分析兩個處理件的差異呢凡橱？其實可以簡單的用上圖的擬合直線的斜率來解釋，如果不同處理之間存在差異亭姥，那么這個擬合線的斜率必定不為零稼钩，也就是與X軸平行。但是這是一種便于理解的方式（雖然你也未必能理解）达罗，實際更加復雜坝撑，考慮因素更多。

注1 負二向分布有兩個參數(shù)粮揉，均值(mean)和離散值（dispersion). 離散值描述方差偏離均值的程度巡李。泊松分布可以認為是負二向分布的離散值為1，也就是均值等于方差（mean=variance）的情況扶认。
注2 這部分涉及大量的統(tǒng)計學知識侨拦，不懂就用維基百科一個個查清楚。

聊完了線性模型和方差分析辐宾，下面的設計矩陣（design matrix）就很好理解了狱从， 其實就是用來告訴不同的差異分析函數(shù)應該如何對待變量。比如說我們要研究的KD和control之間變化螃概，設計矩陣就是

那么比較矩陣（contrast matrix)就是告訴差異分析函數(shù)應該如何對哪個因素進行比較矫夯， 這里就是比較不同處理下表達量的變化。

標準化一二事

其實read count如何標準化的方法有很多吊洼，最常用的是FPKM和RPKM训貌，雖然它們其實是錯的--FPKM/RPKM是錯的。

我推薦閱讀 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同標準化方法之間的差異冒窍。

有一些方法是要求原始數(shù)據(jù)递沪，有一些則要求經過某類標準化后的數(shù)據(jù)，記得區(qū)分综液。

探索性分析一二事

使用DESeq2進行差異基因分析

關于DESeq2分析差異表達基因款慨，其實在https://www.bioconductor.org/help/workflows/rnaseqGene/ 里面介紹的非常清楚了。

我們已經準備好了count matrix谬莹，接下來就是把數(shù)據(jù)導入DESeq2檩奠。DESeq2導入數(shù)據(jù)的方式有如下4種桩了，基本覆蓋了主流read count軟件的結果。
注 DESeq2要求的數(shù)據(jù)是raw count埠戳， 沒必要進行FPKM/TPM/RPFKM/TMM標準化井誉。

本來我們是可以用DESeq2為htseq-count專門提供的 DESeqDataSetFromHTSeq ，然而很尷尬數(shù)據(jù)不夠要自己湊數(shù)整胃，所以只能改用 DESeqDataSetFromMatrix了 :cold_sweat:

導入數(shù)據(jù)颗圣，構建 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數(shù)據(jù)，節(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)
# 出現(xiàn)如下提示信息，說明運行成功
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

用results獲取結果： results的參數(shù)非常的多除嘹，這里不好具體展開 :pensive: 但是你們會自己看的吧

res <- results()

我們可用mcols查看每一項結果的具體含義，比如說log2FoldChange 表示倍數(shù)變化取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的均值。根據(jù)summary可知覆履，low count的比率很高蹋盆，所以大部分基因表達量不高费薄，也就是集中在0的附近（log2(1)=0，也就是變化1倍）.提供了模型預測系數(shù)的分布總覽栖雾。

下圖是沒有經過 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")
})

image

如果經過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")
})

image

當然還有火山圖，不過留給其他方法作圖岩灭，我們先把差異表達的基因找出來拌倍。

res.deseq2 <- subset(res, padj < 0.05)

一般p value 小于0.05就是顯著了, 顯著性不代表結果正確，只用于給后續(xù)的富集分析和GSEA提供排序標準和篩選而已噪径。關于P值的吐槽簡直無數(shù)柱恤， 請多注意。

使用edgeR進行差異基因分析

edgeR在函數(shù)說明中稱其不但可以分析SAGE找爱， CAGE的RNA-Seq梗顺，Tag-RNA，或RNA-seq车摄， 也能分析ChIP-Seq和CRISPR得到的read counts數(shù)據(jù)寺谤。嗯，我信了:confused:吮播！

edgeR使用DGEList函數(shù)讀取count matrix數(shù)據(jù)变屁，也就說你需要提供一個現(xiàn)成的matrix數(shù)據(jù)，而不是指望它能讀取單獨的文件意狠，然后進行合并(當然機智的我發(fā)現(xiàn)粟关，其實可以用 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數(shù)據(jù)环戈。與DESeq2的預過濾不同闷板，DESeq2的預過濾只是為了改善后續(xù)運算性能，在運行過程中依舊會自動處理low count數(shù)據(jù)院塞，edgeR需要在分析前就要排除那些low count數(shù)據(jù)遮晚，而且非常嚴格。從生物學角度拦止，有生物學意義的基因的表達量必須高于某一個閾值县遣。從統(tǒng)計學角度上， low count的數(shù)據(jù)不太可能有顯著性差異创泄，而且在多重試驗矯正階段還會拖后腿艺玲。 綜上所訴，放心大膽的過濾吧鞠抑。

根據(jù)經驗（又是經驗 :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數(shù) /1000000)，keep.lib.size=FALSE表示重新計算文庫大小兴垦。

第三步： 根據(jù)組成偏好(composition bias)標準化徙赢。edgeR的calcNormFactors函數(shù)使用TMM算法對DGEList標準化

genelist.norm <- calcNormFactors(genelist.filted)

注 大部分的mRNA-Seq數(shù)據(jù)分析用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參數(shù)鲤氢。 edgeR的線性模型和差異表達分析需要定義一個實驗設計矩陣搀擂。很直白的就能發(fā)現(xiàn)是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)需要均值和離散值兩個參數(shù)哨颂。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函數(shù)。當不存在實驗設計矩陣(design matrix)的時候蚂子，estimateDisp 等價于 estimateCommonDisp 和 estimateTagwiseDisp 。而當給定實驗設計矩陣(design matrix)時缭黔， estimateDisp 等價于 estimateGLMCommonDisp, estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp食茎。 其中tag與gene同義。

注2 其實這里的第三馏谨， 四别渔， 五步對應的就是DESeq2的DESeq包含的2步，標準化和離散值估測惧互。

第六步: 差異表達檢驗（1）哎媚。這一步主要構建比較矩陣，類似于DESeq2中的results函數(shù)的 contrast 參數(shù)喊儡。

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")

image

第六步：差異表達檢驗（2）淤毛。上面找到的顯著性差異的基因今缚，沒有考慮效應值，也就是具體變化了多少倍低淡。我們也可用找表達量變化比較大的基因姓言，對應的函數(shù)是 glmTreat。

tr <- glmTreat(fit, contrast=B.LvsP, lfc=log2(1.5))
s

image

使用limma進行差異分析

經過上面兩個方法的洗禮蔗蹋，基本上套路你也就知道了何荚，我先簡單小結一下，然后繼續(xù)介紹limma包的 voom 纸颜。

• 導入read count兽泣， 保存為專門的對象用于后續(xù)分析
• 原始數(shù)據(jù)過濾，根據(jù)標準化read count 或者 raw count 作為篩選標準
• raw read count 標準化
• 通過各種算法（如經驗貝葉斯胁孙，EM）預測dispersion離散值
• 廣義線性模型擬合數(shù)據(jù)
• 差異分析唠倦，也就是統(tǒng)計檢驗部分

Limma原先用于處理基因表達芯片數(shù)據(jù)，可是說是這個領域的老大 :sunglasses: 涮较。如果你仔細看edgeR導入界面稠鼻，你就會發(fā)現(xiàn)，edgeR有一部分功能依賴于limma包狂票。Limma采用經驗貝葉斯模型（ Empirical Bayesian model）讓結果更穩(wěn)健候齿。

在處理RNA-Seq數(shù)據(jù)時，raw read count先被轉成log2-counts-per-million (logCPM)闺属，然后對mean-variance關系建模慌盯。建模有兩種方法：

• 精確權重法（precision weights）也就是“voom"
• 經驗貝葉斯先驗趨勢（empirical Bayes prior trend），也就是”limma-trend“

數(shù)據(jù)預處理： 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))

如果分析基因芯片數(shù)據(jù)国瓮，必須好好讀懂LIMMA包灭必。

不同軟件包分析結果比較

基本上每一個包，我都提取了各種的顯著性基因乃摹，比較就需要用韋恩圖了禁漓，但是我偏不 :stuck_out_tongue: 我要用UpSetR.

library(UpSetR)
input <- fromList(list(edgeR=rownames(sig.edger), DESeq2=rownames(sig.deseq2), limma=rownames(sig.limma)))

image

感覺limma的結果有點奇怪，有生之年在折騰吧孵睬。

使用GFOLD進行無重復樣本的差異基因分析

好吧播歼，這部分我鴿了

參考文件

[1] Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data

[2] https://www.bioconductor.org/help/workflows/rnaseqGene/

[3] https://www.bioconductor.org/help/workflows/RnaSeqGeneEdgeRQL/

[4] https://www.bioconductor.org/help/workflows/RNAseq123/