使用limma、Glimma和edgeR，RNA-seq数据分析易如反掌
使用limma、Glimma和edgeR，RNA-seq数据分析易如反掌
23 October 2018
Contents
1摘要
简单且高效地分析RNA测序数据的能力正是Bioconductor的核心优点。 RNA-seq分析一般从基因水平的序列计数开始，涉及到数据预处理，探索性数据分析，差别表达检验以及通路分析，获得的结果可用于指导进一步实验和验证研究。 在这篇工做流程文章中，咱们经过分析来自小鼠乳腺的RNA测序数据，示范了如何使用流行的edgeR包载入、整理、过滤和归一化数据，而后用limma包的voom方法、线性模型和经验贝叶斯调节（empirical Bayes moderation）来评估差别表达并进行基因集检验。经过使用Glimma包，此流程获得了增进，实现告终果的互动探索，使用户得以查看单个样本与基因。 这三个软件包提供的完整分析突出了研究人员可使用Bioconductor轻松地从RNA测序实验的原始计数揭示生物学意义。ios
2背景介绍
RNA测序（RNA-seq）已成为基因表达研究的主要技术。其中，基因组规模的多条件基因差别表达检测是研究者最常探究的问题之一。对于RNA-seq数据，来自Bioconductor项目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用于处理此问题的完善的统计学方法。git
在这篇文章中，咱们描述了一个用于分析RNA-seq数据的edgeR - limma工做流程，使用基因水平计数做为输入，通过预处理和探索性数据处理，而后获得差别表达（DE）基因和基因表达特征（gene signatures）的列表。Glimma包(Su and Ritchie 2016)提供的交互式图表能够同时呈现总体样本和单个基因水平的数据，使得咱们相对静态的R图表而言，能够探索更多的细节。数据库
此工做流程中所分析的实验来自Sheridan等（2015）(Sheridan et al. 2015)，由三个细胞群组成，即基底（basal）、管腔祖细胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。细胞群皆分选自雌性处女小鼠的乳腺，每种都设三个生物学重复。RNA样品分三个批次使用Illumina HiSeq 2000进行测序，获得长为100碱基对的单端序列片断。 本文所描述的分析假设从RNA-seq实验得到的序列片断已经与适当的参考基因组比对，并已经在基因水平上对序列进行了统计计数。在本文条件下，使用Rsubread包提供的基于R的流程将序列片断与小鼠参考基因组（mm20）比对（具体而言，先使用align函数(Liao, Smyth, and Shi 2013)，而后使用featureCounts (Liao, Smyth, and Shi 2014)进行基因水平的总结，利用其内置的mm10基于RefSeq的注释），json
这些样本的计数数据能够从Gene Expression Omnibus (GEO)数据库 http://www.ncbi.nlm.nih.gov/geo/使用GEO序列登记号GSE63310下载。更多关于实验设计和样品制备的信息也能够在GEO使用该登记号查看。api
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样品（请查看下方相关文件名）。markdown
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），九列分别对应此实验中的每一个样品。session
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关联于染色体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那样考虑到基因长度区别的影响，但在咱们的分析中常常会用到它们。虽然也可使用RPKM和FPKM，但CPM和log-CPM只使用计数矩阵便可计算，且已足以知足咱们所感兴趣的比较的须要。假设不一样条件之间剪接异构体（isoform）的使用没有差异，差别表达分析研究同一基因在不一样条件下的表达差别，而不是比较多个基因之间的表达或测定绝对表达量。换而言之，基因长度在咱们感兴趣的比较中始终不变，且任何观测到的差别是来自于条件的变化而不是基因长度的变化。
在此处，使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。其中，在进行对数转换前会给全部基因的计数加上2，以免对零取对数，且可减少低表达基因之间的差别。 若是能够提供基因长度信息，可使用edgeR中的rpkm函数计算RPKM值，就像计算CPM值那样简单。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
5.2删除低表达基因
全部数据集中都混有表达的基因与不表达的基因。尽管咱们想要检测在一种条件中表达但再另外一种条件中不表达的基因，也有一些基因在全部样品中都不表达。实际上，这个数据集中19%的基因在全部九个样品中的计数都是零。
table(rowSums(x$counts==0)==9)
## ## FALSE TRUE ## 22026 5153
咱们应当忽略在任何条件都没有表达到具备生物学意义的程度的基因，这样能够减少基因的子集而保留感兴趣的基因，且能够减少下游在进行差别表达测试时的计算量。 经过检查log-CPM值，能够看出每一个样本中很大一部分基因都是不表达或低表达的（以下图中A部分所示）。在此使用CPM值为1做为标准（至关于log-CPM值为0），将表达量高于此阈值的基因看做表达，反之看做不表达。咱们只保留在至少一个组（或者至少在整个实验的三个样品中）表达的基因用于下游分析。
尽管任何合理的值都能用做表达的阈值，通常咱们在分析中使用CPM值为1，由于它对于大多数数据集都能很好地分隔表达的基因与不表达的基因。在这里，CPM值为1意味着若是一个基因在测序深度最低的样品中（JMS9-P8c, 序列数量约2千万）有至少20个计数，或者在测序深度最高的样品中（JMS8-3，序列数量约7.6千万）有至少76个计数，那么它是“表达的”。若是测序的序列片断是在外显子水平而不是基因水平统计的，且／或实验的测序深度很低，可能须要考虑更低的CPM阈值。
keep.exprs <- rowSums(cpm>1)>=3 x <- x[keep.exprs,, keep.lib.sizes=FALSE] dim(x)
## [1] 14165 9
使用这个标准，基因的数量减小到了开始时数量的大约一半（14165个基因，下图B部分）。须要注意的是，从整个DGEList对象中取子集不只会删除基因计数，同时也删除了相关的基因信息。下方给出的是绘图所用代码。
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.21), las=2, main="", xlab="") title(main="A. Raw data", xlab="Log-cpm") abline(v=0, 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, prior.count=2) plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.21), las=2, main="", xlab="") title(main="B. Filtered data", xlab="Log-cpm") abline(v=0, 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")
Figure 1: 每一个样本过滤前的原始数据（A）和过滤后（B）的数据的log-CPM值密度。竖直虚线标出了过滤步骤中所用阈值即log-CPM为0（至关于CPM值为1）。
5.3归一化基因表达分布
在样品制备或测序过程当中，不具有生物学意义的外部因素会影响单个样品的表达。好比说，在实验中第一批制备的样品会整体上表达高于第二批制备的样品。假设全部样品的表达值的范围和分布都应当类似，须要进行归一化来确保整个实验中每一个样本的表达分布都类似。
密度图和箱线图等展现每一个样品基因表达量分布的图表能够用于判断是否有样品和其余样品分布有差别。在此数据集中，全部样品的log-CPM分布都很类似（上图B部分）。
尽管如此，咱们依然须要使用edgeR中的calcNormFactors函数，用TMM(Robinson and Oshlack 2010)方法进行归一化。此处计算获得的归一化系数被用做文库大小的缩放系数。当咱们使用DGEList对象时，这些归一化系数被自动存储在x$samples$norm.factors。对此数据集而言，TMM归一化的做用比较温和，这体如今全部的缩放因子都相对接近1。
x <- calcNormFactors(x, method = "TMM") x$samples$norm.factors
## [1] 0.896 1.035 1.044 1.041 1.032 0.922 0.984 1.083 0.979
为了更好地可视化表现出归一化的影响，咱们复制了数据并进行了调整，使得第一个样品的计数减小到了其原始值的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.05很是小，而第二个样品的缩放系数6.13很大，它们都并不接近1。
par(mfrow=c(1,2)) lcpm <- cpm(x2, log=TRUE, prior.count=2) boxplot(lcpm, las=2, col=col, main="") title(main="A. Example: Unnormalised data",ylab="Log-cpm") x2 <- calcNormFactors(x2) x2$samples$norm.factors
## [1] 0.0547 6.1306 1.2293 1.1705 1.2149 1.0562 1.1459 1.2613 1.1170
lcpm <- cpm(x2, log=TRUE, prior.count=2) boxplot(lcpm, las=2, col=col, main="") title(main="B. Example: Normalised data",ylab="Log-cpm")
Figure 2: 样例数据：log-CPM值的箱线图展现了未经归一化的数据（A）及归一化后的数据（B）中每一个样本的表达分布。数据集通过调整，样本1和2中的表达计数分别被缩放到其原始值的5%和500%。
5.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使用以测序泳道（批次）定义的色彩组合进行检查。
lcpm <- cpm(x, log=TRUE, prior.count=2) 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")
Figure 3: log-CPM值在维度1和2的MDS图，以样品分组上色并标记（A）和维度3和4的MDS图，以测序道上色并标记（B）。图中的距离对应于最主要的倍数变化（fold change），默认状况下也就是前500个在每对样品之间差别最大的基因的平均（均方根）log2倍数变化。
做为另外一种选择，Glimma包也提供了便于探索多个维度的交互式MDS图。其中的glMDSPlot函数可生成一个html网页（若是设置launch=TRUE，将会在浏览器中打开），其左侧面板含有一张MDS图，而右侧面板包含一张展现了各个维度所解释的方差比例的柱形图。点击柱形图中的柱可切换MDS图绘制时所使用的维度，且将鼠标悬浮于单个点上可显示相应的样本标签。也可切换配色方案，以突显不一样细胞类型或测序泳道（批次）。此数据集的交互式MDS图能够从http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MDS-Plot.html看到。
glMDSPlot(lcpm, labels=paste(group, lane, sep="_"), groups=x$samples[,c(2,5)], launch=FALSE)
6差别表达分析
6.1建立设计矩阵和对比
在此研究中，咱们想知道哪些基因在咱们研究的三组细胞之间以不一样水平表达。在咱们的分析中，假设基础数据是正态分布的，为其拟合一个线性模型。在此以前，须要建立一个包含细胞类型以及测序泳道（批次）信息的设计矩阵。
design <- model.matrix(~0+group+lane) colnames(design) <- gsub("group", "", colnames(design)) design
## Basal LP ML laneL006 laneL008 ## 1 0 1 0 0 0 ## 2 0 0 1 0 0 ## 3 1 0 0 0 0 ## 4 1 0 0 1 0 ## 5 0 0 1 1 0 ## 6 0 1 0 1 0 ## 7 1 0 0 1 0 ## 8 0 0 1 0 1 ## 9 0 1 0 0 1 ## attr(,"assign") ## [1] 1 1 1 2 2 ## attr(,"contrasts") ## attr(,"contrasts")$group ## [1] "contr.treatment" ## ## attr(,"contrasts")$lane ## [1] "contr.treatment"
对于一个给定的实验，一般有几种等价的方法能够建立一个合适的设计矩阵。 好比说，~0+group+lane去除了第一个因子group的截距，但第二个因子lane的截距被保留。 此外也可使用~group+lane，来自group和lane的截距均被保留。 此处的关键是理解如何解释给定模型中估计获得的系数。 咱们在此分析中选取第一种模型，由于在没有group的截距的状况下能更直截了当地设定模型中的对比。用于细胞群之间成对比较的对比能够在limma中用makeContrasts函数设定。
contr.matrix <- makeContrasts( BasalvsLP = Basal-LP, BasalvsML = Basal - ML, LPvsML = LP - ML, levels = colnames(design)) contr.matrix
## Contrasts ## Levels BasalvsLP BasalvsML LPvsML ## Basal 1 1 0 ## LP -1 0 1 ## ML 0 -1 -1 ## laneL006 0 0 0 ## laneL008 0 0 0
limma的线性模型方法的核心优点之一即是其适应任意实验复杂程度的能力。简单的设计，好比此工做流程中关于细胞类型和批次的实验设计，直到更复杂的因子设计和含有交互做用项的模型，都可以被相对简单地处理。当实验或技术效应可被随机效应模型（random effect model）模拟时，limma中的另外一种可能性是使用duplicateCorrelation函数来估计交互做用，这须要在此函数以及lmFit的线性建模步骤均指定一个block参数。
6.2从表达计数数据中删除异方差
据显示对于RNA-seq计数数据而言，当使用原始计数或当其被转换为log-CPM值时，方差并不独立于均值(Law et al. 2014)。使用负二项分布来模拟计数的方法假设均值与方差间具备二次的关系。在limma中，假设log-CPM值符合正态分布，并使用由voom函数计算获得的精确权重来调整均值与方差的关系，从而对log-CPM值进行线性建模。
当操做DGEList对象时，voom从x中自动提取文库大小和归一化因子，以此将原始计数转换为log-CPM值。在voom中，对于log-CPM值额外的归一化能够经过设定normalize.method参数来进行。
下图左侧展现了这个数据集log-CPM值的均值-方差关系。一般而言，方差是测序实验中的技术差别和不一样细胞类型的重复样本之间的生物学差别的结合，而voom图会显示出一个在均值与方差之间递减的趋势。 生物学差别高的实验一般会有更平坦的趋势，其方差值在高表达处稳定。 生物学差别低的实验更倾向于急剧降低的趋势。
不只如此，voom图也提供了对于上游所进行的过滤水平的可视化检测。若是对于低表达基因的过滤不够充分，在图上表达低的一端，受到很是低的表达计数的影响，能够观察到方差水平的降低。若是观察到了这种状况，应当回到最初的过滤步骤并提升用于该数据集的表达阈值。
当前面观察的MDS图中具备明显的样本水平的差别时，能够用voomWithQualityWeights函数来同时合并样本水平的权重和voom(Liu et al. 2015)估算获得的丰度相关的权重。关于此种状况的例子参见Liu等(2016) (Liu et al. 2016)。
par(mfrow=c(1,2)) v <- voom(x, design, plot=TRUE) v
## An object of class "EList" ## $genes ## ENTREZID SYMBOL TXCHROM ## 1 497097 Xkr4 chr1 ## 6 27395 Mrpl15 chr1 ## 7 18777 Lypla1 chr1 ## 9 21399 Tcea1 chr1 ## 10 58175 Rgs20 chr1 ## 14160 more rows ... ## ## $targets ## files group lib.size norm.factors lane ## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 29409426 0.896 L004 ## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 36528591 1.035 L004 ## purep53 GSM1545538_purep53.txt Basal 59598629 1.044 L004 ## JMS8-2 GSM1545539_JMS8-2.txt Basal 53382070 1.041 L006 ## JMS8-3 GSM1545540_JMS8-3.txt ML 78175314 1.032 L006 ## JMS8-4 GSM1545541_JMS8-4.txt LP 55762781 0.922 L006 ## JMS8-5 GSM1545542_JMS8-5.txt Basal 54115150 0.984 L006 ## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 23043111 1.083 L008 ## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19525423 0.979 L008 ## ## $E ## Samples ## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c ## 497097 -4.29 -3.87 2.52 3.30 -4.48 -3.99 3.31 -3.20 -5.29 ## 27395 3.88 4.40 4.52 4.57 4.32 3.79 3.92 4.35 4.13 ## 18777 4.71 5.56 5.40 5.17 5.63 5.08 5.08 5.76 5.15 ## 21399 4.78 4.74 5.37 5.13 4.85 4.94 5.16 5.04 4.99 ## 58175 3.94 3.29 -1.77 -1.88 2.99 3.36 -2.11 3.14 3.52 ## 14160 more rows ... ## ## $weights ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] ## [1,] 1.18 1.18 20.53 20.98 1.77 1.22 21.13 1.18 1.18 ## [2,] 20.88 26.56 31.60 29.66 32.56 26.75 29.79 21.90 17.15 ## [3,] 28.00 33.70 34.85 34.46 35.15 33.55 34.52 31.44 25.23 ## [4,] 27.67 29.60 34.90 34.43 34.84 33.16 34.49 26.14 24.50 ## [5,] 19.74 18.66 3.18 2.63 24.19 24.01 2.65 13.15 14.35 ## 14160 more rows ... ## ## $design ## Basal LP ML laneL006 laneL008 ## 1 0 1 0 0 0 ## 2 0 0 1 0 0 ## 3 1 0 0 0 0 ## 4 1 0 0 1 0 ## 5 0 0 1 1 0 ## 6 0 1 0 1 0 ## 7 1 0 0 1 0 ## 8 0 0 1 0 1 ## 9 0 1 0 0 1 ## attr(,"assign") ## [1] 1 1 1 2 2 ## attr(,"contrasts") ## attr(,"contrasts")$group ## [1] "contr.treatment" ## ## attr(,"contrasts")$lane ## [1] "contr.treatment"
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
Figure 4: 图中绘制了每一个基因的均值（x轴）和方差（y轴），显示了在该数据上使用voom前它们之间的相关性（左），以及当运用voom的精确权重后这种趋势是如何消除的（右）。左侧的图是使用voom函数绘制的，它为进行log-CPM转换后的数据拟合线性模型从而提取残差方差。而后，对方差取平方根（或对标准差取平方根），并相对每一个基因的平均表达做图。均值经过平均计数加上0
5再进行log2转换计算获得。右侧的图使用plotSA绘制了log2残差标准差与log-CPM均值的关系。平均log2残差标准差由水平蓝线标出。在这两幅图中，每一个黑点表示一个基因，红线为对这些点的拟合。
值得注意的是，DGEList对象中存储的另外一个数据框，即基因和样本水平信息所存储之处，保留在了voom建立的EList对象v中。v$genes数据框等同于x$genes，v$targets等同于x$samples，而v$E中所储存的表达值相似于x$counts，尽管它进行了尺度转换。此外，voom的EList对象中还有一个精确权重的矩阵v$weights，而设计矩阵存储于v$design。
6.3拟合线性模型以进行比较
limma的线性建模使用lmFit和contrasts.fit函数进行，它们原先是为微阵列而设计的。这些函数不只能够用于微阵列数据，也能够用于RNA-seq数据。它们会单独为每一个基因的表达值拟合一个模型。而后，经过利用全部基因的信息来进行经验贝叶斯调整，这样能够得到更精确的基因水平的变异程度估计(Smyth 2004)。下一图为此模型的残差关于平均表达值的图。从图中能够看出，方差再也不与表达水平均值相关。
6.4检查DE基因数量
为快速查看差别表达水平，显著上调或下调的基因能够汇总到一个表格中。显著性的判断使用校订p值阈值的默认值5%。在basal与LP的表达水平之间的比较中，发现了4127个在basal中相较于LP下调的基因和4298个在basal中相较于LP上调的基因，即共8425个DE基因。在basal和ML之间发现了一共8510个DE基因（4338个下调基因和4172个上调基因），而在LP和ML中发现了一共5340个DE基因（2895个下调基因和2445个上调基因）。在包括basal细胞类型的比较中皆找到了大量的DE基因，这与咱们在MDS图中观察到的结果相吻合。
summary(decideTests(efit))
## BasalvsLP BasalvsML LPvsML ## Down 4127 4338 2895 ## NotSig 5740 5655 8825 ## Up 4298 4172 2445
一些研究中不只仅须要使用校订p值阈值，更为严格定义的显著性可能须要差别倍数的对数（log-FCs）也高于某个最小值。treat方法(McCarthy and Smyth 2009)能够按照对最小log-FC值的要求，使用通过经验贝叶斯调整的t统计值计算p值。当咱们的检验要求基因的log-FC显著大于1（等同于在本来的尺度上不一样细胞类型之间差两倍）时，差别表达基因的数量获得了降低，basal与LP相比只有3135个DE基因，basal与ML相比只有3270个DE基因，LP与ML相比只有385个DE基因。
tfit <- treat(vfit, lfc=1) dt <- decideTests(tfit) summary(dt)
## BasalvsLP BasalvsML LPvsML ## Down 1417 1512 203 ## NotSig 11030 10895 13780 ## Up 1718 1758 182
在多种比较中皆差别表达的基因能够从decideTests的结果中提取，其中的0表明不差别表达的基因，1表明上调的基因，-1表明下调的基因。共有2409个基因在basal和LP以及basal和ML的比较中都差别表达，其中的20个于下方列出。write.fit函数可用于将三个比较的结果提取并写入到单个输出文件。
de.common <- which(dt[,1]!=0 & dt[,2]!=0) length(de.common)
## [1] 2409
head(tfit$genes$SYMBOL[de.common], n=20)
## [1] "Xkr4" "Rgs20" "Cpa6" "Sulf1" "Eya1" ## [6] "Msc" "Sbspon" "Pi15" "Crispld1" "Kcnq5" ## [11] "Ptpn18" "Arhgef4" "2010300C02Rik" "Aff3" "Npas2" ## [16] "Tbc1d8" "Creg2" "Il1r1" "Il18r1" "Il18rap"
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
Figure 5: 韦恩图展现了仅basal和LP（左）、仅basal和ML（右）的对比的DE基因数量，还有两种对比中共同的DE基因数量（中）。在任何对比中均不差别表达的基因数量标于右下。
write.fit(tfit, dt, file="results.txt")
6.5从上到下检查单个DE基因
使用topTreat函数能够列举出使用treat获得的结果中靠前的DE基因（对于eBayes的结果可使用topTable函数）。默认状况下，topTreat将基因按照校订p值从小到大排列，并为每一个基因给出相关的基因信息、log-FC、平均log-CPM、校订t值、原始及通过多重假设检验校订的p值。列出前多少个基因的数量可由用户指定，若是设为n=Inf则会包括全部的基因。基因Cldn7和Rasef在basal与LP和basal于ML的比较中都位于DE基因的前几名。
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf) basal.vs.ml <- topTreat(tfit, coef=2, n=Inf) head(basal.vs.lp)
## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 12759 12759 Clu chr14 -5.44 8.86 -33.4 3.99e-10 2.7e-06 ## 53624 53624 Cldn7 chr11 -5.51 6.30 -32.9 4.50e-10 2.7e-06 ## 242505 242505 Rasef chr4 -5.92 5.12 -31.8 6.06e-10 2.7e-06 ## 67451 67451 Pkp2 chr16 -5.72 4.42 -30.7 8.01e-10 2.7e-06 ## 228543 228543 Rhov chr2 -6.25 5.49 -29.5 1.11e-09 2.7e-06 ## 70350 70350 Basp1 chr15 -6.07 5.25 -28.6 1.38e-09 2.7e-06
head(basal.vs.ml)
## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 242505 242505 Rasef chr4 -6.51 5.12 -35.5 2.57e-10 1.92e-06 ## 53624 53624 Cldn7 chr11 -5.47 6.30 -32.5 4.98e-10 1.92e-06 ## 12521 12521 Cd82 chr2 -4.67 7.07 -31.8 5.80e-10 1.92e-06 ## 71740 71740 Nectin4 chr1 -5.56 5.17 -31.3 6.76e-10 1.92e-06 ## 20661 20661 Sort1 chr3 -4.91 6.71 -31.2 6.76e-10 1.92e-06 ## 15375 15375 Foxa1 chr12 -5.75 5.63 -28.3 1.49e-09 2.28e-06
6.6差别表达结果的实用图形表示
为可视化地总结全部基因的结果，可以使用plotMD函数绘制均值-差别（MD）图，其中展现了线性模型拟合所获得的log-FC与log-CPM平均值间的关系，而差别表达的基因会被重点标出。
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], xlim=c(-8,13))
Glimma的glMDPlot函数提供了交互式的均值-差别图，拓展了这种图表的功能性。 此函数的输出为一个html页面，左侧面板为结果的总结性图表（与plotMD的输出相似），右侧面板包含各个样本的log-CPM值，下方为结果的表格。 这种交互式展现容许用户使用基因名搜索特定基因，而这在R统计图中是作不到的。 glMDPlot函数不只限于均值-差别图，其默认版本容许读入数据框，而用户能够选择在左侧面板绘图所用的列。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1], side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)
使用Glimma生成的均值-差别图。左侧面板显示了总结性数据（log-FC与log-CPM值的关系），其中选中的基因在每一个样本中的数值显示于右侧面板。下方为结果的表格，其搜索框使用户得以使用可行的注释信息查找某个特定基因，如基因符号Clu。
上方指令生成的均值-差别图能够在线访问（详见http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MD-Plot.html）。 Glimma提供的交互性使得单个图形窗口内可以呈现出额外的信息。 Glimma是以R和Javascript实现的，使用R代码生成数据，并在以后使用Javascript库D3（https://d3js.org）转换为图形，使用Bootstrap库处理界面并生成互动性可搜索的表格的数据表。 这使得图表能够在任何现代的浏览器中查看，对于从Rmarkdown分析报告中将其做为关联文件而附加而言十分方便。
前文所展现的图表中，一些展现了在任意一个条件下表达的全部基因（好比共同DE基因的韦恩图或均值-差别图），而另外一些展现单独的基因（交互性均值-差别图右边面板中所展现的log-CPM值）。而热图使用户得以查看一部分基因的表达。这对于查看单个组或样本的表达颇有用，而不至于在关注于单个基因时失去对于研究总体的注意，也不会形成因为上千个基因所取平均值而致使的失去分辨率。
使用gplots包的heatmap.2函数，咱们为basal与LP的对照中前100个DE基因（按调整p值排序）绘制了一幅热图。热图中正确地将样本按照细胞类型聚类，并从新排序了基因，造成了表达类似的块状。从热图中，咱们观察到对于basal与LP之间的前100个DE基因，ML和LP样本的表达很是类似。
library(gplots) basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100] i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes) mycol <- colorpanel(1000,"blue","white","red") heatmap.2(lcpm[i,], scale="row", labRow=v$genes$SYMBOL[i], labCol=group, col=mycol, trace="none", density.info="none", margin=c(8,6), lhei=c(2,10), dendrogram="column")
Figure 6: 在basal和LP的对比中前100个DE基因log-CPM值的热图。通过缩放调整后，每一个基因（每行）的表达均值为0，而且标准差为1。给定基因相对高表达的样本被标记为红色，相对低表达的样本被标记为蓝色。浅色和白色表明中等表达水平的基因。样本和基因已经过分层聚类的方法从新排序。图中显示有样本聚类的树状图。
7使用camera的基因集检验
在这次分析的最后，咱们要进行一些基因集检验。为此，咱们将camera方法(Wu and Smyth 2012)应用于Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中适应小鼠的c2基因表达特征，这可从http://bioinf.wehi.edu.au/software/MSigDB/以RData对象格式获取。 此外，对于人类和小鼠，来自MSigDB的其余有用的基因集也可今后网站获取，好比表明性（hallmark）基因集。C2基因集的内容收集自在线数据库、出版物以及该领域专家，而标志基因集的内容来自MSigDB，从而对于生物状态或过程的定义足够明确。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123")) idx <- ids2indices(Mm.c2,id=rownames(v)) cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1]) head(cam.BasalvsLP,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_STEM_CELL_UP 739 Up 1.13e-18 5.36e-15 ## LIM_MAMMARY_STEM_CELL_DN 630 Down 1.57e-15 3.71e-12 ## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 163 Up 1.44e-13 2.26e-10 ## SOTIRIOU_BREAST_CANCER_GRADE_1_VS_3_UP 183 Up 2.18e-13 2.58e-10 ## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 87 Down 6.73e-13 6.36e-10
cam.BasalvsML <- camera(v,idx,design,contrast=contr.matrix[,2]) head(cam.BasalvsML,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_STEM_CELL_UP 739 Up 5.09e-23 2.40e-19 ## LIM_MAMMARY_STEM_CELL_DN 630 Down 5.13e-19 1.21e-15 ## LIM_MAMMARY_LUMINAL_MATURE_DN 166 Up 8.88e-16 1.40e-12 ## LIM_MAMMARY_LUMINAL_MATURE_UP 180 Down 6.29e-13 7.43e-10 ## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 163 Up 1.68e-12 1.59e-09
cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3]) head(cam.LPvsML,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_LUMINAL_MATURE_UP 180 Down 8.50e-14 3.40e-10 ## LIM_MAMMARY_LUMINAL_MATURE_DN 166 Up 1.44e-13 3.40e-10 ## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 87 Up 3.84e-11 6.05e-08 ## REACTOME_RESPIRATORY_ELECTRON_TRANSPORT 91 Down 2.66e-08 3.14e-05 ## NABA_CORE_MATRISOME 222 Up 4.43e-08 4.19e-05
camera函数经过比较假设检验来评估一个给定基因集中的基因是否相对于不在集内的基于是言在差别表达基因的排序中更靠前。 它使用limma的线性模型框架，并同时采用设计矩阵和对比矩阵（若是有），且在测试的过程当中会适应来自voom的观测水平权重。 在经过方差膨胀因子（variance inflation factor）调整获得的基因集检验统计值的方差后，取决于基因间相关性（默认设定为0.01，但也可经过数据估计）和基因集的规模，将会返回p值，且根据多重假设检验进行了校订。
此实验是等同于Lim等人(2010)(Lim et al. 2010)的RNA-seq，而他们使用Illumina微阵列分析了相同的分选细胞群，所以当咱们看到该早期文献中的基因表达特征出如今每种对比的列表顶部时大可放心。在LP和ML的对比中，咱们为Lim等人（2010）的成熟管腔基因集（上调及下调）绘制了条码图（barcodeplot）。须要注意的是，因为咱们的对比是将LP与ML相比而不是相反，这些基因集的方向在咱们的数据集中是反过来的（若是将对比反过来，基因集的方向将会与对比一致）。
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
Figure 7: LIM_MAMMARY_LUMINAL_MATURE_UP （红色条形，图表上方）和LIM_MAMMARY_LUMINAL_MATURE_DN（蓝色条形，图表下方）基因集在LP和ML的对比中的条码图，每一个基因集都有一条富集线展现了竖直条形在图表每部分的相对富集程度。Lim等人的实验(Lim et al
2010)很是相似于咱们的，用了相同的分选方式来获取不一样的细胞群，只是他们使用的是微阵列而不是RNA-seq来测定基因表达。须要注意的是，上调基因集发生下调而下调基因集发生上调的逆相关性来自于对比的设定方式（LP相比于ML），若是将其对调，方向性将会吻合。
limma也有其余的基因集检验，好比mroast(Wu et al. 2010)的自包含检验。虽然camera很是适合检验基因集的大型数据库并观察其中哪些相对于其余的在排序上位次更高（如前文所示），自包含检验更善于集中检验一个或少个选中的集合是否自己差别表达。换句话说，camera更适用于搜寻具备意义的基因集，而mroast测试的是已经肯定有意义的基因集的显著性。
8使用到的软件和代码
此RNA-seq工做流程使用了Bioconductor项目3.4版本中的多个软件包，运行于R 3.3.0或更高版本。除了本文中重点提到的软件（limma、Glimma以及edgeR），亦须要一些其余软件包，包括gplots和RColorBrewer还有基因注释包Mus.musculus。 此文档使用knitr编译。全部用到的包的版本号以下所示。进行此分析所需代码都可在Bioconductor工做流程包RNAseq123从http://www.bioconductor.org/help/workflows/RNAseq123/获取。
sessionInfo()
## R Under development (unstable) (2018-10-16 r75448) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 18.04.1 LTS ## ## Matrix products: default ## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so ## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 ## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C ## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] parallel stats4 stats graphics grDevices utils datasets methods ## [9] base ## ## other attached packages: ## [1] gplots_3.0.1 RColorBrewer_1.1-2 ## [3] Mus.musculus_1.3.1 TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.4 ## [5] org.Mm.eg.db_3.7.0 GO.db_3.7.0 ## [7] OrganismDbi_1.25.0 GenomicFeatures_1.35.0 ## [9] GenomicRanges_1.35.0 GenomeInfoDb_1.19.0 ## [11] AnnotationDbi_1.45.0 IRanges_2.17.0 ## [13] S4Vectors_0.21.0 Biobase_2.43.0 ## [15] BiocGenerics_0.29.1 edgeR_3.25.0 ## [17] Glimma_1.11.0 limma_3.39.1 ## [19] knitr_1.20 BiocStyle_2.11.0 ## ## loaded via a namespace (and not attached): ## [1] httr_1.3.1 bit64_0.9-7 jsonlite_1.5 ## [4] R.utils_2.7.0 gtools_3.8.1 assertthat_0.2.0 ## [7] BiocManager_1.30.3 highr_0.7 RBGL_1.59.0 ## [10] blob_1.1.1 GenomeInfoDbData_1.2.0 Rsamtools_1.35.0 ## [13] yaml_2.2.0 progress_1.2.0 RSQLite_2.1.1 ## [16] backports_1.1.2 lattice_0.20-35 digest_0.6.18 ## [19] XVector_0.23.0 htmltools_0.3.6 Matrix_1.2-15 ## [22] R.oo_1.22.0 XML_3.98-1.16 pkgconfig_2.0.2 ## [25] biomaRt_2.39.0 bookdown_0.7 zlibbioc_1.29.0 ## [28] gdata_2.18.0 BiocParallel_1.17.0 SummarizedExperiment_1.13.0 ## [31] magrittr_1.5 crayon_1.3.4 memoise_1.1.0 ## [34] evaluate_0.12 R.methodsS3_1.7.1 graph_1.61.0 ## [37] tools_3.6.0 prettyunits_1.0.2 hms_0.4.2 ## [40] matrixStats_0.54.0 stringr_1.3.1 locfit_1.5-9.1 ## [43] DelayedArray_0.9.0 Biostrings_2.51.0 compiler_3.6.0 ## [46] caTools_1.17.1.1 rlang_0.3.0.1 grid_3.6.0 ## [49] RCurl_1.95-4.11 bitops_1.0-6 rmarkdown_1.10 ## [52] DBI_1.0.0 R6_2.3.0 GenomicAlignments_1.19.0 ## [55] rtracklayer_1.43.0 bit_1.1-14 rprojroot_1.3-2 ## [58] KernSmooth_2.23-15 stringi_1.2.4 Rcpp_0.12.19 ## [61] xfun_0.4
参考文献
Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https://bioconductor.org/packages/release/data/annotation/html/Homo.sapiens.html.
———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https://bioconductor.org/packages/release/data/annotation/html/Mus.musculus.html.
Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber. 2005. “BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis.” Bioinformatics 21:3439–40.
Durinck, S., P. Spellman, E. Birney, and W. Huber. 2009. “Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.” Nature Protocols 4:1184–91.
Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2):115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.” Genome Biology 15:R29.
Liao, Y., G. K. Smyth, and W. Shi. 2013. “The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote.” Nucleic Acids Res 41 (10):e108.
———. 2014. “featureCounts: an Efficient General-Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7):923–30.
Lim, E., D. Wu, B. Pal, T. Bouras, M. L. Asselin-Labat, F. Vaillant, H. Yagita, G. J. Lindeman, G. K. Smyth, and J. E. Visvader. 2010. “Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways.” Breast Cancer Research 12 (2):R21.
Liu, R., K. Chen, N. Jansz, M. E. Blewitt, and M. E. Ritchie. 2016. “Transcriptional Profiling of the Epigenetic Regulator Smchd1.” Genomics Data 7:144–7.
Liu, R., A. Z. Holik, S. Su, N. Jansz, K. Chen, H. S. Leong, M. E. Blewitt, M. L. Asselin-Labat, G. K. Smyth, and M. E. Ritchie. 2015. “Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses.” Nucleic Acids Res 43:e97.
McCarthy, D. J., and G. K. Smyth. 2009. “Testing significance relative to a fold-change threshold is a TREAT.” Bioinformatics 25:765–71.
Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. “limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Res 43 (7):e47.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics26:139–40.
Robinson, M. D., and A. Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-seq data.” Genome Biology 11:R25.
Sheridan, J. M., M. E. Ritchie, S. A. Best, K. Jiang, T. J. Beck, F. Vaillant, K. Liu, et al. 2015. “A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1.” BMC Cancer 15 (1). BioMed Central:221.
Smyth, G. K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Stat Appl Genet Mol Biol 3 (1):Article 3.
Su, S., and M. E. Ritchie. 2016. Glimma: Interactive HTML graphics for RNA-seq data. https://bioconductor.org/packages/devel/bioc/html/Glimma.html.
Subramanian, A., P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc Natl Acad Sci U S A 102 (43):15545–50.
Wu, D., E. Lim, F. Vaillant, M. L. Asselin-Labat, J. E. Visvader, and G. K. Smyth. 2010. “ROAST: rotation gene set tests for complex microarray experiments.” Bioinformatics 26 (17):2176–82.
Wu, D., and G. K. Smyth. 2012. “Camera: a competitive gene set test accounting for inter-gene correlation.” Nucleic Acids Res 40 (17):e133.
- 1. 直通车上10分 易如反掌
- 2. 生物信息学入门 使用 RNAseq counts数据进行差异表达分析（DEG）——edgeR 算法 数据 代码 结果解读
- 3. 如何使用反射来分析PHP的数据结构
- 4. Python探索性数据分析，这样才容易掌握！
- 5. cgmodel简易数据分析
- 6. .Net反射机制分析和使用
- 7. RNAseq数据处理实验记录ydy_sir
- 8. 如何分析用户反馈数据? | 下
- 9. 数据分析和数据分析师
- 10. 如何才能更好使用数据分析工具分析数据？
- 更多相关文章...
- • XSD 如何使用? - XML Schema 教程
- • 使用netwox构造IP数据包 - TCP/IP教程
- • Flink 数据传输及反压详解
- • TiDB 在摩拜单车在线数据业务的应用和实践
-
每一个你不满意的现在，都有一个你没有努力的曾经。