Contents

1 摘要

简单且高效地分析RNA测序数据的能力正是Bioconductor的核心优势。 RNA-seq分析通常从基因水平的序列计数开始,涉及到数据预处理,探索性数据分析,差异表达检验以及通路分析,得到的结果可用于指导进一步实验和验证研究。 在这篇工作流程文章中,我们通过分析来自小鼠乳腺的RNA测序数据,示范了如何使用流行的edgeR包载入、整理、过滤和归一化数据,然后用limma包的voom方法、线性模型和经验贝叶斯调节(empirical Bayes moderation)来评估差异表达并进行基因集检验。通过使用Glimma包,此流程得到了增进,实现了结果的互动探索,使用户得以查看单个样本与基因。 这三个软件包提供的完整分析突出了研究人员可以使用Bioconductor轻松地从RNA测序实验的原始计数揭示生物学意义。

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工作流程,使用基因水平计数作为输入,经过预处理和探索性数据处理,然后得到差异表达(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)进行基因水平的总结,利用其内置的mm10基于RefSeq的注释)。

这些样本的计数数据可以从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

尽管这九个文本文件可以分别读入R然后合并为一个计数矩阵,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 L008

4.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)、染色体名称和位置(chromosome names and locations)、Entrez基因ID(Entrez gene IDs)、Refseq基因ID(Refseq gene IDs)和Ensembl基因ID(Ensembl gene IDs)等。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关联于染色体chr4chr4_JH584294_random,小RNA Mir5098关联于chr2chr5chr8chr11chr17)。 为了处理重复的基因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那样考虑到基因长度区别的影响,但在我们的分析中经常会用到它们。虽然也可以使用RPKM和FPKM,但CPM和log-CPM只使用计数矩阵即可计算,且已足以满足我们所关注的比较的需要。假设不同条件之间剪接异构体(isoform)的使用没有差别,差异表达分析研究同一基因在不同条件下的表达差异,而不是比较多个基因之间的表达或测定绝对表达量。换而言之,基因长度在我们关注的比较中始终不变,且任何观测到的差异是来自于条件的变化而不是基因长度的变化。

在此处,使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。如果可以提供基因长度信息,可以使用edgeR中的rpkm函数计算RPKM值,就像计算CPM值那样简单。

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函数会在进行log2转换前给CPM值加上一个弥补值。默认的弥补值是2/L,其中2是“预先计数”,而L是样本总序列数(以百万计)的平均值,所以log-CPM值是根据CPM值通过log2(CPM + 2/L)计算得到的。这样的计算方式可以确保任意两个具有相同CPM值的序列片段计数的log-CPM值也相同。弥补值的使用可以避免对零取对数,并能使所有样本间的log倍数变化(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.4
summary(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

在下游的线性模型分析中,使用limmavoom函数时也会用到log-CPM值,但voom会默认使用更小的预先计数重新计算自己的log-CPM值。

5.2 删除低表达基因

所有数据集中都混有表达的基因与不表达的基因。尽管我们想要检测在一种条件中表达但再另一种条件中不表达的基因,也有一些基因在所有样品中都不表达。实际上,这个数据集中19%的基因在所有九个样品中的计数都是零。

table(rowSums(x$counts==0)==9)
## 
## FALSE  TRUE 
## 22026  5153

对log-CPM值的分布绘制的图表显示每个样本中很大一部分基因都是不表达或者表达程度相当低的,它们的log-CPM值非常小甚至是负的(下图A部分)。

在任何样本中都没有足够多的序列片段的基因应该从下游分析中过滤掉。这样做的原因有好几个。 从生物学的角度来看,在任何条件下的表达水平都不具有生物学意义的基因都不值得关注,因此最好忽略。 从统计学的角度来看,去除低表达计数基因使数据中的均值 - 方差关系可以得到更精确的估计,并且还减少了在观察差异表达的下游分析中需要进行的统计检验的数量。

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