使用limma、Glimma和edgeR，RNA-seq数据分析易如反掌
2018年12月17日
Contents
1 摘要
简单且高效地分析RNA测序数据的能力正是Bioconductor的核心优势之一。在获得RNA-seq基因表达矩阵后，通常需要对数据进行预处理、探索性数据分析、差异表达检验以及通路分析，以得到可以帮助进一步实验和验证研究的结果。在本工作流程中，我们将通过分析来自小鼠乳腺的RNA测序数据，演示如何使用流行的edgeR包载入、整理、过滤和归一化数据，然后用limma包的voom方法、线性模型和经验贝叶斯调节来评估差异表达并进行基因集检验。通过Glimma包，本流程进一步实现了结果的互动探索，便于用户查看特定样本与基因的分析结果。通过使用这三个Bioconductor包，研究者可以轻松地运行完整的RNA-seq数据分析流程，从原始计数（raw counts）中挖掘出其中蕴含的生物学意义。
2 背景介绍
RNA测序（RNA-seq）是用于研究基因表达的重要技术。其中，在基因组规模下检测多条件之间基因的差异表达是研究者最常探究的问题之一。对于RNA-seq数据，来自Bioconductor项目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用于处理此问题的完善的统计学方法。
在这篇文章中，我们描述了一个用于分析RNA-seq数据的edgeR - limma工作流程，使用基因水平的计数（gene-level counts）作为输入，经过预处理和探索性数据分析，然后得到差异表达（DE）基因和基因表达特征（gene signatures）的列表。Glimma包(Su et al. 2017)提供的交互式图表可以同时呈现整体样本层面与单个基因层面的数据，相对静态的R图表而言，更便于我们探索更多的细节。
此工作流程中我们分析的数据来自Sheridan等人的实验（2015）(Sheridan et al. 2015)，它包含三个细胞群，即基底（basal）、管腔祖细胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。细胞群皆分选自雌性处女小鼠的乳腺，每种都设三个生物学重复。RNA样品分三个批次使用Illumina HiSeq 2000进行测序，得到长为100碱基对的单端序列片段。
本文所述的分析流程假设从RNA-seq实验获得的序列片段已经与适当的参考基因组比对，并已经在基因水平上对序列进行了统计计数。在本文条件下，使用Rsubread包提供的基于R的流程将序列片段与小鼠参考基因组（mm10）比对（具体而言，先使用align函数(Liao, Smyth, and Shi 2013)进行比对，然后使用featureCounts (Liao, Smyth, and Shi 2014)函数，利用其内置的基于RefSeq的mm10注释进行基因水平的总结）。
这些样本的计数数据可以从Gene Expression Omnibus (GEO)数据库 http://www.ncbi.nlm.nih.gov/geo/使用GEO序列登记号GSE63310下载。更多关于实验设计和样品制备的信息也可以在GEO使用该登记号查看。
3 初始配置
library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)4 数据整合
4.1 读入计数数据
首先，从https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file下载文件GSE63310_RAW.tar，并从压缩包中解压出相关的文件。下方的代码将完成此步骤，或者您也可以手动进行这一步并继续后续分析。
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
"GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
R.utils::gunzip(i, overwrite=TRUE)每一个文本文件均为对应样品的原始基因水平计数矩阵。需要注意我们的这次分析仅包含了此实验中的basal、LP和ML样品（可见下方所示文件名）。
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
"GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
"GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
"GSM1545545_JMS9-P8c.txt")
read.delim(files[1], nrow=5)## EntrezID GeneLength Count
## 1 497097 3634 1
## 2 100503874 3259 0
## 3 100038431 1634 0
## 4 19888 9747 0
## 5 20671 3130 1相比于分别读入这九个文本文件然后合并为一个计数矩阵，edgeR提供了更方便的途径，使用readDGE函数即可一步完成。得到的DGEList对象中包含一个计数矩阵，它的27179行分别对应每个基因不重复的Entrez基因ID，九列分别对应此实验中的每个样品。
x <- readDGE(files, columns=c(1,3))
class(x)## [1] "DGEList"
## attr(,"package")
## [1] "edgeR"dim(x)## [1] 27179 9如果数据不是每个样品一个文件的形式，而是一个包含所有样品的计数的文件，则可以先将文件读入R，再使用DGEList函数转换为一个DGEList对象。
4.2 组织样品信息
为进行下游分析，需要将有关实验设计的样品信息与计数矩阵的列关联起来。这里需要包括各种对表达水平有影响的实验变量，无论是生物变量还是技术变量。例如，细胞类型（在这个实验中是basal、LP和ML）、基因型（野生型、敲除）、表型（疾病状态、性别、年龄）、样品处理（用药、对照）和批次信息（如果样品是在不同时间点进行收集和分析的，需要记录进行实验的时间）等。
我们的DGEList对象中包含的samples数据框同时存储了细胞类型（group）和批次（测序泳道lane）信息，每种信息都包含三个不同的水平。在x$samples中，程序会自动计算每个样品的文库大小（即样品的总序列计数），归一化系数会被预先设置为1。
为了方便阅读，我们从DGEList对象x的列名中删去了GEO样品ID（GSM*）。
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
samplenames## [1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4" "JMS8-5"
## [8] "JMS9-P7c" "JMS9-P8c"colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples## files group lib.size norm.factors lane
## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004
## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004
## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004
## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006
## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006
## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006
## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006
## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L0084.3 组织基因注释
我们的DGEList对象中的第二个数据框名为genes，用于存储与计数矩阵的行相关联的基因信息。
为检索这些信息，我们可以使用特定物种的注释包，比如小鼠的Mus.musculus (Bioconductor Core Team 2016b)（或人类的Homo.sapiens (Bioconductor Core Team 2016a)）；或者也可以使用biomaRt 包 (Durinck et al. 2005, 2009)，它通过接入Ensembl genome数据库来进行基因注释。
可以检索的信息类型包括基因符号（gene symbols）、基因名称（gene names）、染色体名称和位置、Entrez基因ID、Refseq基因ID和Ensembl基因ID等。biomaRt主要通过Ensembl基因ID进行检索，而Mus.musculus包含来自不同来源的信息，允许用户从不同基因ID中选择某一种作为检索键。
我们使用Mus.musculus包，利用我们数据集中的Entrez基因ID来检索相关的基因符号和染色体信息。
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
keytype="ENTREZID")
head(genes)## ENTREZID SYMBOL TXCHROM
## 1 497097 Xkr4 chr1
## 2 100503874 Gm19938 <NA>
## 3 100038431 Gm10568 <NA>
## 4 19888 Rp1 chr1
## 5 20671 Sox17 chr1
## 6 27395 Mrpl15 chr1与任何基因ID一样，Entrez基因ID可能不能一对一地匹配我们想获得的基因信息。在处理之前，检查重复的基因ID和弄清楚重复的来源非常重要。我们的基因注释中包含28个能匹配到多个不同染色体的基因（比如基因Gm1987关联于染色体chr4和chr4_JH584294_random，小RNA Mir5098关联于chr2，chr5，chr8，chr11和chr17）。 为了处理重复的基因ID，我们可以合并来自多重匹配基因的所有染色体信息，比如将基因Gm1987分配到chr4 and chr4_JH584294_random，或选取其中一条染色体来代表具有重复注释的基因。为了简单起见，我们选择后者，保留每个基因ID第一次出现的信息。
genes <- genes[!duplicated(genes$ENTREZID),]在此例子中，注释与数据对象中的基因顺序是相同的。如果由于缺失和／或重新排列基因ID导致其顺序不一致，我们可以用match函数来正确排序基因。然后，我们将基因注释的数据框添加到DGEList对象，数据的整合就完成了，此时的数据对象中含有原始计数数据以及相关的样品信息和基因注释。
x$genes <- genes
x## An object of class "DGEList"
## $samples
## files group lib.size norm.factors lane
## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004
## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004
## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004
## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006
## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006
## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006
## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006
## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
##
## $counts
## Samples
## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c
## 497097 1 2 342 526 3 3 535 2 0
## 100503874 0 0 5 6 0 0 5 0 0
## 100038431 0 0 0 0 0 0 1 0 0
## 19888 0 1 0 0 17 2 0 1 0
## 20671 1 1 76 40 33 14 98 18 8
## 27174 more rows ...
##
## $genes
## ENTREZID SYMBOL TXCHROM
## 1 497097 Xkr4 chr1
## 2 100503874 Gm19938 <NA>
## 3 100038431 Gm10568 <NA>
## 4 19888 Rp1 chr1
## 5 20671 Sox17 chr1
## 27174 more rows ...5 数据预处理
5.1 原始数据尺度转换
由于更深的测序总会产生更多的序列片段，在差异表达及相关的分析中，我们很少直接使用序列数。在实际操作时，我们通常将原始的序列数进行归一化，来消除测序深度所导致的差异。通常被使用的方法有基于序列的CPM（counts per million）、log-CPM、FPKM（fragments per kilobase of transcript per million），和基于转录本数目的RPKM（reads per kilobase of transcript per million）。
我们在分析中通常使用CPM和log-CPM转换。虽然RPKM和FPKM可以校正基因长度区别的影响，但CPM和log-CPM只使用计数矩阵即可计算，且已足以满足我们所关注的比较的需要。假设不同条件之间剪接异构体（isoform）的表达比例没有变化，差异表达分析关注的是同一基因在不同条件之间表达水平的相对差异，而不是比较多个基因之间的差异或测定绝对表达量。换而言之，基因长度在我们进行比较的不同组之间是始终不变的，且任何观测到的差异都来自于不同组的条件的变化而不是基因长度的变化。
我们使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。如果可以提供基因长度信息，RPKM值的计算也和CPM值的计算一样简单，只需使用edgeR中的rpkm函数。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)对于一个基因，CPM值为1相当于在本实验测序深度最低的样品中（JMS9-P8c, 文库大小约2千万）有20个计数，或者在测序深度最高的样品中（JMS8-3，文库大小约7.6千万）有76个计数。
log-CPM值将被用于探索性图表中。当设置log=TRUE时，cpm函数会给CPM值加上一个弥补值并进行log2转换。默认的弥补值是2/L，其中2是“预先计数”，而L是样本文库大小（以百万计）的平均值，所以log-CPM值是根据CPM值通过log2(CPM + 2/L)计算得到的。这样的计算方式可以确保任意两个具有相同CPM值的序列片段计数的log-CPM值也相同。弥补值的使用可以避免对零取对数，并能使所有样本间的对数倍数变化（log-fold-change）向0推移而减小低表达基因间微小计数变化带来的巨大的伪差异性，这对于绘制探索性图表很有帮助。在这个数据集中，平均的样本文库大小是4.55千万，所以L约等于45.5，且每个样本中的最小log-CPM值为log2(2/45.5) = -4.51。换而言之，在加上了预先计数弥补值后，此数据集中的零表达计数对应的log-CPM值为-4.51：
L <- mean(x$samples$lib.size) * 1e-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)## [1] 45.5 51.4summary(lcpm)## 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3
## Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51
## 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51
## Median :-0.68 Median :-0.36 Median :-0.10 Median :-0.09 Median :-0.43
## Mean : 0.17 Mean : 0.33 Mean : 0.44 Mean : 0.41 Mean : 0.32
## 3rd Qu.: 4.29 3rd Qu.: 4.56 3rd Qu.: 4.60 3rd Qu.: 4.55 3rd Qu.: 4.58
## Max. :14.76 Max. :13.50 Max. :12.96 Max. :12.85 Max. :12.96
## JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c
## Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51
## 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51
## Median :-0.41 Median :-0.07 Median :-0.17 Median :-0.33
## Mean : 0.25 Mean : 0.40 Mean : 0.37 Mean : 0.27
## 3rd Qu.: 4.32 3rd Qu.: 4.43 3rd Qu.: 4.60 3rd Qu.: 4.44
## Max. :14.85 Max. :13.19 Max. :12.94 Max. :14.01在接下来的线性模型分析中，使用limma的voom函数时也会用到log-CPM值，但voom会默认使用更小的预先计数重新计算自己的log-CPM值。
5.2 删除低表达基因
所有数据集中都混有表达的基因与不表达的基因。我们想要检测在一种条件中表达但在另一种条件中不表达的基因，但也有一些基因在所有样品中都不表达。实际上，这个数据集中19%的基因在所有九个样品中的计数都是零。
table(rowSums(x$counts==0)==9)##
## FALSE TRUE
## 22026 5153log-CPM值的分布图表显示每个样本中很大一部分基因都是不表达或者表达程度相当低的，它们的log-CPM值非常小甚至是负的（图1A）。
在任何样本中都没有足够多的序列片段的基因应该从下游分析中过滤掉。这样做的原因有好几个。 从生物学的角度来看，在任何条件下的表达水平都不具有生物学意义的基因都不值得关注，因此最好忽略。 从统计学的角度来看，去除低表达计数基因使数据中的均值 - 方差关系可以得到更精确的估计，并且还减少了下游的差异表达分析中需要进行的统计检验的数量。
edgeR包中的filterByExpr函数提供了自动过滤基因的方法，可保留尽可能多的有足够表达计数的基因。
keep.exprs <- filterByExpr(x, group=group)
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)## [1] 16624 9此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多计数的基因。实际进行过滤时，使用的是CPM值而不是表达计数，以避免对总序列数大的样本的偏向性。在这个数据集中，总序列数的中位数是5.1千万，且10/51约等于0.2，所以filterByExpr函数保留在至少三个样本中CPM值大于等于0.2的基因。在我们的此次实验中，一个具有生物学意义的基因需要在至少三个样本中表达，因为三种细胞类型组内各有三个重复。过滤的阈值取决于测序深度和实验设计。如果样本总表达计数量增大，那么可以选择更低的CPM阈值，因为更大的总表达计数量提供了更好的分辨率来探究更多表达水平更低的基因。
使用这个标准，基因的数量减少到了16624个，约为开始时数量的60%。过滤后的log-CPM值显示出每个样本的分布基本相同（下图B部分）。需要注意的是，从整个DGEList对象中取子集时同时删除了被过滤的基因的计数和其相关的基因信息。留下的基因相对应的基因信息和计数在过滤后的DGEList对象中被保留。
下方给出的是绘图所用代码。
lcpm.cutoff <- log2(10/M + 2/L)
library(RColorBrewer)
nsamples <- ncol(x)
col <- brewer.pal(nsamples, "Paired")
par(mfrow=c(1,2))
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")
lcpm <- cpm(x, log=TRUE)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")