RNA-seq数据基因水平表达差异分析

~ admin

这篇博文来自于Bioc2014的学习笔记。

使用Bioconductor分析RNA-seq数据，获取基因水平表达差异结果主要分三步：记数(counts) -> 质量控制(diagnostics) -> 统计分析(statistics)

在Bioconductor中，用于RNA-seq数据分析的包很多，这里介绍三种edgeR, DESeq2以及voom。

edgeR

为了避免意外，我们开始一个全新的r session。

rm(list=ls())

载入edgeR。

library(edgeR)
## Loading required package: limma

载入数据

load(url("http://qiubio.com/bioconductor/RNA-seq/ds1.Rdata"))
ls()
## [1] "counts"
class(counts)
## [1] "matrix"
dim(counts)
## [1] 30727    27
head(counts[, 1:7], 3)
##                 ES.07985 DE.07981 GT.66339 FG.08004 PE.07980 FE.66350
## ENSG00000254468        0        0        0        1        0        0
## ENSG00000177951       44       50       24       37       38       41
## ENSG00000188076        0        0        0        0        0        0
##                 ES.66342
## ENSG00000254468        1
## ENSG00000177951      162
## ENSG00000188076        0

为载入的数据分组。这里使用列名提供的信息分组。

grp <- as.factor(substr(colnames(counts), 1, 2)) ##substr截取列名的前两个字母
table(grp)
## grp
## DE ES FE FG GT IS PE PH 
##  3  3  4  3  3  2  6  3
o <- order(grp)
pairs(log2(1+counts[,o[1:7]]), pch=".",lower.panel=NULL) ##这个图给我们样品间关系的一揽印象

1

建立一个DGEList对象，并对数据进行标准化(normalize)

d <- DGEList(counts=counts, group=grp)
d <- calcNormFactors(d)
d$samples
##          group lib.size norm.factors
## ES.07985    ES   769440       0.9721
## DE.07981    DE   947798       1.0560
## GT.66339    GT   548227       0.8668
## FG.08004    FG   726088       1.2161
## PE.07980    PE   688903       1.0602
## FE.66350    FE   521573       1.1357
## ES.66342    ES  2281754       0.9036
## DE.66333    DE  2834208       0.9781
## GT.66341    GT  1596550       0.7663
## FG.66346    FG  2070020       1.1650
## PE.66344    PE  1978504       1.0385
## FE.66331    FE  1517276       1.0618
## PH.66332    PH   840380       1.3760
## PH.66336    PH   852548       1.3735
## PE.66345    PE   915042       1.2394
## PE.66330    PE   938173       1.2296
## FE.66348    FE   629687       1.1131
## FE.66334    FE   647313       1.1129
## ES.66335    ES  2048439       0.7027
## DE.66349    DE  1843746       0.7489
## GT.66337    GT  1761089       0.8636
## FG.66351    FG  1417106       0.9454
## PE.66338    PE  2317354       0.8958
## PH.66340    PH  4087566       1.1389
## PE.66329    PE  3700620       1.0793
## IS.66347    IS 12583285       0.5621
## IS.66343    IS 33286768       0.9297

这里的norm.factors计算方法有3种，分别是TMM, RLE和upperquartile。默认使用TMM方法。

接下来要去除低表达值的数据，这一步可以减少很多计算量，不做也不影响最终的结果。

dim(d)
## [1] 30727    27
cps <- cpm(d) ## count per million
k <- rowSums(cps>=1) > 2
d <- d[k,]
dim(d)
## [1] 21707    27

我们看到数据几乎减少了1/3。

然后进行聚类分析。这一步是质量控制的关键一步。如果同一组的biological replicates不 能相聚在一起，必要时需要考虑重新实验。如果条件允许，这一步就要抛弃一些意外的数据。

cols <- as.numeric(d$samples$group)
plotMDS(d,col=cols)

2

绘制multidimensional scaling plot有很多种手段。LFC就是log2 fold change, BCV就是 biological coefficient of variation. 为了快速绘图，默认情况下只是抽取了top 500个值。

par(mfrow=c(2,2))
plotMDS(d, col=cols, main="500 / lLFC")
plotMDS(d, col=cols, method="bcv", main="500 / BCV")
plotMDS(d, col=cols, top=2000, main="2000 / lLFC")
plotMDS(d, col=cols, top=2000, method="bcv", main="2000 / BCV")

3

而后构建模型。

mm <- model.matrix(~-1+grp)
mm
##    grpDE grpES grpFE grpFG grpGT grpIS grpPE grpPH
## 1      0     1     0     0     0     0     0     0
## 2      1     0     0     0     0     0     0     0
## 3      0     0     0     0     1     0     0     0
## 4      0     0     0     1     0     0     0     0
## 5      0     0     0     0     0     0     1     0
## 6      0     0     1     0     0     0     0     0
## 7      0     1     0     0     0     0     0     0
## 8      1     0     0     0     0     0     0     0
## 9      0     0     0     0     1     0     0     0
## 10     0     0     0     1     0     0     0     0
## 11     0     0     0     0     0     0     1     0
## 12     0     0     1     0     0     0     0     0
## 13     0     0     0     0     0     0     0     1
## 14     0     0     0     0     0     0     0     1
## 15     0     0     0     0     0     0     1     0
## 16     0     0     0     0     0     0     1     0
## 17     0     0     1     0     0     0     0     0
## 18     0     0     1     0     0     0     0     0
## 19     0     1     0     0     0     0     0     0
## 20     1     0     0     0     0     0     0     0
## 21     0     0     0     0     1     0     0     0
## 22     0     0     0     1     0     0     0     0
## 23     0     0     0     0     0     0     1     0
## 24     0     0     0     0     0     0     0     1
## 25     0     0     0     0     0     0     1     0
## 26     0     0     0     0     0     1     0     0
## 27     0     0     0     0     0     1     0     0
## attr(,"assign")
## [1] 1 1 1 1 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$grp
## [1] "contr.treatment"

依据模型估算离散度（dispersion）。你一定会问什么是离散度啊？在McCarthy, DJ, Chen, Y, Smyth, GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40, 4288-4297.一文中定义离散度为生物变异系数（biological coefficient of variation）的平方。这里需要解释一下什么是技术重复（technical replicates），什么是生物重复（biological replicates）。技术重复就是同样的一个生物样品，使用相同的操作技术，两次或者两次以上实验。而生物重复是指不同的生物样品，使用相同的操作技术，两次或者两次以上的实验。而 biological variation在上文中专指排除technical variation以外的因为不同生物样品导致的差异。 \(Total CV^2 = Technical CV^2 + Biological CV^2\)

d <- estimateGLMCommonDisp(d,mm)
d <- estimateGLMTrendedDisp(d,mm)
## Loading required package: splines
d <- estimateGLMTagwiseDisp(d,mm)

其实离散度也可以是质量控制的一部分。我们使用plotBCV来察看同一组内不同样品的不同表达水平值的方差分析。

plotBCV(d)

4

从上图中，可以看到红线所在的Common dispersion有0.34左右。但是当我们调用d$common.dispersion时，它却只有0.12。这是怎么回事呢？这是因为求方根之后才是真正的值。

d$common.dispersion
## [1] 0.1182
sqrt(d$common.dispersion)
## [1] 0.3438

＃＃请无视此段。对于BCV的估计，正常情况下，必须要有很多样品才能比较准确，而生物实验是不可能提供成百上千的重复的，所以通常的手段是将表达值相近的基因等同看待。但是这样还是会对BCV过度估计。对于过度估计，在GLM中可以penalized maximum likelihood estimation来进行修正。

对于plotBCV给出的红线，它在0.2~0.4之间都是容易接受的。但如果超过0.5，则需要考虑实验样品之间的差异是否过大了。

plotMeanVar(d, show.raw=TRUE, show.tagwise=TRUE, show.binned=TRUE)

5

在上图中，很明显地显示出基因水平的差异（variance）随着表达值的升高而升高。这一点会让有些人觉得很困惑，因为很多人会认为当表达值升高时，变异系数应该下降才对。但其实随着表达值升高，变异系数下降只符合技术重复带来的差异变化，而生物重复带来的差异变化却并不是这样。从dispersion的图中我们可以看到，它是一个两头翘的分布。因为生物重复差异在整体的变异系数中占主导地位，所以我们会观察到其基因水平的变异系数反而会随着表达值的升高而升高。

对于基因芯片，我们熟悉了MAplot，对于RNAseq同样可以绘制MAplot。MAplot的y轴的中值应该在0线上。

par(mfrow=c(1,2))
plotSmear(d, pair=c("ES","DE"), ylim=c(-5,5))
plotSmear(d, pair=c("DE","GT"), ylim=c(-5,5))

6

接下来的步骤就是统计分析的步骤了。先做广义线性拟合（Generalized linear models），而后构建比较结构（construct contrast），最后依据比较结构进行比较，给出差异表达的基因。

f <- glmFit(d,mm)
con <- makeContrasts("DE-ES"=grpDE-grpES,levels=colnames(mm))
con
##        Contrasts
## Levels  DE-ES
##   grpDE     1
##   grpES    -1
##   grpFE     0
##   grpFG     0
##   grpGT     0
##   grpIS     0
##   grpPE     0
##   grpPH     0
lrt <- glmLRT(f,contrast=con)
topTags(lrt,20)
## Coefficient:  1*grpDE -1*grpES 
##                  logFC logCPM    LR    PValue       FDR
## ENSG00000095596  8.532  6.658 295.5 3.189e-66 6.923e-62
## ENSG00000076356  7.978  7.297 280.2 6.962e-63 7.556e-59
## ENSG00000164292 10.583  6.440 262.0 6.295e-59 4.555e-55
## ENSG00000185823  7.349  5.685 241.4 1.936e-54 1.051e-50
## ENSG00000125798  9.891  8.044 227.5 2.074e-51 9.003e-48
## ENSG00000132130 11.314  6.398 223.1 1.880e-50 6.801e-47
## ENSG00000126005 -4.501  6.464 218.4 2.059e-49 6.384e-46
## ENSG00000104332  4.707  6.784 217.4 3.278e-49 8.895e-46
## ENSG00000158815 14.932  7.295 202.7 5.407e-46 1.304e-42
## ENSG00000119242  6.266  6.655 198.9 3.606e-45 7.828e-42
## ENSG00000182798  5.062  7.548 185.1 3.694e-42 7.289e-39
## ENSG00000266283 11.071  8.362 181.4 2.338e-41 4.230e-38
## ENSG00000120875  7.807  5.529 178.2 1.193e-40 1.991e-37
## ENSG00000110799  7.639  4.880 177.5 1.693e-40 2.624e-37
## ENSG00000121966  7.529  5.850 174.7 6.890e-40 9.971e-37
## ENSG00000185155  6.432  4.827 171.7 3.206e-39 4.349e-36
## ENSG00000164736 12.848  5.288 171.5 3.560e-39 4.546e-36
## ENSG00000124942  6.854  7.487 170.6 5.400e-39 6.512e-36
## ENSG00000141448 13.186  6.579 166.0 5.478e-38 6.259e-35
## ENSG00000110148  6.790  6.114 162.3 3.517e-37 3.817e-34

最后，得到差异表达基因列表后，我们找出几个来做一个检查，看看结果对不对。

cps <- cpm(d)
o <- order(colnames(counts))
barplot( cps["ENSG00000095596",o], col=cols[o], las=2)

7

输出结果至文件。

tt <- topTags(lrt, n=Inf)$table
write.table(tt, file="LRT1.xls", row.names=FALSE, sep="\t", quote=FALSE)

DESeq2 载入DESeq2

library(DESeq2)

我们载入与上一例相同的数据。

load(url("/bioconductor/RNA-seq/ds1.Rdata"))
ls()
##  [1] "cols"   "con"    "counts" "cps"    "d"      "f"      "grp"   
##  [8] "k"      "lrt"    "mm"     "o"      "tt"
head(counts[, 1:7], 3)
##                 ES.07985 DE.07981 GT.66339 FG.08004 PE.07980 FE.66350
## ENSG00000254468        0        0        0        1        0        0
## ENSG00000177951       44       50       24       37       38       41
## ENSG00000188076        0        0        0        0        0        0
##                 ES.66342
## ENSG00000254468        1
## ENSG00000177951      162
## ENSG00000188076        0
grp <- as.factor(substr(colnames(counts), 1, 2)) ##substr截取列名的前两个字母
dds <- DESeqDataSetFromMatrix(counts, colData=data.frame(grp), design=formula(~-1+grp)) ##这里的colData必须是一个DataFrame或者data.frame。每一行都对应着counts中的一列。design中的公式和limma中的方法一致。
design(dds)
## ~-1 + grp
dds <- DESeq(dds) ##一步到位，不需要太多步骤
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res <- results(dds) ##得到结果
res
## log2 fold change (MAP): grp PH vs DE 
## Wald test p-value: grp PH vs DE 
## DataFrame with 30727 rows and 6 columns
##                  baseMean log2FoldChange     lfcSE       stat    pvalue
##                 <numeric>      <numeric> <numeric>  <numeric> <numeric>
## ENSG00000254468    0.4287       -1.05357    1.6490  -0.638900  0.522888
## ENSG00000177951  101.9852        0.30692    0.2534   1.211289  0.225785
## ENSG00000188076    0.1621       -0.01532    1.7535  -0.008735  0.993031
## ENSG00000177947    2.6249       -2.84301    1.2344  -2.303244  0.021265
## ENSG00000254559  496.3435       -0.93363    0.3140  -2.973341  0.002946
## ...                   ...            ...       ...        ...       ...
## ENSG00000248275  19.84377      1.595e+00    0.4163  3.830e+00 0.0001281
## ENSG00000248103   0.01374      2.614e-05    0.8715  3.000e-05 0.9999761
## ENSG00000238035 127.55955     -3.819e-01    0.2806 -1.361e+00 0.1736234
## ENSG00000250765   1.86059     -4.989e-01    1.0464 -4.768e-01 0.6335268
## ENSG00000230178   0.00751      5.954e-07    0.3614  1.648e-06 0.9999987
##                      padj
##                 <numeric>
## ENSG00000254468        NA
## ENSG00000177951   0.37016
## ENSG00000188076        NA
## ENSG00000177947   0.05871
## ENSG00000254559   0.01137
## ...                   ...
## ENSG00000248275 0.0007656
## ENSG00000248103        NA
## ENSG00000238035 0.3038560
## ENSG00000250765 0.7615653
## ENSG00000230178        NA
mcols(res, use.names=TRUE) ##结果中各例的意义
## DataFrame with 6 rows and 2 columns
##                        type                          description
##                 <character>                          <character>
## baseMean       intermediate          the base mean over all rows
## log2FoldChange      results log2 fold change (MAP): grp PH vs DE
## lfcSE               results         standard error: grp PH vs DE
## stat                results         Wald statistic: grp PH vs DE
## pvalue              results      Wald test p-value: grp PH vs DE
## padj                results                 BH adjusted p-values
plotMA(res)

8

plotDispEsts(dds)

9

write.csv(as.data.frame(res), file="results.csv") ##输出结果

voom

library(limma)
v <- voom(d, mm, plot=TRUE) ##这里的d就是edgeR中的d，下面就是标准的limma流程了。

10

vf <- lmFit(v,mm)            # 'mm' 前文定义了
cf <- contrasts.fit(vf,con)  # 'con' 前文定义了
cf <- eBayes(cf)
topTable(cf)                 # 
##                  logFC AveExpr      t   P.Value adj.P.Val     B
## ENSG00000104332  4.682  5.3343  16.82 4.643e-17 1.008e-12 28.60
## ENSG00000182798  5.049  2.3421  16.02 1.802e-16 1.781e-12 27.44
## ENSG00000126005 -4.485  5.7958 -15.47 4.669e-16 1.781e-12 26.29
## ENSG00000147869  7.395  1.6008  15.35 5.733e-16 1.781e-12 26.23
## ENSG00000076356  7.788  5.7036  15.35 5.745e-16 1.781e-12 25.22
## ENSG00000158815 11.806  1.6305  15.55 4.016e-16 1.781e-12 24.92
## ENSG00000164292  9.515  5.0669  15.50 4.418e-16 1.781e-12 24.89
## ENSG00000164736  9.722  0.6449  15.09 9.219e-16 2.501e-12 24.24
## ENSG00000003056 -3.454  8.4616 -13.94 7.795e-15 1.538e-11 23.84
## ENSG00000125798  8.925  7.0947  14.11 5.645e-15 1.361e-11 22.79

我们比较不同方法给出的头几行，差别还真是大啊。 问题来了，我们应该相信哪种方法呢？我只能说，都算一下，然后每种方法给出的结果里都挑几个跑一下qRT-PCR。有人问是不是取多种办法相同的的部分会验证率高些？我的经验，真不一定。如果data好，用哪种方法验证率都高，如果data不行，那只能看看哪种方法更合适了。

发布者：admin

136 thoughts on “RNA-seq数据基因水平表达差异分析”

向慈 小进行回复 取消回复

电子邮件地址不会被公开。 必填项已用*标注