[人工智能]生信入门（二）——使用limma、Glimma和edgeRRNA-seq数据分析

生信入门（二）——使用limma、Glimma和edgeR,RNA-seq数据分析

一、简介

简单且高效地分析RNA 测序数据的能力是生信的核心优势之一。通常在获得RNA-seq基因表达矩阵后，通常需要对数据进行数据预处理、探索性数据分析、差异表达检验以及通路分析，以得到可以帮助进一步试验和验证研究的结果。未来几篇将介绍使用edgeR包，limma包，Glimma包对RNA-seq数据展开分析，从原始计数（counts）中挖掘生物学意义

二、数据背景

在本篇文章中，描述了一个用于分析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碱基对的单端序列片段。

三、初始配置

library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)

四、数据整合

1、数据下载

数据来源：https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file 下载GSE63310_RAW.tar放在一定的文件夹

2、解压文件

setwd("D:\\RData\\prac002")#设置当前路径，下载数据在该文件夹下
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=TURE)

每一个文本文件均为对应样品的原始基因水平计数矩阵。需要注意我们的这次分析仅包含了此实验中的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)#查看导入数据

运行结果

相比于分别读入这九个文本文件然后合并为一个计数矩阵，edgeR提供了更方便的途径，使用readDGE函数即可一步完成。得到的DGEList对象中包含一个计数矩阵，它的27179行分别对应每个基因不重复的Entrez基因ID，九列分别对应此实验中的每个样品。

class(x)
dim(x)

运行结果

注：如果数据不是每个样本一个文件的形式，而是包含所有样本的技术的文件，则可以先将文件读入R，再使用DGEList函数转换为一个DGEList对象

3、组织样本信息

为进行下游分析，需要将有关实验设计的样品信息与计数矩阵的列关联起来。这里需要包括各种对表达水平有影响的实验变量，无论是生物变量还是技术变量。例如，细胞类型（在这个实验中是basal、LP和ML）、基因型（野生型、敲除）、表型（疾病状态、性别、年龄）、样品处理（用药、对照）和批次信息（如果样品是在不同时间点进行收集和分析的，需要记录进行实验的时间）等。
我们的DGEList对象中包含的samples数据框同时存储了细胞类型（group）和批次（测序泳道lane）信息，每种信息都包含三个不同的水平。在x$samples中，程序会自动计算每个样品的文库大小（即样品的总序列计数），归一化系数会被预先设置为1。 为了方便阅读，我们从DGEList对象x的列名中删去了GEO样品ID（GSM*）。

samplenames<-substring(colnames(x),12,nchar(colnames(x)))
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

运行结果

4.组织基因注释

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中选择某一种作为检索键。

与任何基因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第一次出现的信息。

在此例子中，注释与数据对象中的基因顺序是相同的。如果由于缺失和／或重新排列基因ID导致其顺序不一致，我们可以用match函数来正确排序基因。然后，我们将基因注释的数据框添加到DGEList对象，数据的整合就完成了，此时的数据对象中含有原始计数数据以及相关的样品信息和基因注释。

geneid<-rownames(x)
genes<-select(Mus.musculus,keys=geneid,columns=
c("SYMBOL","TXCHROM"),
              keytype = "ENTREZID")
head(genes)
genes<-genes[!duplicated(genes$ENTREZID),]
x$genes<-genes
x

五、数据预处理

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=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)
summary(lcpm)

结果展示

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

2、删除低表达基因

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

table(rowSums(x$counts==0)==9)

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

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

keep.exprs<-filterByExpr(x,group = group)
x<-x[keep.exprs,,keep.lib.sizes=FALSE]
dim(x)

此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有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")

运行结果

每个样本过滤前的原始数据（A）和过滤后（B）的数据的log-CPM值密度。竖直虚线标出了过滤步骤中所用阈值（相当于CPM值为约0.2）。

3、归一化基因表达分布

在样品制备或测序过程中，不具备生物学意义的外部因素会影响单个样品的表达。比如说，在实验中第一批制备的样品会总体上表达高于第二批制备的样品。差异表达分析假设所有样品表达值的范围和分布都应当相似。我们需要进行归一化来确保整个实验中每个样本的表达分布都相似。

密度图和箱线图等展示每个样品基因表达量分布的图表可以用于判断是否有样品和其他样品分布有差异。在此数据集中，所有样品的log-CPM分布都很相似。

尽管如此，我们依然需要使用edgeR中的calcNormFactors函数，用TMM(Robinson and Oshlack 2010)方法进行归一化。此处计算得到的归一化系数被用作文库大小的缩放系数。当我们使用DGEList对象时，这些归一化系数被自动存储在x$samples$norm.factors。对此数据集而言，TMM归一化的作用比较温和，这体现在所有的缩放因子都相对接近1。

x<-calcNormFactors(x,method="TMM")#TMM归一化处理
x$sample$norm.factors
#为了更好地展示出归一化的效果，我们复制了数据并进行了人工调整，
#使得第一个样品的计数减少到了其原始值的5%，而第二个样品增大到了5倍。
#要注意在实际的数据分析流程中，不应当进行这样的操作。
x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5

#下图显示了没有经过归一化的与经过了归一化的数据的表达分布
#其中归一化前不同样本的分布明显不同，而归一化后比较相似。
#此处，经过我们人工调整的第一个样品的TMM缩放系数0.06非常小
#而第二个样品的缩放系数6.08很大，它们都并不接近1。
par(mfrow=c(1,2))
lcpm<-cpm(x2,log=TRUE)
boxplot(lcpm,las=2,col=col,main="")
title(main="A.Example:Unnormalised data",ylab="Log-cpm")
x2<-calcNormFactors(x2)
x2$samples$norm.factors

lcpm<-cpm(x2,log=TRUE)
boxplot(lcpm,las=2,col=col,main="")
title(main="B.Example:Normalised data",ylab="Log-cpm")

运行结果

样例数据：log-CPM值的箱线图展示了未经归一化的数据（A）及归一化后的数据（B）中每个样本的表达分布。数据集经过调整，样本1和2中的表达计数分别被缩放到其原始值的5%和500%。

4、对样本的无监督聚类

在我们看来，用于检查基因表达分析的最重要的探索性图表之一便是MDS图，或类似的图。这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。理想情况下，样本会在各个实验组内很好的聚类，且我们可以鉴别出远离所属组的样本，并追踪误差或额外方差的来源。如果存在的话，技术重复应当互相非常接近。

这样的图可以用limma中的plotMDS函数绘制。第一个维度表示能够最好地分离样品且解释最大比例的方差的领先倍数变化（leading-fold-change），而后续的维度的影响更小，并与之前的维度正交。当实验设计涉及到多个因子时，建议在多个维度上检查每个因子。如果在其中一些维度上样本可按照某因子聚类，这说明该因子对于表达差异有影响，在线性模型中应当将其包括进去。反之，没有或者仅有微小影响的因子在下游分析时应当被剔除。

在这个数据集中，可以看出样本在维度1和2能很好地按照实验分组聚类，随后在维度3按照测序泳道（样品批次）分离（如下图所示）。由于第一维度解释了数据中最大比例的方差，我们会发现当关注更高维度时，维度上的取值范围会变小。

尽管所有样本都按组聚类，在维度1上最大的转录差异出现在basal和LP以及basal和ML之间。因此，预期在basal样品与其他之间的成对比较中能够得到大量的DE基因，而在ML和LP之间的比较中得到的DE基因数量略少。在其他的数据集中，不按照实验组聚类的样本可能在下游分析中只表现出较小的或不表现出差异表达。

为绘制MDS图，我们为不同的因子设立不同的配色。维度1和2以细胞类型上色，而维度3和4以测序泳道（批次）上色。

#绘制MDS图
lcpm<-cpm(x,log=TRUE)
par(mfrow=c(1,2))
col.group<-group
levels(col.group)<-brewer.pal(nlevels(col.group),"Set1")
col.group<-as.character(col.group)
col.lane<-lane
levels(col.lane)<-brewer.pal(nlevels(col.lane),"Set2")
col.lane<-as.character(col.lane)
plotMDS(lcpm,labels = group,col=col.group)
title(main = "A.Sample groups")
plotMDS(lcpm,labels = lane,col=col.lane,dim=c(3,4))
title(main = "B.Sequencing lanes")

运行结果

使用Glimma中的glMDSplot函数可以得到交互式MDS图。其中的glMDSPlot函数可生成一个html网页（如果设置launch=TRUE参数，将会在生成后直接在浏览器中打开），其左侧面板含有一张MDS图，而右侧面板包含一张展示了各个维度所解释的方差比例的柱形图。点击柱形图中的柱可切换MDS图绘制时所使用的维度，且将鼠标悬浮于单个点上可显示相应的样本标签。也可切换配色方案，以突显不同细胞类型或测序泳道（批次）。

glMDSPlot(lcpm, labels=paste(group, lane, sep="_"), 
          groups=x$samples[,c(2,5)], launch=TRUE)

运行结果

小结：今日告一段落，明日继续~