這個(gè)步驟推薦在R里面做，載入表達(dá)矩陣，然后設(shè)置好分組信息，統(tǒng)一用DEseq2進(jìn)行差異分析，當(dāng)然也可以走走edgeR或者limma的voom流程。

基本任務(wù)是得到差異分析結(jié)果，進(jìn)階任務(wù)是比較多個(gè)差異分析結(jié)果的異同點(diǎn)。

目錄

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

理論基礎(chǔ)：線性模型， 設(shè)計(jì)矩陣和比較矩陣

標(biāo)準(zhǔn)化一二事

探索性分析一二事（沒(méi)寫(xiě)）

使用DESeq2進(jìn)行差異基因分析

使用edgeR進(jìn)行差異基因分析

使用limma進(jìn)行差異基因分析

不同軟件包分析結(jié)果比較

使用GFOLD進(jìn)行無(wú)重復(fù)樣本的差異基因分析（沒(méi)寫(xiě)）

不同差異表達(dá)分析的比較數(shù)據(jù)填坑

原先三個(gè)樣本的HTSeq-count計(jì)數(shù)的數(shù)據(jù)可以在我的GitHub中找到，但是前面已經(jīng)說(shuō)過(guò)Jimmy失誤讓我們分析的人類就只有3個(gè)樣本， 另外一個(gè)樣本需要從另一批數(shù)據(jù)獲取（請(qǐng)注意batch effect)，所以不能保證每一組都有兩個(gè)重復(fù)。

我一直堅(jiān)信”你并不孤獨(dú)“這幾個(gè)字，遇到這種情況的人肯定不止我一個(gè)，于是我找到了幾種解決方法

使用edgeR，指定dispersion值

無(wú)重復(fù)轉(zhuǎn)錄組數(shù)據(jù)推薦用同濟(jì)大學(xué)的GFOLD

以上方法都會(huì)在后續(xù)進(jìn)行介紹，但是我們DESeq2必須得要有重復(fù)的問(wèn)題亟待解決，沒(méi)辦法我只能自己瞎編了。雖然是編，我們也要有模有樣，不能直接復(fù)制一份，要考慮到高通量測(cè)序的read是默認(rèn)符合泊松分布的。我是這樣編的。

計(jì)算KD重復(fù)組的均值差，作為泊松分布的均值

使用概率函數(shù)rpois()隨機(jī)產(chǎn)生一個(gè)數(shù)值，前一步的均值作為lambda，

對(duì)一些read count 低于均值的直接加上對(duì)應(yīng)KD重復(fù)組之間的差值

# import data if sample are small

options(stringsAsFactors = FALSE)control <->'F:/Data/RNA-Seq/matrix/SRR3589956.count',                       sep='\t', col.names = c('gene_id','control'))rep1 <->'F:/Data/RNA-Seq/matrix/SRR3589957.count',                    sep='\t', col.names = c('gene_id','rep1'))rep2 <->'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 <->1:-5,]ENSEMBL <->'(.*?)\\.\\d*?_\\d', '\\1', raw_count_filt$gene_id)row.names(raw_count_filt) <->

## 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=""><->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]="" +="">1,delta_mean)  }}

# add data to raw_count

raw_count_filt$control2 <->

這僅僅是一種填坑的方法而已，更好模擬數(shù)據(jù)的方法需要參閱更加專業(yè)的文獻(xiàn)， 有生之年 我希望能補(bǔ)上這一個(gè)部分。理論基礎(chǔ)：線性模型， 設(shè)計(jì)矩陣和比較矩陣

這部分內(nèi)容最先在 RNA-Seq Data Analysis 的8.5.3節(jié)看到，剛開(kāi)始一點(diǎn)都不理解，但是學(xué)完生物統(tǒng)計(jì)之后，我認(rèn)為這是理解所有差異基因表達(dá)分析R包的關(guān)鍵。

基本上，統(tǒng)計(jì)課都會(huì)介紹如何使用t檢驗(yàn)用來(lái)比較兩個(gè)樣本之間的差異，然后在樣本比較多的時(shí)候使用方差分析確定樣本間是否有差異。當(dāng)然前是樣本來(lái)自于正態(tài)分布的群體，或者隨機(jī)獨(dú)立大量抽樣。

對(duì)于基因芯片的差異表達(dá)分析而言，由于普遍認(rèn)為其數(shù)據(jù)是服從正態(tài)分布，因此差異表達(dá)分析無(wú)非就是用t檢驗(yàn)和或者方差分析應(yīng)用到每一個(gè)基因上。高通量一次性找的基因多，于是就需要對(duì)多重試驗(yàn)進(jìn)行矯正，控制假陽(yáng)性。目前在基因芯片的分析用的最多的就是limma。

但是，高通量測(cè)序(HTS)的read count普遍認(rèn)為是服從泊松分布（當(dāng)然有其他不同意見(jiàn)），不可能直接用正態(tài)分布的t檢驗(yàn)和方差分析。 當(dāng)然我們可以簡(jiǎn)單粗暴的使用對(duì)于的非參數(shù)檢驗(yàn)的方法，但是統(tǒng)計(jì)力不夠，結(jié)果的p值矯正之估計(jì)一個(gè)差異基因都找不到。老板花了一大筆錢(qián)，結(jié)果卻說(shuō)沒(méi)有差異基因，是個(gè)負(fù)結(jié)果，于是好幾千經(jīng)費(fèi)打了水漂，他肯定是不樂(lè)意的。因此，還是得要用參數(shù)檢驗(yàn)的方法，于是就要說(shuō)到方差分析和線性模型之間的關(guān)系了。

線性回歸和方差分析是同一時(shí)期發(fā)展出的兩套方法。在我本科階段的田間統(tǒng)計(jì)學(xué)課程中就介紹用方差分析（ANOVA）分析不同肥料處理后的產(chǎn)量差異，實(shí)驗(yàn)設(shè)計(jì)如下肥料重復(fù)1重復(fù)2重復(fù)3重復(fù)4

A1

…

…

…

…

A2

…

…

…

…

A3

…

…  …

…

…

這是最簡(jiǎn)單的單因素方差分析，每一個(gè)結(jié)果都可以看成 yij = ai + u + eij， 其中u是總體均值，ai是每一個(gè)處理的差異，eij是隨機(jī)誤差。

注：方差分析(Analysis of Variance, ANAOVA)名字聽(tīng)起來(lái)好像是檢驗(yàn)方差，但其實(shí)是為了判斷樣本之間的差異是否真實(shí)存在，為此需要證明不同處理內(nèi)的方差顯著性大于不同處理間的方差。

線性回歸 一般是用于量化的預(yù)測(cè)變量來(lái)預(yù)測(cè)量化的響應(yīng)變量。比如說(shuō)體重與身高的關(guān)系建模：

當(dāng)然線性回歸也可用處理名義型或有序型因子（也就是離散變量）作為預(yù)測(cè)變量，如果要畫(huà)圖的話，就是下面這個(gè)情況。

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

目前認(rèn)為read count之間的差異是符合負(fù)二項(xiàng)分布，也叫g(shù)amma-Possion分布。那么問(wèn)題來(lái)了，如何用GLM或者LM分析兩個(gè)處理件的差異呢？其實(shí)可以簡(jiǎn)單的用上圖的擬合直線的斜率來(lái)解釋，如果不同處理之間存在差異，那么這個(gè)擬合線的斜率必定不為零，也就是與X軸平行。但是這是一種便于理解的方式（雖然你也未必能理解），實(shí)際更加復(fù)雜，考慮因素更多。

注1 負(fù)二向分布有兩個(gè)參數(shù)，均值(mean)和離散值（dispersion). 離散值描述方差偏離均值的程度。泊松分布可以認(rèn)為是負(fù)二向分布的離散值為1，也就是均值等于方差（mean=variance）的情況。

注2 這部分涉及大量的統(tǒng)計(jì)學(xué)知識(shí)，不懂就用維基百科一個(gè)個(gè)查清楚。

聊完了線性模型和方差分析，下面的設(shè)計(jì)矩陣（design matrix）就很好理解了， 其實(shí)就是用來(lái)告訴不同的差異分析函數(shù)應(yīng)該如何對(duì)待變量。比如說(shuō)我們要研究的KD和control之間變化，設(shè)計(jì)矩陣就是樣本處理

sample1

control

sample2

control

sample3

KD

sample4

KD

那么比較矩陣（contrast matrix)就是告訴差異分析函數(shù)應(yīng)該如何對(duì)哪個(gè)因素進(jìn)行比較， 這里就是比較不同處理下表達(dá)量的變化。標(biāo)準(zhǔn)化一二事

其實(shí)read count如何標(biāo)準(zhǔn)化的方法有很多，最常用的是FPKM和RPKM，雖然它們其實(shí)是錯(cuò)的—FPKM/RPKM是錯(cuò)的。

我推薦閱讀 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同標(biāo)準(zhǔn)化方法之間的差異。

有一些方法是要求原始數(shù)據(jù)，有一些則要求經(jīng)過(guò)某類標(biāo)準(zhǔn)化后的數(shù)據(jù)，記得區(qū)分。使用DESeq2進(jìn)行差異基因分析

關(guān)于DESeq2分析差異表達(dá)基因，其實(shí)在https://www./help/workflows/rnaseqGene/ 里面介紹的非常清楚了。

我們已經(jīng)準(zhǔn)備好了count matrix，接下來(lái)就是把數(shù)據(jù)導(dǎo)入DESeq2。DESeq2導(dǎo)入數(shù)據(jù)的方式有如下4種，基本覆蓋了主流read count軟件的結(jié)果。

注 DESeq2要求的數(shù)據(jù)是raw count， 沒(méi)必要進(jìn)行FPKM/TPM/RPFKM/TMM標(biāo)準(zhǔn)化。functionpackageframeworkoutputDESeq2 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

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

導(dǎo)入數(shù)據(jù)，構(gòu)建 DESeq2 所需的 DESeqDataSet 對(duì)象

library(DESeq2)countData <->2:5]condition <->'control','KD','KD','control'))dds <- deseqdatasetfrommatrix(countdata,="" dataframe(condition),="" design="~" condition="">

注： 這一步到下一步之間可以過(guò)濾掉一些low count數(shù)據(jù)，節(jié)省內(nèi)存，提高運(yùn)行速度

nrow(dds)dds <- dds[="" rowsums(counts(dds))=""> 1, ]nrow(dds)

使用DESeq進(jìn)行差異表達(dá)分析： DESeq包含三步，estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons)， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest），可以分布運(yùn)行，也可用一步到位，最后返回 results可用的DESeqDataSet對(duì)象。

dds <->

# 出現(xiàn)如下提示信息，說(shuō)明運(yùn)行成功

estimating size factorsestimating dispersionsgene-wise dispersion estimatesmean-dispersion relationshipfinal dispersion estimatesfitting model and testing

用results獲取結(jié)果： results的參數(shù)非常的多，這里不好具體展開(kāi)

但是你們會(huì)自己看的吧

res <->

我們可用mcols查看每一項(xiàng)結(jié)果的具體含義，比如說(shuō)log2FoldChange 表示倍數(shù)變化取log2結(jié)果，還能畫(huà)個(gè)火山圖。一般簡(jiǎn)單粗暴的用2到3倍作為閾值，但是對(duì)于低表達(dá)的基因，3倍也是噪音，那些高表達(dá)的基因，1.1倍都是生物學(xué)顯著了。更重要的沒(méi)有考慮到組內(nèi)變異，沒(méi)有統(tǒng)計(jì)學(xué)意義。padj 就是用BH對(duì)多重試驗(yàn)進(jìn)行矯正。

mcols(res, use.names = TRUE)DataFrame with 6 rows and 2 columns                       type                                     description                                                   baseMean       intermediate       mean of normalized counts for all sampleslog2FoldChange      results log2 fold change (MLE): condition KD vs controllfcSE               results         standard error: condition KD vs controlstat                results         Wald statistic: condition KD vs controlpvalue              results      Wald test p-value: condition KD vs controlpadj                results                            BH adjusted p-values

用summary看描述性的結(jié)果，大致是上調(diào)的基因占總體的11%，下調(diào)的是7.1%（KD vs control）

summary(res)out of 29469 with nonzero total read countadjusted p-value < 0.1lfc=""> 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="" esults[2]="" see="" 'independentfiltering'="" argument="" of="">

畫(huà)個(gè)MA圖，還能標(biāo)注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，衡量基因表達(dá)量變化，上調(diào)還是下調(diào)。A表示每個(gè)基因的count的均值。根據(jù)summary可知，low count的比率很高，所以大部分基因表達(dá)量不高，也就是集中在0的附近（log2(1)=0，也就是變化1倍）.提供了模型預(yù)測(cè)系數(shù)的分布總覽。

下圖是沒(méi)有經(jīng)過(guò) 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')})

如果經(jīng)過(guò)lfcShrink 收縮log2 fold change， 結(jié)果會(huì)好看很多

res.shrink <- lfcshrink(dds,="" contrast="">'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')})

當(dāng)然還有火山圖，不過(guò)留給其他方法作圖，我們先把差異表達(dá)的基因找出來(lái)。

res.deseq2 <- subset(res,="" padj=""><>0.05)

一般p value 小于0.05就是顯著了, 顯著性不代表結(jié)果正確，只用于給后續(xù)的富集分析和GSEA提供排序標(biāo)準(zhǔn)和篩選而已。關(guān)于P值的吐槽簡(jiǎn)直無(wú)數(shù)， 請(qǐng)多注意。使用edgeR進(jìn)行差異基因分析

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

edgeR使用DGEList函數(shù)讀取count matrix數(shù)據(jù)，也就說(shuō)你需要提供一個(gè)現(xiàn)成的matrix數(shù)據(jù)，而不是指望它能讀取單獨(dú)的文件，然后進(jìn)行合并(當(dāng)然機(jī)智的我發(fā)現(xiàn)，其實(shí)可以用 tximport 或 DESeqDataSetFromHTSeq 讀取單獨(dú)的文件，然后傳遞給DGEList)

第一步： 構(gòu)建DGEList對(duì)象

library(edgeR)group <->'control','KD','KD','control'))genelist <- dgelist(counts="">2:5], group = group)

第二步： 過(guò)濾 low counts數(shù)據(jù)。與DESeq2的預(yù)過(guò)濾不同，DESeq2的預(yù)過(guò)濾只是為了改善后續(xù)運(yùn)算性能，在運(yùn)行過(guò)程中依舊會(huì)自動(dòng)處理low count數(shù)據(jù)，edgeR需要在分析前就要排除那些low count數(shù)據(jù)，而且非常嚴(yán)格。從生物學(xué)角度，有生物學(xué)意義的基因的表達(dá)量必須高于某一個(gè)閾值。從統(tǒng)計(jì)學(xué)角度上， low count的數(shù)據(jù)不太可能有顯著性差異，而且在多重試驗(yàn)矯正階段還會(huì)拖后腿。 綜上所訴，放心大膽的過(guò)濾吧。

根據(jù)經(jīng)驗(yàn)（又是經(jīng)驗(yàn) :dog: ）， 基因至少在某一些文庫(kù)的count超過(guò)10 ~ 15 才被認(rèn)為是表達(dá)。這一步全靠嘗試， 剔除太多就緩緩，剔除太少就嚴(yán)格點(diǎn)。 我們可以簡(jiǎn)單的對(duì)每個(gè)基因的raw count進(jìn)行比較，但是建議用CPM（count-per-million)標(biāo)準(zhǔn)化 后再比較，避免了文庫(kù)大小的影響。

# 簡(jiǎn)單粗暴的方法

keep <- rowsums(genelist$count)=""> 50

# 利用CPM標(biāo)準(zhǔn)化

keep <- rowsums(cpm(genelist)=""> 0.5 ) >=2table(keep)genelist.filted <- genelist[keep,="" ,keep.lib.sizes="">FALSE]

這里的0.5(即閾值）等于 10/(最小的文庫(kù)的 read count數(shù) /1000000)，keep.lib.size=FALSE表示重新計(jì)算文庫(kù)大小。

第三步： 根據(jù)組成偏好(composition bias)標(biāo)準(zhǔn)化。edgeR的calcNormFactors函數(shù)使用TMM算法對(duì)DGEList標(biāo)準(zhǔn)化

genelist.norm <->

注 大部分的mRNA-Seq數(shù)據(jù)分析用TMM標(biāo)準(zhǔn)化就行了，但是也有例外，比如說(shuō)single-cell RNA-Seq(Lun, Bach, and Marioni 2016), 還有就是global differential expression， 基因組一半以上的基因都是差異表達(dá)的，請(qǐng)盡力避免，(D. Wu et al. 2013)， 不然就需要用到內(nèi)參進(jìn)行標(biāo)準(zhǔn)化了(Risso et al. 2014).

第四步： 實(shí)驗(yàn)設(shè)計(jì)矩陣(Design matrix)， 類似于DESeq2中的design參數(shù)。 edgeR的線性模型和差異表達(dá)分析需要定義一個(gè)實(shí)驗(yàn)設(shè)計(jì)矩陣。很直白的就能發(fā)現(xiàn)是1vs0

design <->0+group)colnames(design) <- levels(group)design=""  control="">

1       1  0

2       0  1

3       0  1

4       1  0

第五步： 估計(jì)離散值（Dispersion）。前面已經(jīng)提到負(fù)二項(xiàng)分布（negative binomial，NB)需要均值和離散值兩個(gè)參數(shù)。edgeR對(duì)每個(gè)基因都估測(cè)一個(gè)經(jīng)驗(yàn)貝葉斯穩(wěn)健離散值（mpirical Bayes moderated dispersion），還有一個(gè)公共離散值（common dispersion，所有基因的經(jīng)驗(yàn)貝葉斯穩(wěn)健離散值的均值）以及一個(gè)趨勢(shì)離散值

genelist.Disp <- estimatedisp(genelist.norm,="" design,="" robust="">TRUE)plotBCV(genelist.Disp)

還可以進(jìn)一步通過(guò)quasi-likelihood (QL)擬合NB模型，用于解釋生物學(xué)和技術(shù)性導(dǎo)致的基因特異性變異 (Lund et al. 2012; Lun, Chen, and Smyth 2016).

fit <- glmqlfit(genelist.disp,="" design,="" robust="">TRUE)head(fit$coefficients)

注1 估計(jì)離散值這個(gè)步驟其實(shí)有許多estimate*Disp函數(shù)。當(dāng)不存在實(shí)驗(yàn)設(shè)計(jì)矩陣(design matrix)的時(shí)候，estimateDisp 等價(jià)于 estimateCommonDisp 和 estimateTagwiseDisp 。而當(dāng)給定實(shí)驗(yàn)設(shè)計(jì)矩陣(design matrix)時(shí)， estimateDisp 等價(jià)于 estimateGLMCommonDisp, estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp。 其中tag與gene同義。

注2 其實(shí)這里的第三， 四， 五步對(duì)應(yīng)的就是DESeq2的DESeq包含的2步，標(biāo)準(zhǔn)化和離散值估測(cè)。

第六步: 差異表達(dá)檢驗(yàn)（1）。這一步主要構(gòu)建比較矩陣，類似于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是因?yàn)榍懊嬗昧薵lmQLTFit進(jìn)行擬合，所以需要用QL F-test進(jìn)行檢驗(yàn)。如果前面用的是glmFit，那么對(duì)應(yīng)的就是glmLRT. 作者稱QL F-test更加嚴(yán)格。多重試驗(yàn)矯正用的也是BH方法。

后續(xù)就是提取顯著性差異的基因用作下游分析，做一些圖看看

topTags(res,n=10)is.de <- decidetestsdge(res)summary(is.de)plotmd(res,="" status="is.de," values="">1,-1), col=c('red','blue'),       legend='topright')

第六步：差異表達(dá)檢驗(yàn)（2）。上面找到的顯著性差異的基因，沒(méi)有考慮效應(yīng)值，也就是具體變化了多少倍。我們也可用找表達(dá)量變化比較大的基因，對(duì)應(yīng)的函數(shù)是 glmTreat。

tr <- glmtreat(fit,="" contrast="B.LvsP," lfc="">1.5))plotMD(tr status=is.de, values=c(1,-1), col=c('red','blue'),       legend='topright')

使用limma進(jìn)行差異分析

經(jīng)過(guò)上面兩個(gè)方法的洗禮，基本上套路你也就知道了，我先簡(jiǎn)單小結(jié)一下，然后繼續(xù)介紹limma包的 voom 。

導(dǎo)入read count， 保存為專門(mén)的對(duì)象用于后續(xù)分析

原始數(shù)據(jù)過(guò)濾，根據(jù)標(biāo)準(zhǔn)化read count 或者 raw count 作為篩選標(biāo)準(zhǔn)

raw read count 標(biāo)準(zhǔn)化

通過(guò)各種算法（如經(jīng)驗(yàn)貝葉斯，EM）預(yù)測(cè)dispersion離散值

廣義線性模型擬合數(shù)據(jù)

差異分析，也就是統(tǒng)計(jì)檢驗(yàn)部分

Limma原先用于處理基因表達(dá)芯片數(shù)據(jù)，可是說(shuō)是這個(gè)領(lǐng)域的老大 :sunglasses: 。如果你仔細(xì)看edgeR導(dǎo)入界面，你就會(huì)發(fā)現(xiàn)，edgeR有一部分功能依賴于limma包。Limma采用經(jīng)驗(yàn)貝葉斯模型（ Empirical Bayesian model）讓結(jié)果更穩(wěn)健。

在處理RNA-Seq數(shù)據(jù)時(shí)，raw read count先被轉(zhuǎn)成log2-counts-per-million (logCPM)，然后對(duì)mean-variance關(guān)系建模。建模有兩種方法：

精確權(quán)重法（precision weights）也就是“voom”

經(jīng)驗(yàn)貝葉斯先驗(yàn)趨勢(shì)（empirical Bayes prior trend），也就是”limma-trend“

數(shù)據(jù)預(yù)處理： Limma使用edgeR的DGEList對(duì)象，并且過(guò)濾方法都是一致的，對(duì)應(yīng)edgeR的第一步,第二步， 第三步

library(edgeR)library(limma)group <->'control','KD','KD','control'))genelist <- dgelist(counts="">2:5], group = group)

### filter base  use CPM

keep <- rowsums(cpm(genelist)=""> 0.5 ) >=2table(keep)genelist.filted <- genelist[keep,="" ,keep.lib.sizes="">FALSE]

### normalizaition

x <- calcnormfactors(x,="" method="">'TMM')

差異表達(dá)分析: 使用”limma-trend“

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

差異表達(dá)分析: 使用”limma-voom“

### DGE with voom

v <- voom(genelist.norm,="" design,="" plot="">TRUE)

#v <- voom(counts,="" design,="" plot="">

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="">1.2))topTreat(fit, coef=ncol(design))

如果分析基因芯片數(shù)據(jù)，必須好好讀懂LIMMA包。不同軟件包分析結(jié)果比較

基本上每一個(gè)包，我都提取了各種的顯著性基因，比較就需要用韋恩圖了，但是我偏不，

我要用UpSetR.

library(UpSetR)input <- fromlist(list(edger="rownames(sig.edger)," deseq2="rownames(sig.deseq2)," limma="">

感覺(jué)limma的結(jié)果有點(diǎn)奇怪，有生之年在折騰吧。