From 8f869f14ca1b35a7ab2e55418fe6b0162e56c9d0 Mon Sep 17 00:00:00 2001 From: restore1997 <47097886+restore1997@users.noreply.github.com> Date: Wed, 31 Jul 2019 19:56:01 +0800 Subject: [PATCH 1/3] Update 10RnaSeqGeneEdgeRQL_en.Rmd --- BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd | 941 ++++++++++++++++++- 1 file changed, 940 insertions(+), 1 deletion(-) diff --git a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd index e93d3c6..fdbcdba 100644 --- a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd +++ b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd @@ -2,7 +2,946 @@ title: "From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline" author: - name: Yunshun Chen - affiliation: The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia + affiliation: The Wa--- +题目：「从序列到基因到通路：使用Rsubread包和edgeR包拟似然管道进行RNA-Seq差异表达分析」 +作者： + - 姓名： Yunshun Chen + 附属机构： 沃尔特和伊丽莎霍尔医学研究所，1G Royal Parade，帕克维尔，维多利亚州 3052，墨尔本，澳大利亚; 医学生物学系，墨尔本大学，帕克维尔，维多利亚州 3010，墨尔本，澳大利亚 + - 姓名： Aaron T. L. Lun + 附属机构： 英国剑桥癌症研究所，剑桥大学，李嘉诚中心，罗宾逊路，CB20RE，剑桥，英国 + - 姓名： Gordon K. Smyth + 附属机构： 沃尔特和伊丽莎霍尔医学研究所，1G Royal Parade，帕克维尔，维多利亚州 3052，墨尔本，澳大利亚;数学及统计学系，墨尔本大学，帕克维尔，维多利亚州 3010，墨尔本，澳大利亚 +日期： 2018 年 10 月 28 号 +简介： > + %\简介索引条目{From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline} + %\简介工具{knitr::rmarkdown} + %\简介编码{UTF-8} +参考书目： ref.bib +输出： + 文件类型：：html文件： + 图片标题： 含有 +--- +**R 版本**： `r R.version.string` + +**Bioconductor 版本**： `r BiocManager::version()` + +**包**： `r packageVersion("RnaSeqGeneEdgeRQL")` + +摘要 {#abstract .unnumbered} +======== + +RNA 序列分析 (RNA-seq)广泛应用于在基因表达分析中。RNA-seq分析最常见的目的之一是识别在两种或两种以上生物学条件下差异表达(DE)的基因或分子途径。通过提供小鼠乳腺上皮细胞亚群完整的RNA-seq实验分析，本文演示了从RNA-seq数据中检测DE基因和通路的计算流程。工作流使用来自开源Bioconductor项目的R包，涵盖了分析管道的所有步骤，包括读取序列的比对、数据探索、差异表达分析、可视化和通路分析。使用Rsubread包进行序列比对和计数量化，使用edgeR进行数据分析。使用edgeR包拟似然函数进行差异表达分析。 + + +介绍 {#introduction .unnumbered} +============ + +RNA 序列分析 (RNA-seq) 大量应用于分析生物系统的转录活性。 +RNA-seq分析最常见的目的之一是识别在两种或两种以上生物学条件下差异表达(DE)的基因或分子途径。 +表达的变化可能与生物学上的差异有关，为进一步研究潜在的作用机制提供了途径。 + +通过使用 Bioconductor 包，从原始序列到差异表达和通路分析，本文为提供了一个完整的分析 RNA-seq 的工作流程[@huber2015orchestrating]。 +本文对 RNA-seq 数据进行了全面的分析，这些数据是为了研究怀孕和哺乳对小鼠乳腺腔细胞谱系的影响而收集的[@fu2015egfmcl1]。 +管道使用 Rsubread 包 [@liao2013subread] 读取序列映射并分配给基因，使用 edgeR 包 [@robinson2010edger]进行数据分析 + +此工作流程的期刊文章版本，包括裁判报告和回复，可在以下网址中查看： [F1000Research](http://f1000research.com/articles/5-1438) [@chen2016quasiworkflow]. + +RNA-seq 分析有很多步骤，包括序列比对、序列汇总、差异表达和通路分析。 +这里我们使用 Subread aligner [@liao2013subread] 进行比对，使用 featureCounts [@liao2014featurecounts] 分配序列到基因。 +使用 Rsubread 包不仅快速、高效，还具有原生实现的优势。 +这意味着所有的分析都能在 R 环境中有序控制。 + +使用 edgeR 包的拟似然函数（edgeR-quasi）进行差异分析。 +这种统计方法使用负二项广义线性模型 [@mccarthy2012edgerglm] 涉及F检验s，而非似然比检验 [@lund2012quasiseq]。 +该方法比其他基于负二项的管道具有更严格的出错率控制, 包括传统的 edgeR 包管道 [@robinson2008smallsample; @robinson2007moderated; @mccarthy2012edgerglm] 或者 DESeq2 包 [@love2014deseq2]。 +edgeR-quasi 管道与 QuasiSeq 包 [@lund2012quasiseq] 有相似的统计学方法，其在第三方比较中表现良好（third-party comparisons） [@burden2014error]。 +与 QuasiSeq 包相比,edgeR 包提供了更快、更细致的统计方法 [@lun2016delicious]。 +RNA-seq 管道的 limma 包也提供了极佳的错误率控制 [@law2014voom; @ritchie2015limma]。 +对于大型数据集，鉴于其速度和灵活性，推荐使用 limma 管道，edgeR-quasi 管道在低计数情况下可以提供更好的性能 [@lun2014denovo; @lun2015diffhic]。 +对于这里的数据分析, edgeR-quasi ， limma-voom 和 limma-trend 管道同样适合，输出的结果相似。 + +本文阐述的分析方法可以应用于任何包含一些复制的？RNA-seq研究，但它特别适用于多处理因素和少量生物复制的设计实验。 +该方法假设在两种或两种以上的处理条件下从感兴趣的细胞中提取RNA样本，RNA-seq分析已应用于每个RNA样本，并且至少有一种处理条件存在独立的生物复制。 +本流程中 Rsubread 部分使用原始序列的FASTQ文件作为输入，edgeR 部分使用genewise？序列计数的矩阵作为输入。 + +```{r setup, echo=FALSE, results="hide", warning=FALSE} +knitr::opts_chunk$set(prompt = TRUE, comment = NA) +suppressPackageStartupMessages({library(RnaSeqGeneEdgeRQL)}) +``` + + +生物实验描述 {#description-of-the-biological-experiment .unnumbered}？ +======================================== + +本流程描述了一个完整的RNA-seq 生信分析，该数据可以从GEO数据库的GSE60450系列获取。 +该 RNA-seq 数据用于研究小鼠乳腺腔细胞谱系，特别是在妊娠和哺乳时谱系成员的表达谱如何变化[@fu2015egfmcl1]。 +特别说明，本研究监测了未妊娠、妊娠和哺乳期老鼠乳腺基础干细胞富集细胞(B)和固定腔细胞(L)的表达情况。 +所以共有6组 RNA 样本，每组对应特定的细胞类型和老鼠情况的组合。 +每组收集两个生物学副本。 + +本研究使用 Illumina Hiseq 测序仪对每个样本产生约3000万个100bp的单端读数。 +使用Subread 版本 1.4.4()将读取结果与小鼠mm10基因组进行比对，使用 featureCounts 通过 RefSeq 基因注释将读取结果匹配到Entrez Genes。 +将包含原始序列读取结果的FASTQ文件保存到SRA数据库中，读取数据保存到GEO数据库。 + +本实验设计总结如下表所示，基底细胞型和腔细胞型分别缩写为“B”和“L”。 +每个RNA样本的GEO和SRA标识符也显示如下: +```{r sra} +library(RnaSeqGeneEdgeRQL) +targetsFile <- system.file("extdata", "targets.txt", + package="RnaSeqGeneEdgeRQL") +targets <- read.delim(targetsFile, stringsAsFactors=FALSE) +targets +``` +实验可以看作是六组的单向布局。 +为便于以后使用，我们将处理因子合并为一个单独的分组因子: + +```{r group} +group <- paste(targets$CellType, targets$Status, sep=".") +group <- factor(group) +table(group) +``` + +注意： +本研究在可控的实验条件下，将相同基因的小鼠仔细分类成不同细胞群。 +如下面的分析所示，这里相对较小的样本量(每组n=2)是由低的背景变异性和不同组的表达谱明显不同所证明的。 +具有更大变异性的研究(例如对人类受试者的研究)或具有更小效应的研究很可能需要更多的生物复制来获得可靠的结果。 + + +初步分析 {#preliminary-analysis .unnumbered} +==================== + +下载读取计数 {#downloading-the-read-counts .unnumbered} +--------------------------- + +希望重新生成本文中给出的分析的读者可以从GEO下载读取计数矩阵，或者从原始序列计数重新创建读取计数矩阵。 +我们将首先使用下载的计数矩阵进行分析。 +在本文的最后，我们将给出重新创建这个矩阵所需的R命令。 + +下面的命令为从GEO系列GSE60450下载genewise读取计数。 +压缩的以tab分隔的文本文件“GSE60450_Lactation-GenewiseCounts.txt.gz”将下载到R工作目录: + +```{r checkdownload, echo=FALSE, results="hide", message=FALSE} +if( !file.exists("GSE60450_Lactation-GenewiseCounts.txt.gz") ) { +FileURL <- paste( + "http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE60450", + "format=file", + "file=GSE60450_Lactation-GenewiseCounts.txt.gz", + sep="&") +download.file(FileURL, method="libcurl", "GSE60450_Lactation-GenewiseCounts.txt.gz") +} +``` + +```{r download, eval=FALSE} +FileURL <- paste( + "http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE60450", + "format=file", + "file=GSE60450_Lactation-GenewiseCounts.txt.gz", + sep="&") +download.file(FileURL, "GSE60450_Lactation-GenewiseCounts.txt.gz") +``` + +可以将这些计数放在R的数据框中： + +```{r readcounts} +GenewiseCounts <- read.delim("GSE60450_Lactation-GenewiseCounts.txt.gz", + row.names="EntrezGeneID") +colnames(GenewiseCounts) <- substring(colnames(GenewiseCounts),1,7) +dim(GenewiseCounts) +head(GenewiseCounts) +``` +GenewiseCounts 的行名是 Entrez 基因标识符（Entrez Gene Identifiers）。 +首列包含每个基因的长度，表示该基因外显子和超转录序列中碱基的总数。 +其余12列包含阅读计数,与“targets”行对应。 + + +edgeR 包将数据储存在简单的基于列表的数据对象中，称为“DGEList”。 +这个对象很容易使用，因为它可以像R中的普通列表一样操作，也可以像矩阵一样作为子集。 +DGEList 对象包含的主要内容是阅读计数的矩阵，在“数据框”格式里的样本信息和可选的基因注释。 +我们在 edgeR 包中使用DGEList功能j将计数输入到一个DGEList对象中。 + +```{r DGEList, message=FALSE} +library(edgeR) +y <- DGEList(GenewiseCounts[,-1], group=group, + genes=GenewiseCounts[,1,drop=FALSE]) +options(digits=3) +y$samples +``` + +增加基因注释 {#adding-gene-annotation .unnumbered} +---------------------- + +Entrez 基因标识符与NCBI数据库的基因信息相关联。 +可以使用 org.Mm.eg.db 包补充基因注释信息。 +在这里将基因名这一列添加到“y$gene”中。 + +```{r symbols, message=FALSE} +library(org.Mm.eg.db) +y$genes$Symbol <- mapIds(org.Mm.eg.db, rownames(y), + keytype="ENTREZID", column="SYMBOL") +head(y$genes) +``` +不再具有官方基因符号的Entrez b标识符将从分析中删除。 +整个“DGEList”对象，包括注释和计数，可以像矩阵一样按行做子集： + +```{r dropNAsymbols} +y <- y[!is.na(y$genes$Symbol), ] +dim(y) +``` + + +设计矩阵 {#design-matrix .unnumbered} +------------- + +在 edgeR 包中，需要特殊的设计矩阵来建立线性模型和差异表达分析。 +设计矩阵记录了每个样本的处理条件，并定义了线性模型中实验效果的参数化。 +本研究的实验设计可以看作是单向布局，设计矩阵可以简单直观地构造: + +```{r design} +design <- model.matrix(~0+group) +colnames(design) <- levels(group) +design +``` + +这个设计矩阵只是将每个组链接到属于它的示例。 +设计矩阵的每一行对应一个样本，而每一列表示对应于六个组中的一个组的系数。 + + +过滤和删除低计数 {#filtering-to-remove-low-counts .unnumbered} +------------------------------ + +所有文库中计数非常低的基因都应该在下游分析之前删除。 +这在生物学和统计学上都是合理的。 +从生物学的角度来看，一个基因在被翻译成蛋白质或被认为具有重要的生物学意义之前，必须在某个最低水平上表达。 +从统计学的角度来看，持续低计数的基因不太可能被评估为显著差异表达，因为低计数不能提供足够的统计证据来做出可靠的判断。 +因此这类基因可以在不丢失任何信息的情况下从分析中移除。 + +根据经验，我们要求一个基因至少在一些文库中有10个左右的数量，然后才会被认为在研究中表达出来。 +在实践中，筛选实际上是基于每百万次计数(count-per-million,CPM)的值，以避免偏向于在较大库中表达的基因，而不是在较小库中表达的基因。 +在 edgeR 包中，过滤可以使用一个考虑库大小和实验设计的函数来完成: +```{r keep} +keep <- filterByExpr(y, design) +table(keep) +``` +然后，“DGEList”对象被子集化，只保留未过滤的基因: +```{r filter} +y <- y[keep, , keep.lib.sizes=FALSE] +``` +使用`keep.lib.sizes=FALSE`函数，在过滤后重新计算库的大小。 +这通常是推荐的，尽管对下游分析的影响一般很小。 + +在分析中，“filterByExpr”函数保留感兴趣基因的最大数量，另外其他合理的过滤标准也是可行的。 +例如，“keep <- rowSums(y$counts)>50”是一个非常简单的标准，它可以保持总读取数大于50的基因。 +这将为该数据集提供与实际使用的过滤类似的下游结果。 + +另一个合理的过滤规则是计算每个基因的平均log-CPM，并选择截断值: +```{r aveLogCPM, fig.cap="Histogram of average log2 CPM."} +AveLogCPM <- aveLogCPM(y) +hist(AveLogCPM) +``` +同样，这可能会使该数据的下游结果与实际使用的过滤结果类似。 + +无论过滤规则是什么，它都应该独立于目标文件中的信息。 +它不应该参考哪个RNA文库属于哪个组，因为这样做会对后续的差异表达分析产生偏倚。 + + +成分偏倚标准化 {#normalization-for-composition-bias .unnumbered} +---------------------------------- + +使用“calcNormFactors”函数对M值的截尾均值(TMM)[@robinson2010tmm]进行标准化，该函数返回“DGEList”参数，其中只包含“norm.factors”。 +它计算一组标准化因子，每个样本一个，以消除库之间的成分偏差。 +这些因素和库大小的乘积定义了有效库大小，它在所有下游分析中替代了原始库大小。 + +```{r norm} +y <- calcNormFactors(y) +y$samples +``` + +所有库的归一化因子相乘得到统一。 +低于1的标准化因子表明，少数高计数基因垄断了测序，导致其他基因的计数低于通常给定库大小的情况。 +因此，该示例的有效库大小将按比例缩小。 +由此可见，L细胞哺乳期的样品的归一化因子较低。 +这是表明这些样本包含了一些高度上调的基因。 + +注释 +通常，我们发现TMM归一化对于几乎所有设计良好的mRNA基因表达实验都是令人满意的。 +单细胞转录组测序是个例外，其需要特殊的标准化途径[@lun2016pooling]。 +另一种不太常见的，需要特殊处理的研究类型是总体差异表达，超过一半的基因组在不同实验条件下以相同的方向差异表达[@wu2013mirna]。 +在设计良好的实验中，总体差异基因表达应该被避免。 +当无法避免时，则需要在实验中加入一些标准化参考，如spike-ins，以便进行可靠的标准化[@risso2014ruvseq]。 + + + +探索库之间的差异 {#exploring-differences-between-libraries .unnumbered} +--------------------------------------- + +RNA样本可以通过多维尺度(MDS)图在二维空间中聚类。 +这既是一个分析步骤，也是一个质量控制步骤，目的是探索不同样本表达谱之间的总体差异。 +这里我们描绘了MDS图来表示细胞群： + +```{r mdsplot, fig.width=7, fig.height=7, fig.cap="MDS plot showing distances between expression profiles. This provides a type of unsupervised clustering of the samples. In this case, the first dimension separates samples by cell type while the second dimension corresponds to mouse status."} +pch <- c(0,1,2,15,16,17) +colors <- rep(c("darkgreen", "red", "blue"), 2) +plotMDS(y, col=colors[group], pch=pch[group]) +legend("topleft", legend=levels(group), pch=pch, col=colors, ncol=2) +``` + +在MDS图中，每对样本之间的距离可以解释为最能区分这对样本的基因在样本之间的主要log-FC。 +默认情况下，主要log-FC定义为该对样本之间最大的500个log2-FC的均方根。 +从上图中可以看出，来自同一组的复制样本聚在一起，而来自不同组的样本则明显地分离。 +换句话说，组间的差异要比组内的差异大得多，这意味着组间可能存在统计上显著的差异。 +左侧基底细胞与右侧管腔细胞在x轴上的距离约为6个单位，对应两种细胞类型之间的z主要log-FC约为64倍。 +与基底细胞相比，未妊娠、妊娠期和哺乳期细胞表达谱的差异在腔细胞中似乎被放大了。 + +利用均值差(MD)图可以更接近地探索单个样本的表达谱。 +MD图可视化了库之间库大小调整的log-FC和库之间平均log表达。？？ +下面的命令生成一个MD图，将示例1与由所有其他示例的平均值构建的人工参考库进行比较: +```{r mdplot1, fig.width=7, fig.height=7, fig.cap="Mean-difference plot of log2-expression in sample 1 versus the average log2-expression across all other samples. Each point represents a gene, and the red line indicates a log-ratio of zero. The majority of points cluster around the red line."} +plotMD(y, column=1) +abline(h=0, col="red", lty=2, lwd=2) +``` +大部分基因以零对数倍变化（log-FC）线为中心。 +图左下角的对角线对应于第一个样本中计数为0、1、2等的基因。 + +将所有样品制作MD图作为质量检验是一种很好的做法。 +我们现在看到的腔哺乳期样品中，有一个标准化因子很低: +```{r mdplot11, fig.width=7, fig.height=7, fig.cap="Mean-difference plot of log2-expression in sample 11 versus the average log2-expression across all other samples. The plot shows a number of genes that are both highly expressed and highly up-regulated."} +plotMD(y, column=11) +abline(h=0, col="red", lty=2, lwd=2) +``` +对于这个样本，log-ratio显示出明显的正偏态，有许多高度上调的基因。 +特别地，在图的右上角有许多点，对应于与其他样本相比，该样本中表达量高且高度上调的基因。 +这些基因解释了为什么这个样本的标准化因子远低于1。 +相比之下，样本1的log-ratio略有负偏态，对应于大于1的标准化因子。 + + +分散估计 {#dispersion-estimation .unnumbered} +--------------------- + +edgeR 包使用负二项分布(NB)来模拟每个样本中每个基因的读计数。 +NB分布的分散参数解释了生物复制之间的差异 [@mccarthy2012edgerglm]. +edgeR 包估计每个基因的经验贝叶斯调节分散度。 +它还估计了一个共同的离散度，这是一个对所有基因的平均总体离散度估计，以及一个趋势离散度，其中一个基因的离散度是由其丰富度预测的。 +分散估计更容易从“estimateDisp”函数获得： + +```{r estimateDisp} +y <- estimateDisp(y, design, robust=TRUE) +``` +这将返回一个带有附加内容 (`common.dispersion`, `trended.dispersion` and `tagwise.dispersion`)的“DGEList”对象，以保持估计分散。 +在这里，“robust=TRUE”被用来保护经验贝叶斯估计不受异常大或个体离散度小的异常基因可能性的影响 [@phipson2016robust]. + +可以使用“plotBCV”函数来可视化分散估计： +```{r plotBCV, fig.width=7, fig.height=7, fig.cap="Scatterplot of the biological coefficient of variation (BCV) against the average abundance of each gene. The plot shows the square-root estimates of the common, trended and tagwise NB dispersions."} +plotBCV(y) +``` + +“plotBCV”图的纵轴表示平方根分散，也称为生物变异系数(BCV)[@mccarthy2012edgerglm]。 +在RNA-seq研究中，计数非常低的基因NB的分散度往往更高。 +对于数量较多的基因，其分散趋势随丰度的增加而平稳下降，并逐渐趋近于常数。 +从我们过去的经验来看，对于基因相同的小鼠或细胞系，BCV的渐近值往往在0.05到0.2之间，而对于人类受试者，其值略大(>0.3)。 + +NB模型可以用拟似然(QL)方法进行扩展，以考虑来自生物学和技术来源的基因特异性变异。 [@lund2012quasiseq; @lun2016delicious]. +在QL框架下，NB的分散趋势被用来描述所有基因的整体生物变异性，而在总体水平上下的基因特异性变异性被QL分散所捕捉。 +在QL方法中，不使用单个(tagwise) NB分散体。 + +使用“glmQLFit”函数来估计QL分散： +```{r glmQLFit} +fit <- glmQLFit(y, design, robust=TRUE) +head(fit$coefficients) +``` +这将返回一个“DGEGLM”对象，其中包含每个基因的GLM系数的估计值。 +它还包含了一些经验贝叶斯(EB)统计数据，包括QL离散趋势、压缩QL离散估计和先验自由度(df)。 +使用“plotQLDisp”函数来可视化QL分布： + +```{r QLDisp, fig.width=7, fig.height=7, fig.cap="A plot of the quarter-root QL dispersion against the average abundance of each gene. Estimates are shown for the raw (before EB moderation), trended and squeezed (after EB moderation) dispersions. Note that the QL dispersions and trend shown here are relative to the NB dispersion trend show in the previous figure produced using `plotBCV()`."} +plotQLDisp(fit) +``` + +QL函数以“limma”包调节方差的方式调节genewise QL离散度估计 [@smyth2004ebayes]. +原始QL离散度估计被压缩成一个总体趋势，这种适度降低了估计的不确定性，并提高了测试能力。 +压缩的程度由数据估计的先前df值决定。 +先前较大的df估计表明QL在基因间的分散性较小，这意味着应该进行较强的EB调节。 +先前较小的df估计值表明，真实的未知离散度是高度可变的，因此较弱的趋势缓和是适当的。 + +```{r df.prior} +summary(fit$df.prior) +``` + +通常建议在“glmQLFit”中设置“robust=TRUE”[@phipson2016robust] +这允许特定于基因的先验df估计，较低的离群基因值和较高的基因主体值。 +这减少了从具有极高或极低原始分散的基因中获得假阳性的机会，同时增加了检测基因主体差异表达的统计能力。 + +差异表达分析 {#differential-expression-analysis .unnumbered} +================================ + +差异表达测试 {#testing-for-differential-expression .unnumbered} +----------------------------------- + +下一步是在实验组中测试差异基因的表达。 +最有趣的比较之一是在基础怀孕组和哺乳期组之间。 +使用“makecontrast”函数可以方便地构造与任何指定比较对应的对比度: + +```{r B.PvsL} +B.LvsP <- makeContrasts(B.lactating-B.pregnant, levels=design) +``` + +在随后的结果中，正logFC表明与妊娠小鼠相比，哺乳期小鼠基因表达上调，而负logFC表明基因在妊娠小鼠中高表达。 + +我们将使用QL f检验而不是更常见的似然比检验(LRT)，因为它们通过考虑离散度估计中的不确定性来提供更严格的错误率控制: + +```{r glmQLFTest} +res <- glmQLFTest(fit, contrast=B.LvsP) +``` + +使用“topTags”函数来看top DE基因： + +```{r topTags} +topTags(res) +``` + +为了控制伪发现率(FDR)，采用benjamin - hochberg方法进行多次测试校正。 +Top DE基因“Csn1s2b”有一个较大的正logFC，这表明它在哺乳期的基底细胞中的表达远远高于怀孕小鼠。 +这种基因确实是牛奶中蛋白质的主要来源。 + +可以用testsdge来表示具有5% FDR的DE基因的总数。 +事实上，本比较中有超过5000个DE 基因： + +```{r decideTests} +is.de <- decideTestsDGE(res) +summary(is.de) +``` + +差分表达变化的大小可以用拟合的模型MD plot来表示: + +```{r plotMDfit, fig.width=7, fig.height=7, fig.cap="MD plot showing the log-fold change and average abundance of each gene. Significantly up and down DE genes are highlighted in red and blue, respectively."} +plotMD(res, status=is.de, values=c(1,-1), col=c("red","blue"), + legend="topright") +``` +每个基因的logFc与log2-CPM的平均丰度一起绘制，即“logCPM”在表的上方。 +高度差异表达的基因被重点显示。 + + +超过FC阈值的差异表达 {#differential-expression-above-a-fold-change-threshold .unnumbered} +----------------------------------------------------- + +无论差异有多小，“glmQLFTest”都会根据统计显著性来识别差异表达。 +出于某些目的，我们可能只对表达变化较大的基因感兴趣。 +上述分析发现，在基础妊娠组和哺乳期组之间有超过5000个DE基因。 +由于DE基因的数量如此之多，我们有必要把范围缩小到生物学上更有意义的基因。 + +一种常用的方法是同时应用FDR和logFC 来截断。 +然而，这往往有利于低表达基因，也不能正确控制FDR。 +更好和更严格的方法是修改统计测试，以便检测大于指定阈值的表达式变化。 +在“edgeR”包中，可以使用“glmTreat”函数来实现该功能。 +这个函数类似于微阵列的处理方法[@mccarthy2009treat]，但适用于NB框架。 +这里我们测试差分表达FC是否显著大于1.5，即logfc是否显著大于$\log_2(1.5)$: + +```{r treat} +tr <- glmTreat(fit, contrast=B.LvsP, lfc=log2(1.5)) +topTags(tr) +``` +注意，参数\code{lfc}是' ' log-fold-change ' '的缩写。 +以5%为FDR截断值，大约1100个基因被检测出，其FC显著大于1.5。 + + +```{r treatdecideTests} +is.de <- decideTestsDGE(tr) +summary(is.de) +``` +The p-values from `glmTreat` are larger than those from `glmQLFTest`, and the number of significantly DE genes is fewer, because it is testing an interval null hypothesis and requires stronger evidence for differential expression than does a conventional test. +It provides greater specificity for identifying the most important genes with large fold changes. + +The test results can be visualized in an MD plot: + +```{r plotMDtreat, fig.width=7, fig.height=7, fig.cap="MD plot showing the log-fold change and average abundance of each gene. Genes with fold-changes significantly greater than 1.5 are highlighted."} +plotMD(tr, status=is.de, values=c(1,-1), col=c("red","blue"), + legend="topright") +``` + +The \code{glmTreat} method evaluates variability as well as the magnitude of change of expression values and therefore is not equivalent to a simple fold change cutoff. +Nevertheless, all the statistically significant expression changes have logFC greater than 0.8 and almost all (97%) are greater than 0.9. +These values compare to the threshold value of $\log_2(1.5) = 0.58$. +In general, an estimated logFC must exceed the TREAT threshold by a number of standard errors for it to be called significant. +In other words, the whole confidence interval for the logFC must clear the threshold rather than just the estimated value itself. +It is better to interpret the threshold as *the FC below which we are definitely not interested in the gene* rather than *the FC above which we are interested in the gene*. + +The value of the FC threshold can be varied depending on the dataset. +In the presence of a huge number of DE genes, a relatively large FC threshold may be appropriate to narrow down the search to genes of interest. +In the absence of DE genes, on the other hand, a small or even no FC threshold shall be used. +If the threshold level is set to zero, then `glmTreat` becomes equivalent to `glmQLFTest` in the workflow shown here. + +In general, using `glmTreat` to reduce the number of DE genes is better than simply reducing the FDR cutoff, because `glmTreat` prioritizes genes with larger changes that are likely to be more biologically significant. +`glmTreat` can also be used with *edgeR* pipelines other than quasi-likelihood, although we don't demonstrate that here. + + +Heat map clustering {#heat-map-clustering .unnumbered} +------------------- + +Heatmaps are a popular way to display differential expression results for publication purposes. +To create a heatmap, we first convert the read counts into log2-counts-per-million (logCPM) values. +This can be done with the `cpm` function: + +```{r cpm} +logCPM <- cpm(y, prior.count=2, log=TRUE) +rownames(logCPM) <- y$genes$Symbol +colnames(logCPM) <- paste(y$samples$group, 1:2, sep="-") +``` +The introduction of `prior.count` is to avoid undefined values and to reduce the variability of the logCPM values for genes with low counts. +Larger values for `prior.count` shrink the logFCs for low count genes towards zero. + +We will create a heatmap to visualize the top 30 DE genes according to the TREAT test between `B.lactating` and `B.pregnant`. +The advantage of a heatmap is that it can display the expression pattern of the genes across all the samples. +Visualization of the results is aided by clustering together genes that have correlated expression patterns. +First we select the logCPM values for the 30 top genes: + +```{r order} +o <- order(tr$table$PValue) +logCPM <- logCPM[o[1:30],] +``` + +A heat map can then be produced by the `coolmap` function: +```{r heatmap, message=FALSE, fig.width=8, fig.height=12, fig.cap="Heat map across all the samples using the top 30 most DE genes between the basal lactating group and the basal pregnancy group."} +coolmap(logCPM, margins=c(7,7), lhei=c(1,6), lwid=c(1,3)) +``` +By default, `coolmap` clusters genes by correlation (highly correlated genes are closest) and clusters samples based on Euclidean distance between the expression values. +The coolmap has placed the basal samples on the right and the luminal samples on the left. +In each case, the virgin and pregnant samples are more like each other than the lactacting samples. +As expected, the two replicate samples from each group are clustered together. + +There a number of other R functions that can create heatmaps, but `coolmap` is the simplest because the default clustering metrics are chosen to be appropriate for expression data. + + +Analysis of deviance {#analysis-of-deviance .unnumbered} +-------------------- + +The differential expression analysis comparing two groups can be easily extended to comparisons between three or more groups. +This is done by creating a matrix of independent contrasts. +In this manner, users can perform a one-way analysis of deviance (ANODEV) for each gene [@mccullagh1989glms]. + +Suppose we want to compare the three groups in the luminal population, i.e., virgin, pregnant and lactating. +An appropriate contrast matrix can be created as shown below, to make pairwise comparisons between all three groups: + +```{r makeContrasts} +con <- makeContrasts( + L.PvsL = L.pregnant - L.lactating, + L.VvsL = L.virgin - L.lactating, + L.VvsP = L.virgin - L.pregnant, levels=design) +``` + +The QL F-test is then applied to identify genes that are DE between the three groups. +This combines the three pairwise comparisons into a single F-statistic and p-value. +The top set of significant genes can be displayed with `topTags`: + +```{r anovaQLFtest} +res <- glmQLFTest(fit, contrast=con) +topTags(res) +``` + +Note that the three contrasts of pairwise comparisons are linearly dependent. +Constructing the contrast matrix with any two of the contrasts would be sufficient for an ANODEV test. +If the contrast matrix contains all three possible pairwise comparisons, then only the log-fold changes of the first two contrasts are shown in the output of `topTags`. + +Complicated contrasts {#complicated-contrasts .unnumbered} +--------------------- + +The flexibility of the GLM framework makes it possible to specify arbitrary contrasts for differential expression tests. +Suppose we are interested in testing whether the change in expression between lactating and pregnant mice is the same for basal cells as it is for luminal cells. +In statistical terminology, this is the interaction effect between mouse status and cell type. +The contrast corresponding to this testing hypothesis can be made as follows. + +```{r complicatedContrasts} +con <- makeContrasts( + (L.lactating-L.pregnant)-(B.lactating-B.pregnant), + levels=design) +``` + +Then the QL F-test is conducted to identify genes that are DE under this contrast. +The top set of DE genes are viewed with `topTags`. + +```{r complicatedQLTest} +res <- glmQLFTest(fit, contrast=con) +topTags(res) +``` + +Pathway analysis {#pathway-analysis .unnumbered} +================ + +Gene ontology analysis {#gene-ontology-analysis .unnumbered} +---------------------- + + +We now consider the problem of interpreting the differential expression results in terms of higher order biological processes or molecular pathways. +One of the most common used resources is gene ontology (GO) databases, which annotate genes according to a dictionary of annotation terms. +A simple and often effective way to interpret the list of DE genes is to count the number of DE genes that are annotated with each possible GO term. +GO terms that occur frequently in the list of DE genes are said to be over-represented or enriched. + +In *edgeR*, GO analyses can be conveniently conducted using the `goana` function. +Here were apply `goana` to the output of the TREAT analysis comparing `B.lactating` to `B.pregant`. +The top most significantly enriched GO terms can be viewed with `topGO`. + +```{r goana} +go <- goana(tr, species="Mm") +topGO(go, n=15) +``` +The `goana` function automatically extracts DE genes from the `tr` object, and conducts overlap tests for the up- and down-regulated DE genes separately. +By default, an FDR cutoff of 5% is used when extracting DE genes, but this can be varied. +The row names of the output are the universal identifiers of the GO terms and the `Term` column gives the human-readable names of the terms. +The `Ont` column shows the ontology domain that each GO term belongs to. +The three domains are: biological process (BP), cellular component (CC) and molecular function (MF). +The `N` column represents the total number of genes annotated with each GO term. +The `Up` and `Down` columns indicate the number of genes within the GO term that are significantly up- and down-regulated in this differential expression comparison, respectively. +The `P.Up` and `P.Down` columns contain the p-values for over-representation of the GO term in the up- and down-regulated genes, respectively. +Note that the p-values are not adjusted for multiple testing---we would usually ignore GO terms with p-values greater than about $10^{-5}$. + +By default the output table from `topGO` is sorted by the minimum of `P.Up` and `P.Down`. +Other options are available. +For example, `topGO(go, sort="up")` lists the top GO terms that are over-represented in the up-regulated genes. +The domain of the enriched GO terms can also be specified by users. +For example, `topGO(go, ontology="BP")` restricts to the top GO terms belonging to the biological process domain +while `topGO(go, ontology="MF")` restricts to molecular function terms. + +The `goana` function uses the NCBI RefSeq annotation and requires the use of Entrez Gene Identifiers. + +KEGG pathway analysis {#kegg-pathway-analysis .unnumbered} +--------------------- + + +Another popular annotation database is the Kyoto Encyclopedia of Genes and Genomes (KEGG). +Much smaller than GO, this is a curated database of molecular pathways and disease signatures. +A KEGG analysis can be done exactly as for GO, but using the `kegga` function: + +```{r kegga} +keg <- kegga(tr, species="Mm") +topKEGG(keg, n=15, truncate=34) +``` + +The output from `topKEGG` is the same as from `topGO` except that row names become KEGG pathway IDs, `Term` becomes `Pathway` and there is no `Ont` column. +Both the GO and KEGG analyses show that the cell cycle pathway is strongly down-regulated upon lactation in mammary stem cells. + +By default, the `kegga` function automatically reads the latest KEGG annotation from the Internet each time it is run. +The KEGG database uses Entrez Gene Ids, and the `kegga` function assumes these are available as the row names of `tr`. + +FRY gene set tests {#fry-gene-set-tests .unnumbered} +------------------ + + +The GO and KEGG analyses shown above are relatively simple analyses that rely on a list of DE genes. +The list of DE genes is overlapped with the various GO and KEGG annotation terms. +The results will depend on the significance threshold that is used to assess differential expression. + +If the aim is to test for particular gene expression signatures or particular pathways, a more nuanced approach is to conduct a `roast` or `fry` gene set test [@wu2010roast]. +These functions test whether a set of genes is DE, assessing the whole set of genes as a whole. +Gene set tests consider all the genes in the specified set and do not depend on any pre-emptive significance cutoff. +The set of genes can be chosen to be representative of any pathway or phenotype of interest. + +`roast` gives p-values using random rotations of the residual space. +In the *edgeR* context, `fry` is generally recommended over `roast`. +`fry` gives an accurate analytic approximation to the results that `roast` would give, with default settings, if an extremely large number of rotations was used. + +Here, suppose we are interested in two GO terms related to cytokinesis. +Each GO term is used to define a set of genes annotated with that term. +The names of these terms are shown below: + +```{r select, message=FALSE} +library(GO.db) +cyt.go <- c("GO:0032465", "GO:0000281") +term <- select(GO.db, keys=cyt.go, columns="TERM") +term +``` + +The first step is to extract the genes associated with each GO term from the GO database. +This produces a list of two components, one for each GO term. +Each component is a vector of Entrez Gene IDs for that GO term: + +```{r GO2ALLEGS} +Rkeys(org.Mm.egGO2ALLEGS) <- cyt.go +cyt.go.genes <- as.list(org.Mm.egGO2ALLEGS) +``` + +Suppose the comparison of interest is between the virgin and lactating groups in the basal population. +We can use `fry` to test whether the cytokinesis GO terms are DE for this comparison: + +```{r fry} +B.VvsL <- makeContrasts(B.virgin-B.lactating, levels=design) +fry(y, index=cyt.go.genes, design=design, contrast=B.VvsL) +``` + +Each row of the output corresponds to a gene set. +The `NGenes` column provides the number of genes in each set. +The `Direction` column indicates the net direction of change. +The `PValue` column gives the two-sided p-value for testing whether the set is DE as a whole, either up or down. +The `PValue.Mixed` column gives a p-value for testing whether genes in the set tend to be DE, without regard to direction. +The `PValue` column is appropriate when genes in the set are expected to be co-regulated, all or most changing expression in the same direction. +The `PValue.Mixed` column is appropriate when genes in the set are not necessarily co-regulated or may be regulated in different directions for the contrast in question. +FDRs are calculated from the corresponding p-values across all sets. + +The results of a gene set test can be viewed in a barcode plot produced by the `barcodeplot` function. +Suppose visualization is performed for the gene set defined by the GO term GO:0032465: + +```{r barcode, fig.width=8, fig.height=5, fig.cap="Barcode plot showing enrichment of the GO term GO:0032465 in the basal virgin group compared to the basal lactating group. X-axis shows logFC for B.virgin vs B.lactating. Black bars represent genes annotated with the GO term. The worm shows relative enrichment."} +res <- glmQLFTest(fit, contrast=B.VvsL) +index <- rownames(fit) %in% cyt.go.genes[[1]] +barcodeplot(res$table$logFC, index=index, labels=c("B.lactating","B.virgin"), + main=cyt.go[1]) +``` + +In the plot, all genes are ranked from left to right by increasing log-fold change for the contrast and the genes within the gene set are represented by vertical bars, forming the barcode-like pattern. +The curve (or *worm*) above the barcode shows the relative local enrichment of the bars in each part of the plot. +The dotted horizontal line indicates neutral enrichment; the worm above the dotted line shows enrichment while the worm below the dotted line shows depletion. +In this particular barcode plot the worm shows enrichment on the right for positive logFCs, and depletion on the left for negative logFCs. +The conclusion is that genes associated with this GO term tend to be up-regulated in the basal cells of virgin mice compared to lactating mice, confirming the result of the `fry` test above. + +Camera gene set enrichment analysis {#camera-gene-set-enrichment-analysis .unnumbered} +----------------------------------- + +Finally we demonstrate a gene set enrichment style analysis using the Molecular Signatures Database (MSigDB) [@subramanian2005gsea]. +We will use the C2 collection of the MSigDB, which is a collection of nearly 5000 curated gene sets, each representing the molecular signature of a particular biological process or phenotype. +The MSigDB itself is purely human, but the Walter and Eliza Hall Institute (WEHI) maintains a mouse version of the database. +We load the mouse version of the C2 collection from the WEHI website: +```{r camera} +load(url("http://bioinf.wehi.edu.au/software/MSigDB/mouse_c2_v5p1.rdata")) +``` +This will load `Mm.c2`, which is a list of gene sets, each a vector of Entrez Ids. +This can be converted to a list of index numbers: + +```{r id2indices.Mm.c2} +idx <- ids2indices(Mm.c2,id=rownames(y)) +``` + +First we compare virgin stem cells to virgin luminal cells: + +```{r cam.BVvsLV} +BvsL.v <- makeContrasts(B.virgin - L.virgin, levels=design) +cam <- camera(y, idx, design, contrast=BvsL.v, inter.gene.cor=0.01) +options(digits=2) +head(cam,14) +``` +With a large gene set collection, setting `inter.gene.cor = 0.01` gives a good compromise between biological interpretability and FDR control. +As expected, the mammary stem cell and mammary luminal cell signatures from Lim et al [@lim2010transcriptome] are top-ranked, and in the expected directions. + +We can visualize the top signature, combining the up and down mammary stem cell signatures to make a bi-directional signature set: + +```{r barcode2, fig.width=8, fig.height=6.4, fig.cap="Barcode plot showing strong enrichment of mammary stem cell signature in the stem cell vs luminal cell comparison. Red bars show up signature genes, blue bars show down genes. The worms show relative enrichment."} +res <- glmQLFTest(fit, contrast=BvsL.v) +barcodeplot(res$table$logFC, + index=idx[["LIM_MAMMARY_STEM_CELL_UP"]], + index2=idx[["LIM_MAMMARY_STEM_CELL_DN"]], + labels=c("L.virgin","B.virgin"), + main="LIM_MAMMARY_STEM_CELL", + alpha=1) +``` + + +Packages used {#packages-used .unnumbered} +============= + + +This workflow depends on various packages from version 3.3 of the Bioconductor project, running on R version 3.3.0 or higher. +The complete list of the packages used for this workflow are shown below: + +```{r info} +sessionInfo() +``` + + +Read alignment and quantification {#read-alignment-and-quantification .unnumbered} +================================= + +Download raw sequence files from the SRA {#download-raw-sequence-files-from-the-sra .unnumbered} +---------------------------------------- + +We now revisit the question of recreating the matrix of read counts from the raw sequence reads. +Unlike the above workflow, which works for any version of R, read alignment requires Unix or Mac OS and, in practice, a high performance Unix server is recommended. +Read alignment and read counting require only one Bioconductor package, *Rsubread*. +However the `fastq-dump` utility from the SRA Toolkit is also required to convert from SRA to FASTQ format. +This can be downloaded from the NCBI website () and installed on any Unix system. + +The first task is to download the raw sequence files, which are stored in SRA format on the SRA repository. +The SRA files need to be unpacked into FASTQ format using the `fastq-dump` utility. +The following R code makes a system call to `fastq-dump` to download each SRA file and convert it to FASTQ format: + + for (sra in targets$SRA) system(paste("fastq-dump", sra)) + +The `fastq-dump` utility automatically downloads the specified SRA data set from the internet. +The above code will produce 12 FASTQ files, in the current working directory, with file names given by the following vector: + + all.fastq <- paste0(targets$SRA, ".fastq") + + +Accuracy of base-calling {#accuracy-of-base-calling .unnumbered} +------------------------ + + +Sequencers typically store base-calling quality scores for each read in the FASTQ files. +*Rsubread*'s `qualityScores` function can be used to extract these scores from any particular file: + + QS <- qualityScores("SRR1552444.fastq") + +The boxplot function provides a compact way to view the quality scores by position across all reads: + + boxplot(QS, ylab="Quality score", xlab="Base position", + main="SRR1552444.fastq", cex=0.25, col="orange") + + +
+Boxplots of quality scores by base position for the first FASTQ file. +

+Figure 12: Boxplots of quality scores by base position for the first FASTQ file.
+

+
+ + +The vertical axis shows the Phred quality score, equal to $-10\log_{10}(p)$ where $p$ is the probability of an erroneous call. +The maximum possible value is 40, and all values above 10 correspond to extremely small error probabilities. +The horizontal axis shows position within a read. +The file contains 100bp single-end reads, so the scale is from 1 to 100. +The plot displays a compact boxplot at each base position. +As is very commonly observed, the quality scores are best in the middle of the reads and decrease slightly towards the start and end of the reads. +However the quality remains generally good even near the ends of the reads: the scores would need to be very much lower than this before they would cause problems for the alignment. +Similar plots can be made for each of the FASTQ files. + + +Build a genome index {#build-a-genome-index .unnumbered} +-------------------- + +Before the sequence reads can be aligned, we need to build an index for the GRCm38/mm10 (Dec 2011) build of the mouse genome. +Most laboratories that use *Rsubread* regularly will already have an index file prepared, as this is a once-off operation for each genome release. +If you are using *Rsubread* for mouse for the first time, then the latest mouse genome build can be downloaded from the NCBI webpage +(). +The full URL is: + + ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz + +(Note that this link is for patch 6 of the RefSeq mm10 build, which is valid at the time of writing in May 2018. +The link will change as new patches are released periodically.) +An index can then be built by: + + library(Rsubread) + buildindex(basename = "mm10", + reference = "GCF_000001635.26_GRCm38.p6_genomic.fna.gz") + +Aligning reads {#aligning-reads .unnumbered} +-------------- + +The sequence reads can now be aligned to the mouse genome using the `align` function: + + all.bam <- sub(".fastq", ".bam", all.fastq) + align(index="mm10", readfile1=all.fastq, input_format="FASTQ", + output_file=all.bam) + +This produces a set of BAM files containing the read alignments for each RNA library. +The mapping proportions can be summarized by the `propmapped` function: + + > propmapped(all.bam) + Samples NumTotal NumMapped PropMapped + 1 SRR1552450.bam 30109290 26577308 0.883 + 2 SRR1552451.bam 28322351 24794251 0.875 + 3 SRR1552452.bam 31688348 27937620 0.882 + 4 SRR1552453.bam 29614284 26074034 0.880 + 5 SRR1552454.bam 27225012 24381742 0.896 + 6 SRR1552455.bam 25433157 22813815 0.897 + 7 SRR1552444.bam 27919481 23927833 0.857 + 8 SRR1552445.bam 29731031 25487822 0.857 + 9 SRR1552446.bam 29879070 25500318 0.853 + 10 SRR1552447.bam 29245388 25187577 0.861 + 11 SRR1552448.bam 31425424 27326500 0.870 + 12 SRR1552449.bam 31276061 27204156 0.870 + +Ideally, the proportion of mapped reads should be above 80%. +By default, only reads with unique mapping locations are reported by *Rsubread* as being successfully mapped. +Restricting to uniquely mapped reads is recommended, as it avoids spurious signal from non-uniquely mapped reads derived from, e.g., repeat regions. + + +Quantifying read counts for each gene {#quantifying-read-counts-for-each-gene .unnumbered} +------------------------------------- + +The read counts for each gene can be quantified using the `featureCounts` function in *Rsubread*. +Conveniently, the *Rsubread* package includes inbuilt NCBI RefSeq annotation of the mouse and human genomes. +`featureCounts` generates a matrix of read counts for each gene in each sample: + + fc <- featureCounts(all.bam, annot.inbuilt="mm10") + +The output is a simple list, containing the matrix of counts (`counts`), a data frame of gene characteristics (`annotation`), a vector of file names (`targets`) and summary mapping statistics (`stat`): + + > names(fc) + [1] "counts" "annotation" "targets" "stat" + +The row names of `fc$counts` are the Entrez gene identifiers for each gene. +The column names are the output file names from `align`, which we simplify here for brevity: + + > colnames(fc$counts) <- rownames(targets) + +The first six rows of the counts matrix are shown below. + + > head(fc$counts) + MCL1.DG MCL1.DH MCL1.DI MCL1.DJ MCL1.DK MCL1.DL MCL1.LA + 497097 438 299 65 237 354 287 0 + 100503874 1 0 1 1 0 4 0 + 100038431 0 0 0 0 0 0 0 + 19888 1 1 0 0 0 0 10 + 20671 106 181 82 104 43 83 16 + 27395 309 232 339 290 291 270 558 + MCL1.LB MCL1.LC MCL1.LD MCL1.LE MCL1.LF + 497097 0 0 0 0 0 + 100503874 0 0 0 0 0 + 100038431 0 0 0 0 0 + 19888 3 10 2 0 0 + 20671 25 18 8 3 10 + 27395 468 488 332 312 344 + +Finally, a `DGEList` object can be assembled by: + + y <- DGEList(fc$counts, group=group) + y$genes <- fc$annotation[, "Length", drop=FALSE] + + +Data and software availability {#data-and-software-availability .unnumbered} +============================== + +This html document was generated as a vignette of the RnaSeqGeneEdgeRQL package, which can be installed by running + + install.packages("BiocManager") + BiocManager::install("RnaSeqGeneEdgeRQL") + +at the R prompt. +Except for the targets file `targets.txt`, all data analyzed in the workflow is read automatically from public websites as part of the code. +All software used is publicly available as part of Bioconductor (Release 3.3 or later), except for the `fastq-dump` utility, which can be downloaded from NCBI website as described in the text. +The article includes the complete code necessary to reproduce the analyses shown. + +The original Bioconductor 3.3 version of this workflow, including a pdf article version and a corresponding Rnw file, is available from and [@chen2016quasiworkflow]. + + +Author contributions {#author-contributions .unnumbered} +==================== + +All authors developed and tested the code workflow. +All authors wrote the article. + +Grant information {#grant-information .unnumbered} +================= + +This work was supported by the National Health and Medical Research Council +(Fellowship 1058892 and Program 1054618 to G.K.S, Independent Research Institutes Infrastructure Support to the Walter and Eliza Hall Institute) +and by a Victorian State Government Operational Infrastructure Support Grant. + +Acknowledgments {#acknowledgments .unnumbered} +=============== + +The authors thank Wei Shi and Yang Liao for advice with *Rsubread* and Yifang Hu for creating the mouse version of the MSigDB. +We also wish to acknowledge the early developers of the quasi-likelihood testing approach. +Davis McCarthy authored the first versions of the quasi-likelihood functions in *edgeR* in February 2011 and Steve Lund and Dan Nettleton worked independently on *QuasiSeq* around the same time. +We thank Steve Lund and Dan Nettleton for a valuable and enjoyable collaboration that lead to Lund et al [@lund2012quasiseq]. + +How to cite {#how-to-cite .unnumbered} +==================== + +Please cite this workflow article as + +Chen, Yunshun, A. T. L. Lun, and G. K. Smyth. 2016. "From Reads to Genes to Pathways: Differential Expression Analysis of RNA-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline." *F1000Research* 5:1438. + +References {#references .unnumbered} +========== + +lter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia - name: Aaron T. L. Lun affiliation: Cancer Research UK Cambridge Institute, University of Cambridge, Li Ka Shing Centre, Robinson Way, CB20RE, Cambridge, United Kingdom - name: Gordon K. Smyth From 6d651a2eb60545410ae4e2c88f7a3071d695cf9f Mon Sep 17 00:00:00 2001 From: restore1997 <47097886+restore1997@users.noreply.github.com> Date: Thu, 1 Aug 2019 23:01:22 +0800 Subject: [PATCH 2/3] Update 10RnaSeqGeneEdgeRQL_en.Rmd --- BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd | 1347 +++--------------- 1 file changed, 204 insertions(+), 1143 deletions(-) diff --git a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd index fdbcdba..f2ae507 100644 --- a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd +++ b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd @@ -1,8 +1,4 @@ --- -title: "From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline" -author: - - name: Yunshun Chen - affiliation: The Wa--- 题目：「从序列到基因到通路：使用Rsubread包和edgeR包拟似然管道进行RNA-Seq差异表达分析」 作者： - 姓名： Yunshun Chen @@ -447,78 +443,78 @@ topTags(tr) is.de <- decideTestsDGE(tr) summary(is.de) ``` -The p-values from `glmTreat` are larger than those from `glmQLFTest`, and the number of significantly DE genes is fewer, because it is testing an interval null hypothesis and requires stronger evidence for differential expression than does a conventional test. -It provides greater specificity for identifying the most important genes with large fold changes. +来自“glmTreat”的p值大于来自“glmQLFTest”的p值，显著DE基因的数量也更少，因为它测试的是区间零假设，需要比传统测试更强的证据来证明差异表达。 +它为识别具有较大FC的最重要基因提供了更大的特异性。 -The test results can be visualized in an MD plot: +测试结果可由一个MD图观察到： ```{r plotMDtreat, fig.width=7, fig.height=7, fig.cap="MD plot showing the log-fold change and average abundance of each gene. Genes with fold-changes significantly greater than 1.5 are highlighted."} plotMD(tr, status=is.de, values=c(1,-1), col=c("red","blue"), legend="topright") ``` -The \code{glmTreat} method evaluates variability as well as the magnitude of change of expression values and therefore is not equivalent to a simple fold change cutoff. -Nevertheless, all the statistically significant expression changes have logFC greater than 0.8 and almost all (97%) are greater than 0.9. -These values compare to the threshold value of $\log_2(1.5) = 0.58$. -In general, an estimated logFC must exceed the TREAT threshold by a number of standard errors for it to be called significant. -In other words, the whole confidence interval for the logFC must clear the threshold rather than just the estimated value itself. -It is better to interpret the threshold as *the FC below which we are definitely not interested in the gene* rather than *the FC above which we are interested in the gene*. +\code{glmTreat}方法评估变异性以及表达值改变的大小，因此并不等同于简单的FC截断。 +然而，所有有统计学意义的表达变化logFC大于0.8，且几乎所有(97%)大于0.9。 +这些值与$\log_2(1.5) = 0.58$的阈值进行比较。 +一般来说，估计的logFC必须超过许多标准错误的治疗阈值，才能称为显著。 +换句话说，logFC的整个置信区间必须清除阈值，而不仅仅是估计值本身。 +It is better to interpret the threshold as *the FC below which we are definitely not interested in the gene* rather than *the FC above which we are interested in the gene*.最好将阈值解释为“我们对该基因绝对不感兴趣的FC”，而不是*我们对该基因感兴趣的FC *。 -The value of the FC threshold can be varied depending on the dataset. -In the presence of a huge number of DE genes, a relatively large FC threshold may be appropriate to narrow down the search to genes of interest. -In the absence of DE genes, on the other hand, a small or even no FC threshold shall be used. -If the threshold level is set to zero, then `glmTreat` becomes equivalent to `glmQLFTest` in the workflow shown here. +FC阈值可以根据数据集的不同而变化。 +在存在大量DE基因的情况下，一个相对较大的FC阈值可能适合将搜索范围缩小到感兴趣的基因。 +另一方面，在缺乏DE基因的情况下，应使用小的甚至没有FC阈值。 +如果阈值水平设置为零，那么在这里显示的工作流中，“glmTreat”就等同于“glmQLFTest”。 -In general, using `glmTreat` to reduce the number of DE genes is better than simply reducing the FDR cutoff, because `glmTreat` prioritizes genes with larger changes that are likely to be more biologically significant. -`glmTreat` can also be used with *edgeR* pipelines other than quasi-likelihood, although we don't demonstrate that here. +总的来说，使用“glmTreat”来减少DE基因的数量比简单地减少FDR截断要好，因为“glmTreat”优先考虑那些可能具有更大生物学意义的变化的基因。 +“glmTreat”也可以用于除准似然外的*edgeR*管道，不过我们在这里没有演示。 -Heat map clustering {#heat-map-clustering .unnumbered} +热图聚类 {#heat-map-clustering .unnumbered} ------------------- -Heatmaps are a popular way to display differential expression results for publication purposes. -To create a heatmap, we first convert the read counts into log2-counts-per-million (logCPM) values. -This can be done with the `cpm` function: +出于发表的目的，人们经常使用热图来表示差异表达基因的结果。 +要创建热图，首先将读取计数转换为log2-count -per-million (logCPM)值。 +这可以通过“cpm”函数实现: ```{r cpm} logCPM <- cpm(y, prior.count=2, log=TRUE) rownames(logCPM) <- y$genes$Symbol colnames(logCPM) <- paste(y$samples$group, 1:2, sep="-") ``` -The introduction of `prior.count` is to avoid undefined values and to reduce the variability of the logCPM values for genes with low counts. -Larger values for `prior.count` shrink the logFCs for low count genes towards zero. +“prior.count”的引入是为了避免未定义值以及减少低计数基因的logCPM值的变异性。 +较大的“prior.count”值将低计数基因的logFCs缩小到0。 -We will create a heatmap to visualize the top 30 DE genes according to the TREAT test between `B.lactating` and `B.pregnant`. -The advantage of a heatmap is that it can display the expression pattern of the genes across all the samples. -Visualization of the results is aided by clustering together genes that have correlated expression patterns. -First we select the logCPM values for the 30 top genes: +我们将根据“B.lactating”和“B.lactating”间的治疗测试创建一个热图，以可视化top 30 DE 。 +heatmap的优势在于它可以显示所有样本中基因的表达谱。 +将具有相关表达谱的基因聚在一起，有助于结果的可视化。 +首先我们选择30个 top 基因的logCPM值： ```{r order} o <- order(tr$table$PValue) logCPM <- logCPM[o[1:30],] ``` -A heat map can then be produced by the `coolmap` function: +然后使用“coolmap”函数来生成热图： ```{r heatmap, message=FALSE, fig.width=8, fig.height=12, fig.cap="Heat map across all the samples using the top 30 most DE genes between the basal lactating group and the basal pregnancy group."} coolmap(logCPM, margins=c(7,7), lhei=c(1,6), lwid=c(1,3)) ``` -By default, `coolmap` clusters genes by correlation (highly correlated genes are closest) and clusters samples based on Euclidean distance between the expression values. -The coolmap has placed the basal samples on the right and the luminal samples on the left. -In each case, the virgin and pregnant samples are more like each other than the lactacting samples. -As expected, the two replicate samples from each group are clustered together. +默认情况下，“coolmap”通过相关性(高度相关的基因是最接近的)聚类基因，并基于表达值之间的欧氏距离聚类样本。 +Coolmap将基底细胞样本置于右侧，腔细胞样本置于左侧。 +在每种情况下，未妊娠和怀孕的样本都比泌乳的样本更相似。 +正如预期的那样，来自每个组的两个复制样本聚集在一起。 -There a number of other R functions that can create heatmaps, but `coolmap` is the simplest because the default clustering metrics are chosen to be appropriate for expression data. +还有许多其他的R函数可以创建heatmap，但是“coolmap”是最简单的，因为其默认的聚类指标被选择为适合于表达式数据。 -Analysis of deviance {#analysis-of-deviance .unnumbered} +异常分析 {#analysis-of-deviance .unnumbered} -------------------- -The differential expression analysis comparing two groups can be easily extended to comparisons between three or more groups. -This is done by creating a matrix of independent contrasts. -In this manner, users can perform a one-way analysis of deviance (ANODEV) for each gene [@mccullagh1989glms]. +比较两组之间的差异表达分析可以很容易地扩展到三个或更多组之间的比较。 +这是通过创建一个独立的对比矩阵来实现的。 +通过这种方式，用户可以对每个基因执行单向异常分析(ANODEV) [@mccullagh1989glms]。 -Suppose we want to compare the three groups in the luminal population, i.e., virgin, pregnant and lactating. -An appropriate contrast matrix can be created as shown below, to make pairwise comparisons between all three groups: +假设我们要比较腔细胞种群中的三组，即未妊娠，怀孕和哺乳期。 +如下图所示，可以创建一个合适的对比矩阵，对三组进行两两比较: ```{r makeContrasts} con <- makeContrasts( @@ -527,26 +523,26 @@ con <- makeContrasts( L.VvsP = L.virgin - L.pregnant, levels=design) ``` -The QL F-test is then applied to identify genes that are DE between the three groups. -This combines the three pairwise comparisons into a single F-statistic and p-value. -The top set of significant genes can be displayed with `topTags`: +然后使用QL f检验来识别三组之间的DE基因。 +这将三个成对的比较组合成一个单f统计量和p值。 +最显著的一组基因可以用“topTags”显示: ```{r anovaQLFtest} res <- glmQLFTest(fit, contrast=con) topTags(res) ``` -Note that the three contrasts of pairwise comparisons are linearly dependent. -Constructing the contrast matrix with any two of the contrasts would be sufficient for an ANODEV test. -If the contrast matrix contains all three possible pairwise comparisons, then only the log-fold changes of the first two contrasts are shown in the output of `topTags`. +注意，成对比较的三个对比是线性相关的。 +用任意两个对比构建对比矩阵就足以进行ANODEV测试。 +如果对比矩阵包含所有三种可能的成对比较，那么在“topTags”的输出中只显示前两个对比的log-FC。 -Complicated contrasts {#complicated-contrasts .unnumbered} +复杂对比 {#complicated-contrasts .unnumbered} --------------------- -The flexibility of the GLM framework makes it possible to specify arbitrary contrasts for differential expression tests. -Suppose we are interested in testing whether the change in expression between lactating and pregnant mice is the same for basal cells as it is for luminal cells. -In statistical terminology, this is the interaction effect between mouse status and cell type. -The contrast corresponding to this testing hypothesis can be made as follows. +GLM框架的灵活性使得为差异表达式测试指定任意对比成为可能。 +假设我们有兴趣测试哺乳和怀孕小鼠的基底细胞和腔细胞的表达变化是否相同。 +在统计学术语中，这是小鼠状态和细胞类型之间的交互作用。 +与此检验假设相对应的对比如下: ```{r complicatedContrasts} con <- makeContrasts( @@ -554,92 +550,92 @@ con <- makeContrasts( levels=design) ``` -Then the QL F-test is conducted to identify genes that are DE under this contrast. -The top set of DE genes are viewed with `topTags`. +然后进行QL f检验，以确定在这种对比下DE的基因。 +Top DE基因组用“topTags”来查看。 ```{r complicatedQLTest} res <- glmQLFTest(fit, contrast=con) topTags(res) ``` -Pathway analysis {#pathway-analysis .unnumbered} +通路分析 {#pathway-analysis .unnumbered} ================ -Gene ontology analysis {#gene-ontology-analysis .unnumbered} +基因本体论分析 {#gene-ontology-analysis .unnumbered} ---------------------- -We now consider the problem of interpreting the differential expression results in terms of higher order biological processes or molecular pathways. -One of the most common used resources is gene ontology (GO) databases, which annotate genes according to a dictionary of annotation terms. -A simple and often effective way to interpret the list of DE genes is to count the number of DE genes that are annotated with each possible GO term. -GO terms that occur frequently in the list of DE genes are said to be over-represented or enriched. +我们现在考虑的问题是如何用更高阶的生物学过程或分子途径来解释差异表达的结果。 +基因本体(GO)数据库根据注释术语的字典对基因进行注释,是最常用的资源之一。 +解释DE基因列表的一个简单且通常有效的方法是计算每个可能的GO术语注释的DE基因的数量。 +常出现在DE基因列表中的GO术语被认为是过度表达或富集的。 -In *edgeR*, GO analyses can be conveniently conducted using the `goana` function. -Here were apply `goana` to the output of the TREAT analysis comparing `B.lactating` to `B.pregant`. -The top most significantly enriched GO terms can be viewed with `topGO`. +在“edgeR”包中，可以使用“goana”函数方便地进行GO分析。 +这里我们应用“goana”函数到“B.lactating”和“B.lactating”比较治疗分析的输出结果中。 +可以通过“topGO”查看最显著富集的GO术语。 ```{r goana} go <- goana(tr, species="Mm") topGO(go, n=15) ``` -The `goana` function automatically extracts DE genes from the `tr` object, and conducts overlap tests for the up- and down-regulated DE genes separately. -By default, an FDR cutoff of 5% is used when extracting DE genes, but this can be varied. -The row names of the output are the universal identifiers of the GO terms and the `Term` column gives the human-readable names of the terms. -The `Ont` column shows the ontology domain that each GO term belongs to. -The three domains are: biological process (BP), cellular component (CC) and molecular function (MF). -The `N` column represents the total number of genes annotated with each GO term. -The `Up` and `Down` columns indicate the number of genes within the GO term that are significantly up- and down-regulated in this differential expression comparison, respectively. -The `P.Up` and `P.Down` columns contain the p-values for over-representation of the GO term in the up- and down-regulated genes, respectively. -Note that the p-values are not adjusted for multiple testing---we would usually ignore GO terms with p-values greater than about $10^{-5}$. - -By default the output table from `topGO` is sorted by the minimum of `P.Up` and `P.Down`. -Other options are available. -For example, `topGO(go, sort="up")` lists the top GO terms that are over-represented in the up-regulated genes. -The domain of the enriched GO terms can also be specified by users. -For example, `topGO(go, ontology="BP")` restricts to the top GO terms belonging to the biological process domain -while `topGO(go, ontology="MF")` restricts to molecular function terms. - -The `goana` function uses the NCBI RefSeq annotation and requires the use of Entrez Gene Identifiers. - -KEGG pathway analysis {#kegg-pathway-analysis .unnumbered} +goana函数自动从tr对象中提取DE基因，并分别对上调和下调的DE基因进行重叠测试。 +默认情况下，提取DE基因时使用5% FDR作为截断值，该值也可以变化。 +输出的行名称是GO术语的通用标识符，而“Term”列给出了人类可读的术语名称。 +“Ont”列显示了每个GO术语所属的本体域。 +这三个领域分别是:生物过程(BP)、细胞成分(CC)和分子功能(MF)。 +“N”列表示每个GO项注释的基因总数。 +“上”和“下”列分别表示在这个差异表达比较中，GO项中显著上调和下调的基因数量。 +“P.Up”和“P.Down”列分别包含了在上调基因和下调基因中过表达GO术语的p值。 +注意，p值没有针对多重测试进行调整——我们通常会忽略p值大于$10^(-5)$的GO术语。 + +默认情况下，“topGO”的输出表按照“P.Up”和“P.Down”的最小值排列。 +还有其他选项。 +比如，“topGO(go, sort="up")”列举了过度富集上调基因的top GO术语。 +丰富的GO术语的域也可以由用户指定。 +比如，“topGO(go, ontology="BP")”限制top GO术语在生物过程领域内。 +而“topGO(go, ontology="MF")”限制在分子功能领域内。 + +“goana”函数使用NCBI RefSeq注释，并要求使用Entrez基因标识符。 + +KEGG 通路分析 {#kegg-pathway-analysis .unnumbered} --------------------- -Another popular annotation database is the Kyoto Encyclopedia of Genes and Genomes (KEGG). -Much smaller than GO, this is a curated database of molecular pathways and disease signatures. -A KEGG analysis can be done exactly as for GO, but using the `kegga` function: +另一个流行的注释数据库是《京都基因与基因组百科全书》(KEGG)。 +这是一个分子途径和疾病特征的数据库,比GO小很多。 +KEGG分析可以像GO一样完成，但是使用' kegga '函数: ```{r kegga} keg <- kegga(tr, species="Mm") topKEGG(keg, n=15, truncate=34) ``` -The output from `topKEGG` is the same as from `topGO` except that row names become KEGG pathway IDs, `Term` becomes `Pathway` and there is no `Ont` column. -Both the GO and KEGG analyses show that the cell cycle pathway is strongly down-regulated upon lactation in mammary stem cells. +“topKEGG”的输出与“topGO”相同，只是行名变成了KEGG通路ID，“Term”变成了“pathway”，没有“Ont”列。 +GO和KEGG分析均表明，乳腺干细胞泌乳后细胞周期通路出现明显下调。 -By default, the `kegga` function automatically reads the latest KEGG annotation from the Internet each time it is run. -The KEGG database uses Entrez Gene Ids, and the `kegga` function assumes these are available as the row names of `tr`. +默认情况下，“kegga”函数每次运行时都会自动从互联网上读取最新的KEGG注释。 +KEGG数据库使用Entrez基因ID，而' kegga '函数假定这些ID可用作' tr '的行名。 -FRY gene set tests {#fry-gene-set-tests .unnumbered} +FRY 基因集测试 {#fry-gene-set-tests .unnumbered} ------------------ -The GO and KEGG analyses shown above are relatively simple analyses that rely on a list of DE genes. -The list of DE genes is overlapped with the various GO and KEGG annotation terms. -The results will depend on the significance threshold that is used to assess differential expression. +上面显示的GO和KEGG分析是相对简单的分析，依赖于DE基因列表。 +DE基因列表与各种GO和KEGG注释术语重叠。 +结果将取决于用于评估差异表达的显著性阈值。 -If the aim is to test for particular gene expression signatures or particular pathways, a more nuanced approach is to conduct a `roast` or `fry` gene set test [@wu2010roast]. -These functions test whether a set of genes is DE, assessing the whole set of genes as a whole. -Gene set tests consider all the genes in the specified set and do not depend on any pre-emptive significance cutoff. -The set of genes can be chosen to be representative of any pathway or phenotype of interest. +如果目标是测试特定的基因表达特征或特定的途径，一个更微妙的方法是进行“roast”或“fry”基因集测试 [@wu2010roast]。 +这些功能测试一组基因是否DE，评估整个基因集作为一个整体。 +基因集测试考虑指定集中的所有基因，不依赖于任何先前制定的显著性截断。 +这组基因可以被选择为代表任何感兴趣的途径或表型。 -`roast` gives p-values using random rotations of the residual space. -In the *edgeR* context, `fry` is generally recommended over `roast`. -`fry` gives an accurate analytic approximation to the results that `roast` would give, with default settings, if an extremely large number of rotations was used. +“roast”使用剩余空间的随机旋转给出p值。 +在“edgeR”包的文章中，“fry”通常比“roast”更受欢迎。 +“fry”给出了一个精确的近似分析结果，在默认设置下，如果使用了大量的旋转，“roast”将给出一个精确的近似分析结果。 -Here, suppose we are interested in two GO terms related to cytokinesis. -Each GO term is used to define a set of genes annotated with that term. -The names of these terms are shown below: +这里，假设我们对与细胞分裂有关的两个GO术语感兴趣。 +每个GO术语都用来定义一组用该术语标注的基因。 +这些术语的名称如下: ```{r select, message=FALSE} library(GO.db) @@ -648,34 +644,34 @@ term <- select(GO.db, keys=cyt.go, columns="TERM") term ``` -The first step is to extract the genes associated with each GO term from the GO database. -This produces a list of two components, one for each GO term. -Each component is a vector of Entrez Gene IDs for that GO term: +第一步是从GO数据库中提取与每个GO术语相关的基因。 +这将生成一个包含两个成分的列表，每个成分对应一个GO项。 +每个成分都是该GO术语的Entrez基因ID的向量: ```{r GO2ALLEGS} Rkeys(org.Mm.egGO2ALLEGS) <- cyt.go cyt.go.genes <- as.list(org.Mm.egGO2ALLEGS) ``` -Suppose the comparison of interest is between the virgin and lactating groups in the basal population. -We can use `fry` to test whether the cytokinesis GO terms are DE for this comparison: +假设兴趣的比较是在基细胞群中的未妊娠组和哺乳期组之间。 +我们可以用“fry”来测试细胞分裂的GO术语在这个比较中是否为DE: ```{r fry} B.VvsL <- makeContrasts(B.virgin-B.lactating, levels=design) fry(y, index=cyt.go.genes, design=design, contrast=B.VvsL) ``` -Each row of the output corresponds to a gene set. -The `NGenes` column provides the number of genes in each set. -The `Direction` column indicates the net direction of change. -The `PValue` column gives the two-sided p-value for testing whether the set is DE as a whole, either up or down. -The `PValue.Mixed` column gives a p-value for testing whether genes in the set tend to be DE, without regard to direction. -The `PValue` column is appropriate when genes in the set are expected to be co-regulated, all or most changing expression in the same direction. -The `PValue.Mixed` column is appropriate when genes in the set are not necessarily co-regulated or may be regulated in different directions for the contrast in question. -FDRs are calculated from the corresponding p-values across all sets. +每一行输出对应一个基因集。 +“NGenes”列提供了每组基因的数量。 +“Direction”列表示变化的网格？方向。 +“PValue”列给出了双面p值，用于测试集合是作为一个整体DE，上调还是下调。 +“PValue.Mixed”列给出了一个p值，用于测试集合中的基因是否倾向于DE，而不考虑方向。 +当集合中的基因被期望是共同调控的，所有或大部分的表达都在同一个方向上改变时,“PValue”列是合适的，。 +当集合中的基因不一定是共同调控的，或者可能是在不同的方向上调控，以形成对比时，“PValue.Mixed”列是合适的。 +FDRs由所有集合中对应的p值计算得到。 -The results of a gene set test can be viewed in a barcode plot produced by the `barcodeplot` function. -Suppose visualization is performed for the gene set defined by the GO term GO:0032465: +基因集测试的结果可以在“barcodeplot”函数生成的条形码图中查看。 +假设对GO术语GO:0032465定义的基因集进行可视化: ```{r barcode, fig.width=8, fig.height=5, fig.cap="Barcode plot showing enrichment of the GO term GO:0032465 in the basal virgin group compared to the basal lactating group. X-axis shows logFC for B.virgin vs B.lactating. Black bars represent genes annotated with the GO term. The worm shows relative enrichment."} res <- glmQLFTest(fit, contrast=B.VvsL) @@ -684,30 +680,30 @@ barcodeplot(res$table$logFC, index=index, labels=c("B.lactating","B.virgin"), main=cyt.go[1]) ``` -In the plot, all genes are ranked from left to right by increasing log-fold change for the contrast and the genes within the gene set are represented by vertical bars, forming the barcode-like pattern. -The curve (or *worm*) above the barcode shows the relative local enrichment of the bars in each part of the plot. -The dotted horizontal line indicates neutral enrichment; the worm above the dotted line shows enrichment while the worm below the dotted line shows depletion. -In this particular barcode plot the worm shows enrichment on the right for positive logFCs, and depletion on the left for negative logFCs. -The conclusion is that genes associated with this GO term tend to be up-regulated in the basal cells of virgin mice compared to lactating mice, confirming the result of the `fry` test above. +在图中，通过增加对比的logFC，将所有基因从左到右排列，基因集中的基因用竖条表示，形成条形码状图案？。 +条形码上方的曲线(或*worm*)显示了图中每个部分的相对局部富集。 +虚线表示中性富集，虚线以上的worm表示富集，虚线以下的worm表示下降？。 +在这个特定的条形码图中，worm在右侧显示正logfc的富集，在左侧显示负logfc的depletion（耗尽？）。 +结论是，与泌乳小鼠相比，与GO相关的基因在初生小鼠的基底细胞中趋向于上调，证实了上述“fry”试验的结果。 -Camera gene set enrichment analysis {#camera-gene-set-enrichment-analysis .unnumbered} +Camera 基因富集分析 {#camera-gene-set-enrichment-analysis .unnumbered} ----------------------------------- -Finally we demonstrate a gene set enrichment style analysis using the Molecular Signatures Database (MSigDB) [@subramanian2005gsea]. -We will use the C2 collection of the MSigDB, which is a collection of nearly 5000 curated gene sets, each representing the molecular signature of a particular biological process or phenotype. -The MSigDB itself is purely human, but the Walter and Eliza Hall Institute (WEHI) maintains a mouse version of the database. -We load the mouse version of the C2 collection from the WEHI website: + 最后我们使用分子特征数据库（MSigDB）演示了基因集富集样式分析[@subramanian2005gsea]。 +我们将使用MSigDB的C2集合，这是一个由近5000个精心选择的基因集组成的集合，每个集合代表特定生物学过程或表型的分子特征。 +MSigDB本身是纯人类的，但是Walter and Eliza Hall Institute (WEHI)维护一个小鼠版本的数据库。 +我们加载来自WEHI网站c2数据集的小鼠样本： ```{r camera} load(url("http://bioinf.wehi.edu.au/software/MSigDB/mouse_c2_v5p1.rdata")) ``` -This will load `Mm.c2`, which is a list of gene sets, each a vector of Entrez Ids. -This can be converted to a list of index numbers: +加载“Mm.c2”，这是一个基因集列表，每个都是Entrez id的向量。 +转换为索引编号列表: ```{r id2indices.Mm.c2} idx <- ids2indices(Mm.c2,id=rownames(y)) ``` -First we compare virgin stem cells to virgin luminal cells: +首先我们比较未妊娠小鼠的干细胞和腔细胞： ```{r cam.BVvsLV} BvsL.v <- makeContrasts(B.virgin - L.virgin, levels=design) @@ -715,10 +711,10 @@ cam <- camera(y, idx, design, contrast=BvsL.v, inter.gene.cor=0.01) options(digits=2) head(cam,14) ``` -With a large gene set collection, setting `inter.gene.cor = 0.01` gives a good compromise between biological interpretability and FDR control. -As expected, the mammary stem cell and mammary luminal cell signatures from Lim et al [@lim2010transcriptome] are top-ranked, and in the expected directions. +对一个大的基因集收集，设置“inter.gene.cor=0.01”在生物学解释能力和FDR控制之间提供了一个很好的折衷。 +正如预期的那样，来自Lim等的乳腺干细胞和腔细胞都以预期的方向排名靠前[@lim2010transcriptome]。 -We can visualize the top signature, combining the up and down mammary stem cell signatures to make a bi-directional signature set: +我们可以看到顶部的标志，结合上下乳腺干细胞的标志，形成一个双向标志集: ```{r barcode2, fig.width=8, fig.height=6.4, fig.cap="Barcode plot showing strong enrichment of mammary stem cell signature in the stem cell vs luminal cell comparison. Red bars show up signature genes, blue bars show down genes. The worms show relative enrichment."} res <- glmQLFTest(fit, contrast=BvsL.v) @@ -731,52 +727,52 @@ barcodeplot(res$table$logFC, ``` -Packages used {#packages-used .unnumbered} +使用的包 {#packages-used .unnumbered} ============= -This workflow depends on various packages from version 3.3 of the Bioconductor project, running on R version 3.3.0 or higher. -The complete list of the packages used for this workflow are shown below: +本工作流依赖于Bioconductor项目3.3版的各种包，运行在R 3.3.0或更高版本上。 +该工作流程使用的包的完整列表如下: ```{r info} sessionInfo() ``` -Read alignment and quantification {#read-alignment-and-quantification .unnumbered} +Read 比对和量化 {#read-alignment-and-quantification .unnumbered} ================================= -Download raw sequence files from the SRA {#download-raw-sequence-files-from-the-sra .unnumbered} +从SRA下载原始序列文件 {#download-raw-sequence-files-from-the-sra .unnumbered} ---------------------------------------- -We now revisit the question of recreating the matrix of read counts from the raw sequence reads. -Unlike the above workflow, which works for any version of R, read alignment requires Unix or Mac OS and, in practice, a high performance Unix server is recommended. -Read alignment and read counting require only one Bioconductor package, *Rsubread*. -However the `fastq-dump` utility from the SRA Toolkit is also required to convert from SRA to FASTQ format. -This can be downloaded from the NCBI website () and installed on any Unix system. +现在我们重新讨论从原始序列读取中重新创建读取计数矩阵的问题。 +与上面的工作流(适用于任何版本的R)不同，read对齐需要Unix或Mac OS，实际上，建议使用高性能Unix服务器。 +Read比对和计数只需要一个 “Rsubread” Bioconductor包。 +并且，也需要来自SRA工具包的“fastq-dump”将SRA转换为FASTQ格式。 +可以从NCBI网站()下载，并安装在任何Unix系统上。 -The first task is to download the raw sequence files, which are stored in SRA format on the SRA repository. -The SRA files need to be unpacked into FASTQ format using the `fastq-dump` utility. -The following R code makes a system call to `fastq-dump` to download each SRA file and convert it to FASTQ format: +第一个任务是下载原始序列文件，这些文件以SRA格式存储在SRA存储库中。 +需要使用“FASTQ -dump”程序将SRA文件解压缩为FASTQ格式。 +以下R代码调用系统“FASTQ -dump”下载每个SRA文件并将其转换为FASTQ格式: for (sra in targets$SRA) system(paste("fastq-dump", sra)) -The `fastq-dump` utility automatically downloads the specified SRA data set from the internet. -The above code will produce 12 FASTQ files, in the current working directory, with file names given by the following vector: +“fastq-dump”程序自动从internet下载指定的SRA数据集。 +上述代码将在当前工作目录下生成12个FASTQ文件，文件名如下向量所示: all.fastq <- paste0(targets$SRA, ".fastq") -Accuracy of base-calling {#accuracy-of-base-calling .unnumbered} +base-calling 精度 {#accuracy-of-base-calling .unnumbered} ------------------------ -Sequencers typically store base-calling quality scores for each read in the FASTQ files. -*Rsubread*'s `qualityScores` function can be used to extract these scores from any particular file: +排序器通常在FASTQ文件中存储每次读取的 base-calling 质量分数。 +可以使用“Rsubread”包的“qualityScores”函数从任何形式的文件中提取这些分数： QS <- qualityScores("SRR1552444.fastq") -The boxplot function provides a compact way to view the quality scores by position across all reads: +boxplot函数提供了一种紧凑的方法来查看所有读取的位置的质量分数: boxplot(QS, ylab="Quality score", xlab="Base position", main="SRR1552444.fastq", cex=0.25, col="orange") @@ -785,986 +781,51 @@ The boxplot function provides a compact way to view the quality scores by positi
Boxplots of quality scores by base position for the first FASTQ file.

-Figure 12: Boxplots of quality scores by base position for the first FASTQ file.
+图 12: Boxplots of quality scores by base position for the first FASTQ file.

-The vertical axis shows the Phred quality score, equal to $-10\log_{10}(p)$ where $p$ is the probability of an erroneous call. -The maximum possible value is 40, and all values above 10 correspond to extremely small error probabilities. -The horizontal axis shows position within a read. -The file contains 100bp single-end reads, so the scale is from 1 to 100. -The plot displays a compact boxplot at each base position. -As is very commonly observed, the quality scores are best in the middle of the reads and decrease slightly towards the start and end of the reads. -However the quality remains generally good even near the ends of the reads: the scores would need to be very much lower than this before they would cause problems for the alignment. -Similar plots can be made for each of the FASTQ files. +纵轴显示Phred质量分数，等于$-10\log_{10}(p)$，其中$p$是错误调用的概率。 +最大可能值是40，所有大于10的值对应的错误概率都非常小。 +横轴显示读取中的位置。 +该文件包含100bp的单端读取，所以范围是从1到100。 +该图在每个基本位置显示一个紧凑的箱形图。 +正如我们经常观察到的，质量分数在read的中间最好，在read的开始和结束时略有下降。 +然而，即使在接近末端的reads，质量通常仍然很好:分数需要比这低得多，才会对对齐造成问题。 +可以为每个FASTQ文件绘制类似的图。 -Build a genome index {#build-a-genome-index .unnumbered} +建立基因组索引 {#build-a-genome-index .unnumbered} -------------------- -Before the sequence reads can be aligned, we need to build an index for the GRCm38/mm10 (Dec 2011) build of the mouse genome. -Most laboratories that use *Rsubread* regularly will already have an index file prepared, as this is a once-off operation for each genome release. -If you are using *Rsubread* for mouse for the first time, then the latest mouse genome build can be downloaded from the NCBI webpage -(). -The full URL is: +在序列读取对齐之前，我们需要为小鼠基因组构建的GRCm38/mm10(2011年12月)建立索引。 +大多数经常使用“Rsubread”包的实验室已经准备好了索引文件，因为这是每次基因组释放的一次性操作。 +若第一次使用针对小鼠的“Rsubread”包，可以从NCBI网站()下载最新的小时基因构建文件。 - ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz - -(Note that this link is for patch 6 of the RefSeq mm10 build, which is valid at the time of writing in May 2018. -The link will change as new patches are released periodically.) -An index can then be built by: - - library(Rsubread) - buildindex(basename = "mm10", - reference = "GCF_000001635.26_GRCm38.p6_genomic.fna.gz") - -Aligning reads {#aligning-reads .unnumbered} --------------- - -The sequence reads can now be aligned to the mouse genome using the `align` function: - - all.bam <- sub(".fastq", ".bam", all.fastq) - align(index="mm10", readfile1=all.fastq, input_format="FASTQ", - output_file=all.bam) - -This produces a set of BAM files containing the read alignments for each RNA library. -The mapping proportions can be summarized by the `propmapped` function: - - > propmapped(all.bam) - Samples NumTotal NumMapped PropMapped - 1 SRR1552450.bam 30109290 26577308 0.883 - 2 SRR1552451.bam 28322351 24794251 0.875 - 3 SRR1552452.bam 31688348 27937620 0.882 - 4 SRR1552453.bam 29614284 26074034 0.880 - 5 SRR1552454.bam 27225012 24381742 0.896 - 6 SRR1552455.bam 25433157 22813815 0.897 - 7 SRR1552444.bam 27919481 23927833 0.857 - 8 SRR1552445.bam 29731031 25487822 0.857 - 9 SRR1552446.bam 29879070 25500318 0.853 - 10 SRR1552447.bam 29245388 25187577 0.861 - 11 SRR1552448.bam 31425424 27326500 0.870 - 12 SRR1552449.bam 31276061 27204156 0.870 - -Ideally, the proportion of mapped reads should be above 80%. -By default, only reads with unique mapping locations are reported by *Rsubread* as being successfully mapped. -Restricting to uniquely mapped reads is recommended, as it avoids spurious signal from non-uniquely mapped reads derived from, e.g., repeat regions. - - -Quantifying read counts for each gene {#quantifying-read-counts-for-each-gene .unnumbered} -------------------------------------- - -The read counts for each gene can be quantified using the `featureCounts` function in *Rsubread*. -Conveniently, the *Rsubread* package includes inbuilt NCBI RefSeq annotation of the mouse and human genomes. -`featureCounts` generates a matrix of read counts for each gene in each sample: - - fc <- featureCounts(all.bam, annot.inbuilt="mm10") - -The output is a simple list, containing the matrix of counts (`counts`), a data frame of gene characteristics (`annotation`), a vector of file names (`targets`) and summary mapping statistics (`stat`): - - > names(fc) - [1] "counts" "annotation" "targets" "stat" - -The row names of `fc$counts` are the Entrez gene identifiers for each gene. -The column names are the output file names from `align`, which we simplify here for brevity: - - > colnames(fc$counts) <- rownames(targets) - -The first six rows of the counts matrix are shown below. - - > head(fc$counts) - MCL1.DG MCL1.DH MCL1.DI MCL1.DJ MCL1.DK MCL1.DL MCL1.LA - 497097 438 299 65 237 354 287 0 - 100503874 1 0 1 1 0 4 0 - 100038431 0 0 0 0 0 0 0 - 19888 1 1 0 0 0 0 10 - 20671 106 181 82 104 43 83 16 - 27395 309 232 339 290 291 270 558 - MCL1.LB MCL1.LC MCL1.LD MCL1.LE MCL1.LF - 497097 0 0 0 0 0 - 100503874 0 0 0 0 0 - 100038431 0 0 0 0 0 - 19888 3 10 2 0 0 - 20671 25 18 8 3 10 - 27395 468 488 332 312 344 - -Finally, a `DGEList` object can be assembled by: - - y <- DGEList(fc$counts, group=group) - y$genes <- fc$annotation[, "Length", drop=FALSE] - - -Data and software availability {#data-and-software-availability .unnumbered} -============================== - -This html document was generated as a vignette of the RnaSeqGeneEdgeRQL package, which can be installed by running - - install.packages("BiocManager") - BiocManager::install("RnaSeqGeneEdgeRQL") - -at the R prompt. -Except for the targets file `targets.txt`, all data analyzed in the workflow is read automatically from public websites as part of the code. -All software used is publicly available as part of Bioconductor (Release 3.3 or later), except for the `fastq-dump` utility, which can be downloaded from NCBI website as described in the text. -The article includes the complete code necessary to reproduce the analyses shown. - -The original Bioconductor 3.3 version of this workflow, including a pdf article version and a corresponding Rnw file, is available from and [@chen2016quasiworkflow]. - - -Author contributions {#author-contributions .unnumbered} -==================== - -All authors developed and tested the code workflow. -All authors wrote the article. - -Grant information {#grant-information .unnumbered} -================= - -This work was supported by the National Health and Medical Research Council -(Fellowship 1058892 and Program 1054618 to G.K.S, Independent Research Institutes Infrastructure Support to the Walter and Eliza Hall Institute) -and by a Victorian State Government Operational Infrastructure Support Grant. - -Acknowledgments {#acknowledgments .unnumbered} -=============== - -The authors thank Wei Shi and Yang Liao for advice with *Rsubread* and Yifang Hu for creating the mouse version of the MSigDB. -We also wish to acknowledge the early developers of the quasi-likelihood testing approach. -Davis McCarthy authored the first versions of the quasi-likelihood functions in *edgeR* in February 2011 and Steve Lund and Dan Nettleton worked independently on *QuasiSeq* around the same time. -We thank Steve Lund and Dan Nettleton for a valuable and enjoyable collaboration that lead to Lund et al [@lund2012quasiseq]. - -How to cite {#how-to-cite .unnumbered} -==================== - -Please cite this workflow article as - -Chen, Yunshun, A. T. L. Lun, and G. K. Smyth. 2016. "From Reads to Genes to Pathways: Differential Expression Analysis of RNA-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline." *F1000Research* 5:1438. - -References {#references .unnumbered} -========== - -lter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia - - name: Aaron T. L. Lun - affiliation: Cancer Research UK Cambridge Institute, University of Cambridge, Li Ka Shing Centre, Robinson Way, CB20RE, Cambridge, United Kingdom - - name: Gordon K. Smyth - affiliation: The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia -date: 28 Oct 2018 -vignette: > - %\VignetteIndexEntry{From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} -bibliography: ref.bib -output: - BiocStyle::html_document: - fig_caption: true ---- -**R version**: `r R.version.string` - -**Bioconductor version**: `r BiocManager::version()` - -**Package**: `r packageVersion("RnaSeqGeneEdgeRQL")` - -Abstract {#abstract .unnumbered} -======== - -RNA sequencing (RNA-seq) has become a very widely used technology for profiling gene expression. One of the most common aims of RNA-seq profiling is to identify genes or molecular pathways that are differentially expressed (DE) between two or more biological conditions. This article demonstrates a computational workflow for the detection of DE genes and pathways from RNA-seq data by providing a complete analysis of an RNA-seq experiment profiling epithelial cell subsets in the mouse mammary gland. The workflow uses R software packages from the open-source Bioconductor project and covers all steps of the analysis pipeline, including alignment of read sequences, data exploration, differential expression analysis, visualization and pathway analysis. Read alignment and count quantification is conducted using the Rsubread package and the statistical analyses are performed using the edgeR package. The differential expression analysis uses the quasi-likelihood functionality of edgeR. - - -Introduction {#introduction .unnumbered} -============ - -RNA sequencing (RNA-seq) has become a very widely used technology for profiling transcriptional activity in biological systems. -One of the most common aims of RNA-seq profiling is to identify genes or molecular pathways that are differentially expressed (DE) between two or more biological conditions. -Changes in expression can then be associated with differences in biology, providing avenues for further investigation into potential mechanisms of action. - -This article provides a detailed workflow for analyzing an RNA-seq study all the way from the raw reads through to differential expression and pathway analysis using Bioconductor packages [@huber2015orchestrating]. -The article gives a complete analysis of RNA-seq data that were collected to study the effects of pregnancy and lactation on the luminal cell lineage in the mouse mammary gland [@fu2015egfmcl1]. -The pipeline uses the *Rsubread* package [@liao2013subread] for mapping reads and assigning them to genes and the *edgeR* package [@robinson2010edger] for statistical analyses. - -A journal article version of this workflow, including referee reports and responses, is available at [F1000Research](http://f1000research.com/articles/5-1438) [@chen2016quasiworkflow]. - -RNA-seq analysis involves a number of steps, including read alignment, read summarization, differential expression and pathway analysis. -Here we use the Subread aligner [@liao2013subread] for mapping and featureCounts [@liao2014featurecounts] for assigning reads to genes. -As well as being fast and efficient, these algorithms have the advantage of having native implementations as R functions in the *Rsubread* package. -This means that the entire analysis can be conducted efficiently within the R environment. - -The workflow uses *edgeR*'s quasi-likelihood pipeline (edgeR-quasi) for differential expression. -This statistical methodology uses negative binomial generalized linear models [@mccarthy2012edgerglm] but with F-tests instead of likelihood ratio tests [@lund2012quasiseq]. -This method provides stricter error rate control than other negative binomial based pipelines, including the traditional *edgeR* pipelines [@robinson2008smallsample; @robinson2007moderated; @mccarthy2012edgerglm] or *DESeq2* [@love2014deseq2]. -The edgeR-quasi pipeline is based on similar statistical methodology to that of the *QuasiSeq* package [@lund2012quasiseq], which has performed well in third-party comparisons [@burden2014error]. -Compared to *QuasiSeq*, the *edgeR* functions offer speed improvements and some additional statistical refinements [@lun2016delicious]. -The RNA-seq pipelines of the *limma* package also offer excellent error rate control [@law2014voom; @ritchie2015limma]. -While the *limma* pipelines are recommended for large-scale datasets, because of their speed and flexibility, the edgeR-quasi pipeline gives better performance in low-count situations [@lun2014denovo; @lun2015diffhic]. -For the data analyzed here, the edgeR-quasi, limma-voom and limma-trend pipelines are all equally suitable and give similar results. - -The analysis approach illustrated in this article can be applied to any RNA-seq study that includes some replication, but it is especially appropriate for designed experiments with multiple treatment factors and with small numbers of biological replicates. -The approach assumes that RNA samples have been extracted from cells of interest under two or more treatment conditions, that RNA-seq profiling has been applied to each RNA sample and that there are independent biological replicates for at least one of the treatment conditions. -The *Rsubread* part of the workflow takes FASTQ files of raw sequence reads as input, while the *edgeR* part of the pipeline takes a matrix of genewise read counts as input. - -```{r setup, echo=FALSE, results="hide", warning=FALSE} -knitr::opts_chunk$set(prompt = TRUE, comment = NA) -suppressPackageStartupMessages({library(RnaSeqGeneEdgeRQL)}) -``` - - -Description of the biological experiment {#description-of-the-biological-experiment .unnumbered} -======================================== - -This workflow demonstrates a complete bioinformatics analysis of an RNA-seq study that is available from the GEO repository as series GSE60450. -The RNA-seq data were collected to study the lineage of luminal cells in the mouse mammary gland and in particular how the expression profiles of the members of the lineage change upon pregnancy and lactation [@fu2015egfmcl1]. -Specifically, the study examined the expression profiles of basal stem-cell enriched cells (B) and committed luminal cells (L) in the mammary glands of virgin, pregnant and lactating mice. -There are therefore six groups of RNA samples, one for each combination of cell type and mouse status. -Two biological replicates were collected for each group. - -This study used an Illumina Hiseq sequencer to generate about 30 million 100bp single-end reads for each sample. -Subread version 1.4.4 () was used to align the reads to the mouse mm10 genome and featureCounts was used to assign reads to Entrez Genes using RefSeq gene annotation. -The FASTQ files containing the raw sequence reads were deposited to the Sequence Read Archive (SRA) repository and the read counts were deposited to GEO. - -This experimental design is summarized in the table below, where the basal and luminal cell types are abbreviated with `B` and `L` respectively. -The GEO and SRA identifiers for each RNA sample are also shown: -```{r sra} -library(RnaSeqGeneEdgeRQL) -targetsFile <- system.file("extdata", "targets.txt", - package="RnaSeqGeneEdgeRQL") -targets <- read.delim(targetsFile, stringsAsFactors=FALSE) -targets -``` -The experiment can be viewed as a one-way layout with six groups. -For later use, we combine the treatment factors into a single grouping factor: - -```{r group} -group <- paste(targets$CellType, targets$Status, sep=".") -group <- factor(group) -table(group) -``` - -**Note**. -This study isolated carefully sorted cell populations obtained from genetically indentical mice under controlled laboratory conditions. -As will be shown during the analysis below, the relatively small sample size ($n=2$ in each group) here is justified by the low background variability and by the fact that the expression profiles of the different groups are distinctly different. -Studies with more variability (for example on human subjects) or with smaller effect sizes may well require more biological replicates for reliable results. - - -Preliminary analysis {#preliminary-analysis .unnumbered} -==================== - -Downloading the read counts {#downloading-the-read-counts .unnumbered} ---------------------------- - -Readers wishing to reproduce the analysis presented in this article can either download the matrix of read counts from GEO or recreate the read count matrix from the raw sequence counts. -We will present first the analysis using the downloaded matrix of counts. -At the end of this article we will present the R commands needed to recreate this matrix. - -The following commands download the genewise read counts for the GEO series GSE60450. -The zipped tab-delimited text file `GSE60450_Lactation-GenewiseCounts.txt.gz` will be downloaded to the working R directory: - -```{r checkdownload, echo=FALSE, results="hide", message=FALSE} -if( !file.exists("GSE60450_Lactation-GenewiseCounts.txt.gz") ) { -FileURL <- paste( - "http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE60450", - "format=file", - "file=GSE60450_Lactation-GenewiseCounts.txt.gz", - sep="&") -download.file(FileURL, method="libcurl", "GSE60450_Lactation-GenewiseCounts.txt.gz") -} -``` - -```{r download, eval=FALSE} -FileURL <- paste( - "http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE60450", - "format=file", - "file=GSE60450_Lactation-GenewiseCounts.txt.gz", - sep="&") -download.file(FileURL, "GSE60450_Lactation-GenewiseCounts.txt.gz") -``` - -The counts can then be read into a data.frame in R: - -```{r readcounts} -GenewiseCounts <- read.delim("GSE60450_Lactation-GenewiseCounts.txt.gz", - row.names="EntrezGeneID") -colnames(GenewiseCounts) <- substring(colnames(GenewiseCounts),1,7) -dim(GenewiseCounts) -head(GenewiseCounts) -``` -The row names of `GenewiseCounts` are the Entrez Gene Identifiers. -The first column contains the length of each gene, being the total number of bases in exons and UTRs for that gene. -The remaining 12 columns contain read counts and correspond to rows of `targets`. - - -The *edgeR* package stores data in a simple list-based data object called a `DGEList`. -This object is easy to use as it can be manipulated like an ordinary list in R, and it can also be subsetted like a matrix. -The main components of a `DGEList` object are a matrix of read counts, sample information in the `data.frame` format and optional gene annotation. -We enter the counts into a `DGEList` object using the function `DGEList` in *edgeR*: - -```{r DGEList, message=FALSE} -library(edgeR) -y <- DGEList(GenewiseCounts[,-1], group=group, - genes=GenewiseCounts[,1,drop=FALSE]) -options(digits=3) -y$samples -``` - -Adding gene annotation {#adding-gene-annotation .unnumbered} ----------------------- - -The Entrez Gene Ids link to gene information in the NCBI database. -The *org.Mm.eg.db* package can be used to complement the gene annotation information. -Here, a column of gene symbols is added to `y$genes`: - -```{r symbols, message=FALSE} -library(org.Mm.eg.db) -y$genes$Symbol <- mapIds(org.Mm.eg.db, rownames(y), - keytype="ENTREZID", column="SYMBOL") -head(y$genes) -``` -Entrez Ids that no longer have official gene symbols are dropped from the analysis. -The whole `DGEList` object, including annotation as well as counts, can be subsetted by rows as if it was a matrix: - -```{r dropNAsymbols} -y <- y[!is.na(y$genes$Symbol), ] -dim(y) -``` - - -Design matrix {#design-matrix .unnumbered} -------------- - -Linear modeling and differential expression analysis in *edgeR* requires a design matrix to be specified. -The design matrix records which treatment conditions were applied to each samples, and it also defines how the experimental effects are parametrized in the linear models. -The experimental design for this study can be viewed as a one-way layout and the design matrix can be constructed in a simple and intuitive way by: - -```{r design} -design <- model.matrix(~0+group) -colnames(design) <- levels(group) -design -``` - -This design matrix simply links each group to the samples that belong to it. -Each row of the design matrix corresponds to a sample whereas each column represents a coefficient corresponding to one of the six groups. - - -Filtering to remove low counts {#filtering-to-remove-low-counts .unnumbered} ------------------------------- - -Genes that have very low counts across all the libraries should be removed prior to downstream analysis. -This is justified on both biological and statistical grounds. -From biological point of view, a gene must be expressed at some minimal level before it is likely to be translated into a protein or to be considered biologically important. -From a statistical point of view, genes with consistently low counts are very unlikely be assessed as significantly DE because low counts do not provide enough statistical evidence for a reliable judgement to be made. -Such genes can therefore be removed from the analysis without any loss of information. - -As a rule of thumb, we require that a gene have a count of at least 10 or so in at least some libraries before it is considered to be expressed in the study. -In practice, the filtering is actually based on count-per-million (CPM) values so as to avoid favoring genes that are expressed in larger libraries over those expressed in smaller libraries. -In *edgeR*, the filtering can be accomplished using a function that takes into account the library sizes and the experimental design: -```{r keep} -keep <- filterByExpr(y, design) -table(keep) -``` -The `DGEList` object is then subsetted to retain only the non-filtered genes: -```{r filter} -y <- y[keep, , keep.lib.sizes=FALSE] -``` -The option `keep.lib.sizes=FALSE` causes the library sizes to be recomputed after the filtering. -This is generally recommended, although the effect on the downstream analysis is usually small. - -The `filterByExpr` function attempts to keep the maximum number of interesting genes in the analysis, but other sensible filtering criteria are also possible. -For example `keep <- rowSums(y$counts) > 50` is a very simple criterion that would keep genes with a total read count of more than 50. -This would give similar downstream results for this dataset to the filtering actually used. - -Another sensible filter rule would be to compute the average-log-CPM for each gene, and to choose a cutoff value heuristically: -```{r aveLogCPM, fig.cap="Histogram of average log2 CPM."} -AveLogCPM <- aveLogCPM(y) -hist(AveLogCPM) -``` -Again this might give similar downstream results for this data to the filtering actually used. - -Whatever the filtering rule, it should be independent of the information in the targets file. -It should not make any reference to which RNA libraries belong to which group, because doing so would bias the subsequent differential expression analysis. - - -Normalization for composition bias {#normalization-for-composition-bias .unnumbered} ----------------------------------- - -Normalization by trimmed mean of M values (TMM) [@robinson2010tmm] is performed by using the `calcNormFactors` function, which returns the `DGEList` argument with only the `norm.factors` changed. -It calculates a set of normalization factors, one for each sample, to eliminate composition biases between libraries. -The product of these factors and the library sizes defines the effective library size, which replaces the original library size in all downstream analyses. - -```{r norm} -y <- calcNormFactors(y) -y$samples -``` - -The normalization factors of all the libraries multiply to unity. -A normalization factor below one indicates that a small number of high count genes are monopolizing the sequencing, causing the counts for other genes to be lower than would be usual given the library size. -As a result, the effective library size will be scaled down for that sample. -Here we see that the luminal-lactating samples have low normalization factors. -This is a sign that these samples contain a number of very highly upregulated genes. - -**Note**. -In general, we find TMM normalization to be satisfactory for almost all well-designed mRNA gene expression experiments. -Single-cell RNA-seq is an exception, for which specialized normalization methods are needed [@lun2016pooling]. -Another, less common, type of study requiring special treatment is that with global differential expression, with more than half of the genome differentially expressed between experimental conditions in the same direction [@wu2013mirna]. -Global differential expression should generally be avoided in well designed experiments. -When it can't be avoided, then some normalization reference such as spike-ins needs to be built into the experiment for reliable normalization to be done [@risso2014ruvseq]. - - - -Exploring differences between libraries {#exploring-differences-between-libraries .unnumbered} ---------------------------------------- - -The RNA samples can be clustered in two dimensions using multi-dimensional scaling (MDS) plots. -This is both an analysis step and a quality control step to explore the overall differences between the expression profiles of the different samples. -Here we decorate the MDS plot to indicate the cell groups: - -```{r mdsplot, fig.width=7, fig.height=7, fig.cap="MDS plot showing distances between expression profiles. This provides a type of unsupervised clustering of the samples. In this case, the first dimension separates samples by cell type while the second dimension corresponds to mouse status."} -pch <- c(0,1,2,15,16,17) -colors <- rep(c("darkgreen", "red", "blue"), 2) -plotMDS(y, col=colors[group], pch=pch[group]) -legend("topleft", legend=levels(group), pch=pch, col=colors, ncol=2) -``` - -In the MDS plot, the distance between each pair of samples can be interpreted as the leading log-fold change between the samples for the genes that best distinguish that pair of samples. -By default, leading fold-change is defined as the root-mean-square of the largest 500 log2-fold changes between that pair of samples. -The figure above shows that replicate samples from the same group cluster together while samples from different groups are well separated. -In other words, differences between groups are much larger than those within groups, meaning that there are likely to be statistically significant differences between the groups. -The distance between basal cells on the left and luminal cells on the right is about six units on the x-axis, corresponding to a leading fold change of about 64-fold between the two cell types. -The differences between the virgin, pregnant and lactating expression profiles appear to be magnified in luminal cells compared to basal. - -The expression profiles of individual samples can be explored more closely with mean-difference (MD) plots. -An MD plot visualizes the library size-adjusted log-fold change between two libraries (the difference) against the average log-expression across those libraries (the mean). -The following command produces an MD plot that compares sample 1 to an artificial reference library constructed from the average of all the other samples: -```{r mdplot1, fig.width=7, fig.height=7, fig.cap="Mean-difference plot of log2-expression in sample 1 versus the average log2-expression across all other samples. Each point represents a gene, and the red line indicates a log-ratio of zero. The majority of points cluster around the red line."} -plotMD(y, column=1) -abline(h=0, col="red", lty=2, lwd=2) -``` -The bulk of the genes are centered around the line of zero log-fold change. -The diagonal lines in the lower left of the plot correspond to genes with counts of 0, 1, 2 and so on in the first sample. - -It is good practice to make MD plots for all the samples as a quality check. -We now look at one of the luminal-lactating samples that was observed have a low normalization factor: -```{r mdplot11, fig.width=7, fig.height=7, fig.cap="Mean-difference plot of log2-expression in sample 11 versus the average log2-expression across all other samples. The plot shows a number of genes that are both highly expressed and highly up-regulated."} -plotMD(y, column=11) -abline(h=0, col="red", lty=2, lwd=2) -``` -For this sample, the log-ratios show noticeable positive skew, with a number of very highly upregulated genes. -In particular, there are a number of points in the upper right of the plot, corresponding to genes that are both highly expressed and highly up-regulated in this sample compared to others. -These genes explain why the normalization factor for this sample is well below one. -By contrast, the log-ratios for sample 1 were somewhat negatively skewed, corresponding to a normalization factor above one. - - -Dispersion estimation {#dispersion-estimation .unnumbered} ---------------------- - -*edgeR* uses the negative binomial (NB) distribution to model the read counts for each gene in each sample. -The dispersion parameter of the NB distribution accounts for variability between biological replicates [@mccarthy2012edgerglm]. -*edgeR* estimates an empirical Bayes moderated dispersion for each individual gene. -It also estimates a common dispersion, which is a global dispersion estimate averaged over all genes, and a trended dispersion where the dispersion of a gene is predicted from its abundance. -Dispersion estimates are most easily obtained from the `estimateDisp` function: - -```{r estimateDisp} -y <- estimateDisp(y, design, robust=TRUE) -``` -This returns a `DGEList` object with additional components (`common.dispersion`, `trended.dispersion` and `tagwise.dispersion`) added to hold the estimated dispersions. -Here `robust=TRUE` has been used to protect the empirical Bayes estimates against the possibility of outlier genes with exceptionally large or small individual dispersions [@phipson2016robust]. - -The dispersion estimates can be visualized with `plotBCV`: -```{r plotBCV, fig.width=7, fig.height=7, fig.cap="Scatterplot of the biological coefficient of variation (BCV) against the average abundance of each gene. The plot shows the square-root estimates of the common, trended and tagwise NB dispersions."} -plotBCV(y) -``` - -The vertical axis of the `plotBCV` plot shows square-root dispersion, also known as *biological coefficient of variation* (BCV) [@mccarthy2012edgerglm]. -For RNA-seq studies, the NB dispersions tend to be higher for genes with very low counts. -The dispersion trend tends to decrease smoothly with abundance and to asymptotic to a constant value for genes with larger counts. -From our past experience, the asymptotic value for the BCV tends to be in range from 0.05 to 0.2 for genetically identical mice or cell lines, whereas somewhat larger values ($>0.3$) are observed for human subjects. - -The NB model can be extended with quasi-likelihood (QL) methods to account for gene-specific variability from both biological and technical sources [@lund2012quasiseq; @lun2016delicious]. -Under the QL framework, the NB dispersion trend is used to describe the overall biological variability across all genes, and gene-specific variability above and below the overall level is picked up by the QL dispersion. -In the QL approach, the individual (tagwise) NB dispersions are not used. - -The estimation of QL dispersions is performed using the `glmQLFit` function: -```{r glmQLFit} -fit <- glmQLFit(y, design, robust=TRUE) -head(fit$coefficients) -``` -This returns a `DGEGLM` object with the estimated values of the GLM coefficients for each gene. -It also contains a number of empirical Bayes (EB) statistics including the QL dispersion trend, the squeezed QL dispersion estimates and the prior degrees of freedom (df). -The QL dispersions can be visualized by `plotQLDisp`: - -```{r QLDisp, fig.width=7, fig.height=7, fig.cap="A plot of the quarter-root QL dispersion against the average abundance of each gene. Estimates are shown for the raw (before EB moderation), trended and squeezed (after EB moderation) dispersions. Note that the QL dispersions and trend shown here are relative to the NB dispersion trend show in the previous figure produced using `plotBCV()`."} -plotQLDisp(fit) -``` - -The QL functions moderate the genewise QL dispersion estimates in the same way that the *limma* package moderates variances [@smyth2004ebayes]. -The raw QL dispersion estimates are squeezed towards a global trend, and this moderation reduces the uncertainty of the estimates and improves testing power. -The extent of the squeezing is governed by the value of the prior df estimated from the data. -Large prior df estimates indicate that the QL dispersions are less variable between genes, meaning that strong EB moderation should be performed. -Smaller prior df estimates indicate that the true unknown dispersions are highly variable, so weaker moderation towards the trend is appropriate. - -```{r df.prior} -summary(fit$df.prior) -``` - -Setting `robust=TRUE` in `glmQLFit` is usually recommended [@phipson2016robust]. -This allows gene-specific prior df estimates, with lower values for outlier genes and higher values for the main body of genes. -This reduces the chance of getting false positives from genes with extremely high or low raw dispersions, while at the same time increasing statistical power to detect differential expression for the main body of genes. - -Differential expression analysis {#differential-expression-analysis .unnumbered} -================================ - -Testing for differential expression {#testing-for-differential-expression .unnumbered} ------------------------------------ - -The next step is to test for differential expression between the experimental groups. -One of the most interesting comparisons is that between the basal pregnant and lactating groups. -The contrast corresponding to any specified comparison can be constructed conveniently using the `makeContrasts` function: - -```{r B.PvsL} -B.LvsP <- makeContrasts(B.lactating-B.pregnant, levels=design) -``` - -In subsequent results, a positive $\log_2$-fold-change (logFC) will indicate a gene up-regulated in lactating mice relative to pregnant, whereas a negative logFC will indicate a gene more highly expressed in pregnant mice. - -We will use QL F-tests instead of the more usual likelihood ratio tests (LRT) as they give stricter error rate control by accounting for the uncertainty in dispersion estimation: - -```{r glmQLFTest} -res <- glmQLFTest(fit, contrast=B.LvsP) -``` - -The top DE genes can be viewed with `topTags`: - -```{r topTags} -topTags(res) -``` - -In order to control the false discovery rate (FDR), multiple testing correction is performed using the Benjamini-Hochberg method. -The top DE gene *Csn1s2b* has a large positive logFC, showing that it is far more highly expressed in the basal cells of lactating than pregnant mice. -This gene is indeed known to be a major source of protein in milk. - -The total number of DE genes identified at an FDR of 5% can be shown with `decideTestsDGE`. -There are in fact more than 5000 DE genes in this comparison: - -```{r decideTests} -is.de <- decideTestsDGE(res) -summary(is.de) -``` - -The magnitude of the differential expression changes can be visualized with a fitted model MD plot: - -```{r plotMDfit, fig.width=7, fig.height=7, fig.cap="MD plot showing the log-fold change and average abundance of each gene. Significantly up and down DE genes are highlighted in red and blue, respectively."} -plotMD(res, status=is.de, values=c(1,-1), col=c("red","blue"), - legend="topright") -``` -The logFC for each gene is plotted against the average abundance in log2-CPM, i.e., `logCPM` in the table above. -Genes that are significantly DE are highlighted. - - -Differential expression above a fold-change threshold {#differential-expression-above-a-fold-change-threshold .unnumbered} ------------------------------------------------------ - -`glmQLFTest` identifies differential expression based on statistical significance regardless of how small the difference might be. -For some purposes we might be interested only in genes with reasonably large expression changes. -The above analysis found more than 5000 DE genes between the basal pregnant and lactating groups. -With such a large number of DE genes, it makes sense to narrow down the list to genes that are more biologically meaningful. - -A commonly used approach is to apply FDR and logFC cutoffs simultaneously. -However this tends to favor lowly expressed genes, and also fails to control the FDR correctly. -A better and more rigorous approach is to modify the statistical test so as to detect expression changes greater than a specified threshold. -In *edgeR*, this can be done using the `glmTreat` function. -This function is analogous to the TREAT method for microarrays [@mccarthy2009treat] but is adapted to the NB framework. -Here we test whether the differential expression fold changes are significantly greater than 1.5, that is, whether the logFCs are significantly greater than $\log_2(1.5)$: - -```{r treat} -tr <- glmTreat(fit, contrast=B.LvsP, lfc=log2(1.5)) -topTags(tr) -``` -Note that the argument \code{lfc} is an abbreviation for ``log-fold-change". -About 1100 genes are detected as DE with a FC significantly above 1.5 at an FDR cut-off of 5%. - - -```{r treatdecideTests} -is.de <- decideTestsDGE(tr) -summary(is.de) -``` -The p-values from `glmTreat` are larger than those from `glmQLFTest`, and the number of significantly DE genes is fewer, because it is testing an interval null hypothesis and requires stronger evidence for differential expression than does a conventional test. -It provides greater specificity for identifying the most important genes with large fold changes. - -The test results can be visualized in an MD plot: - -```{r plotMDtreat, fig.width=7, fig.height=7, fig.cap="MD plot showing the log-fold change and average abundance of each gene. Genes with fold-changes significantly greater than 1.5 are highlighted."} -plotMD(tr, status=is.de, values=c(1,-1), col=c("red","blue"), - legend="topright") -``` - -The \code{glmTreat} method evaluates variability as well as the magnitude of change of expression values and therefore is not equivalent to a simple fold change cutoff. -Nevertheless, all the statistically significant expression changes have logFC greater than 0.8 and almost all (97%) are greater than 0.9. -These values compare to the threshold value of $\log_2(1.5) = 0.58$. -In general, an estimated logFC must exceed the TREAT threshold by a number of standard errors for it to be called significant. -In other words, the whole confidence interval for the logFC must clear the threshold rather than just the estimated value itself. -It is better to interpret the threshold as *the FC below which we are definitely not interested in the gene* rather than *the FC above which we are interested in the gene*. - -The value of the FC threshold can be varied depending on the dataset. -In the presence of a huge number of DE genes, a relatively large FC threshold may be appropriate to narrow down the search to genes of interest. -In the absence of DE genes, on the other hand, a small or even no FC threshold shall be used. -If the threshold level is set to zero, then `glmTreat` becomes equivalent to `glmQLFTest` in the workflow shown here. - -In general, using `glmTreat` to reduce the number of DE genes is better than simply reducing the FDR cutoff, because `glmTreat` prioritizes genes with larger changes that are likely to be more biologically significant. -`glmTreat` can also be used with *edgeR* pipelines other than quasi-likelihood, although we don't demonstrate that here. - - -Heat map clustering {#heat-map-clustering .unnumbered} -------------------- - -Heatmaps are a popular way to display differential expression results for publication purposes. -To create a heatmap, we first convert the read counts into log2-counts-per-million (logCPM) values. -This can be done with the `cpm` function: - -```{r cpm} -logCPM <- cpm(y, prior.count=2, log=TRUE) -rownames(logCPM) <- y$genes$Symbol -colnames(logCPM) <- paste(y$samples$group, 1:2, sep="-") -``` -The introduction of `prior.count` is to avoid undefined values and to reduce the variability of the logCPM values for genes with low counts. -Larger values for `prior.count` shrink the logFCs for low count genes towards zero. - -We will create a heatmap to visualize the top 30 DE genes according to the TREAT test between `B.lactating` and `B.pregnant`. -The advantage of a heatmap is that it can display the expression pattern of the genes across all the samples. -Visualization of the results is aided by clustering together genes that have correlated expression patterns. -First we select the logCPM values for the 30 top genes: - -```{r order} -o <- order(tr$table$PValue) -logCPM <- logCPM[o[1:30],] -``` - -A heat map can then be produced by the `coolmap` function: -```{r heatmap, message=FALSE, fig.width=8, fig.height=12, fig.cap="Heat map across all the samples using the top 30 most DE genes between the basal lactating group and the basal pregnancy group."} -coolmap(logCPM, margins=c(7,7), lhei=c(1,6), lwid=c(1,3)) -``` -By default, `coolmap` clusters genes by correlation (highly correlated genes are closest) and clusters samples based on Euclidean distance between the expression values. -The coolmap has placed the basal samples on the right and the luminal samples on the left. -In each case, the virgin and pregnant samples are more like each other than the lactacting samples. -As expected, the two replicate samples from each group are clustered together. - -There a number of other R functions that can create heatmaps, but `coolmap` is the simplest because the default clustering metrics are chosen to be appropriate for expression data. - - -Analysis of deviance {#analysis-of-deviance .unnumbered} --------------------- - -The differential expression analysis comparing two groups can be easily extended to comparisons between three or more groups. -This is done by creating a matrix of independent contrasts. -In this manner, users can perform a one-way analysis of deviance (ANODEV) for each gene [@mccullagh1989glms]. - -Suppose we want to compare the three groups in the luminal population, i.e., virgin, pregnant and lactating. -An appropriate contrast matrix can be created as shown below, to make pairwise comparisons between all three groups: - -```{r makeContrasts} -con <- makeContrasts( - L.PvsL = L.pregnant - L.lactating, - L.VvsL = L.virgin - L.lactating, - L.VvsP = L.virgin - L.pregnant, levels=design) -``` - -The QL F-test is then applied to identify genes that are DE between the three groups. -This combines the three pairwise comparisons into a single F-statistic and p-value. -The top set of significant genes can be displayed with `topTags`: - -```{r anovaQLFtest} -res <- glmQLFTest(fit, contrast=con) -topTags(res) -``` - -Note that the three contrasts of pairwise comparisons are linearly dependent. -Constructing the contrast matrix with any two of the contrasts would be sufficient for an ANODEV test. -If the contrast matrix contains all three possible pairwise comparisons, then only the log-fold changes of the first two contrasts are shown in the output of `topTags`. - -Complicated contrasts {#complicated-contrasts .unnumbered} ---------------------- - -The flexibility of the GLM framework makes it possible to specify arbitrary contrasts for differential expression tests. -Suppose we are interested in testing whether the change in expression between lactating and pregnant mice is the same for basal cells as it is for luminal cells. -In statistical terminology, this is the interaction effect between mouse status and cell type. -The contrast corresponding to this testing hypothesis can be made as follows. - -```{r complicatedContrasts} -con <- makeContrasts( - (L.lactating-L.pregnant)-(B.lactating-B.pregnant), - levels=design) -``` - -Then the QL F-test is conducted to identify genes that are DE under this contrast. -The top set of DE genes are viewed with `topTags`. - -```{r complicatedQLTest} -res <- glmQLFTest(fit, contrast=con) -topTags(res) -``` - -Pathway analysis {#pathway-analysis .unnumbered} -================ - -Gene ontology analysis {#gene-ontology-analysis .unnumbered} ----------------------- - - -We now consider the problem of interpreting the differential expression results in terms of higher order biological processes or molecular pathways. -One of the most common used resources is gene ontology (GO) databases, which annotate genes according to a dictionary of annotation terms. -A simple and often effective way to interpret the list of DE genes is to count the number of DE genes that are annotated with each possible GO term. -GO terms that occur frequently in the list of DE genes are said to be over-represented or enriched. - -In *edgeR*, GO analyses can be conveniently conducted using the `goana` function. -Here were apply `goana` to the output of the TREAT analysis comparing `B.lactating` to `B.pregant`. -The top most significantly enriched GO terms can be viewed with `topGO`. - -```{r goana} -go <- goana(tr, species="Mm") -topGO(go, n=15) -``` -The `goana` function automatically extracts DE genes from the `tr` object, and conducts overlap tests for the up- and down-regulated DE genes separately. -By default, an FDR cutoff of 5% is used when extracting DE genes, but this can be varied. -The row names of the output are the universal identifiers of the GO terms and the `Term` column gives the human-readable names of the terms. -The `Ont` column shows the ontology domain that each GO term belongs to. -The three domains are: biological process (BP), cellular component (CC) and molecular function (MF). -The `N` column represents the total number of genes annotated with each GO term. -The `Up` and `Down` columns indicate the number of genes within the GO term that are significantly up- and down-regulated in this differential expression comparison, respectively. -The `P.Up` and `P.Down` columns contain the p-values for over-representation of the GO term in the up- and down-regulated genes, respectively. -Note that the p-values are not adjusted for multiple testing---we would usually ignore GO terms with p-values greater than about $10^{-5}$. - -By default the output table from `topGO` is sorted by the minimum of `P.Up` and `P.Down`. -Other options are available. -For example, `topGO(go, sort="up")` lists the top GO terms that are over-represented in the up-regulated genes. -The domain of the enriched GO terms can also be specified by users. -For example, `topGO(go, ontology="BP")` restricts to the top GO terms belonging to the biological process domain -while `topGO(go, ontology="MF")` restricts to molecular function terms. - -The `goana` function uses the NCBI RefSeq annotation and requires the use of Entrez Gene Identifiers. - -KEGG pathway analysis {#kegg-pathway-analysis .unnumbered} ---------------------- - - -Another popular annotation database is the Kyoto Encyclopedia of Genes and Genomes (KEGG). -Much smaller than GO, this is a curated database of molecular pathways and disease signatures. -A KEGG analysis can be done exactly as for GO, but using the `kegga` function: - -```{r kegga} -keg <- kegga(tr, species="Mm") -topKEGG(keg, n=15, truncate=34) -``` - -The output from `topKEGG` is the same as from `topGO` except that row names become KEGG pathway IDs, `Term` becomes `Pathway` and there is no `Ont` column. -Both the GO and KEGG analyses show that the cell cycle pathway is strongly down-regulated upon lactation in mammary stem cells. - -By default, the `kegga` function automatically reads the latest KEGG annotation from the Internet each time it is run. -The KEGG database uses Entrez Gene Ids, and the `kegga` function assumes these are available as the row names of `tr`. - -FRY gene set tests {#fry-gene-set-tests .unnumbered} ------------------- - - -The GO and KEGG analyses shown above are relatively simple analyses that rely on a list of DE genes. -The list of DE genes is overlapped with the various GO and KEGG annotation terms. -The results will depend on the significance threshold that is used to assess differential expression. - -If the aim is to test for particular gene expression signatures or particular pathways, a more nuanced approach is to conduct a `roast` or `fry` gene set test [@wu2010roast]. -These functions test whether a set of genes is DE, assessing the whole set of genes as a whole. -Gene set tests consider all the genes in the specified set and do not depend on any pre-emptive significance cutoff. -The set of genes can be chosen to be representative of any pathway or phenotype of interest. - -`roast` gives p-values using random rotations of the residual space. -In the *edgeR* context, `fry` is generally recommended over `roast`. -`fry` gives an accurate analytic approximation to the results that `roast` would give, with default settings, if an extremely large number of rotations was used. - -Here, suppose we are interested in two GO terms related to cytokinesis. -Each GO term is used to define a set of genes annotated with that term. -The names of these terms are shown below: - -```{r select, message=FALSE} -library(GO.db) -cyt.go <- c("GO:0032465", "GO:0000281") -term <- select(GO.db, keys=cyt.go, columns="TERM") -term -``` - -The first step is to extract the genes associated with each GO term from the GO database. -This produces a list of two components, one for each GO term. -Each component is a vector of Entrez Gene IDs for that GO term: - -```{r GO2ALLEGS} -Rkeys(org.Mm.egGO2ALLEGS) <- cyt.go -cyt.go.genes <- as.list(org.Mm.egGO2ALLEGS) -``` - -Suppose the comparison of interest is between the virgin and lactating groups in the basal population. -We can use `fry` to test whether the cytokinesis GO terms are DE for this comparison: - -```{r fry} -B.VvsL <- makeContrasts(B.virgin-B.lactating, levels=design) -fry(y, index=cyt.go.genes, design=design, contrast=B.VvsL) -``` - -Each row of the output corresponds to a gene set. -The `NGenes` column provides the number of genes in each set. -The `Direction` column indicates the net direction of change. -The `PValue` column gives the two-sided p-value for testing whether the set is DE as a whole, either up or down. -The `PValue.Mixed` column gives a p-value for testing whether genes in the set tend to be DE, without regard to direction. -The `PValue` column is appropriate when genes in the set are expected to be co-regulated, all or most changing expression in the same direction. -The `PValue.Mixed` column is appropriate when genes in the set are not necessarily co-regulated or may be regulated in different directions for the contrast in question. -FDRs are calculated from the corresponding p-values across all sets. - -The results of a gene set test can be viewed in a barcode plot produced by the `barcodeplot` function. -Suppose visualization is performed for the gene set defined by the GO term GO:0032465: - -```{r barcode, fig.width=8, fig.height=5, fig.cap="Barcode plot showing enrichment of the GO term GO:0032465 in the basal virgin group compared to the basal lactating group. X-axis shows logFC for B.virgin vs B.lactating. Black bars represent genes annotated with the GO term. The worm shows relative enrichment."} -res <- glmQLFTest(fit, contrast=B.VvsL) -index <- rownames(fit) %in% cyt.go.genes[[1]] -barcodeplot(res$table$logFC, index=index, labels=c("B.lactating","B.virgin"), - main=cyt.go[1]) -``` - -In the plot, all genes are ranked from left to right by increasing log-fold change for the contrast and the genes within the gene set are represented by vertical bars, forming the barcode-like pattern. -The curve (or *worm*) above the barcode shows the relative local enrichment of the bars in each part of the plot. -The dotted horizontal line indicates neutral enrichment; the worm above the dotted line shows enrichment while the worm below the dotted line shows depletion. -In this particular barcode plot the worm shows enrichment on the right for positive logFCs, and depletion on the left for negative logFCs. -The conclusion is that genes associated with this GO term tend to be up-regulated in the basal cells of virgin mice compared to lactating mice, confirming the result of the `fry` test above. - -Camera gene set enrichment analysis {#camera-gene-set-enrichment-analysis .unnumbered} ------------------------------------ - -Finally we demonstrate a gene set enrichment style analysis using the Molecular Signatures Database (MSigDB) [@subramanian2005gsea]. -We will use the C2 collection of the MSigDB, which is a collection of nearly 5000 curated gene sets, each representing the molecular signature of a particular biological process or phenotype. -The MSigDB itself is purely human, but the Walter and Eliza Hall Institute (WEHI) maintains a mouse version of the database. -We load the mouse version of the C2 collection from the WEHI website: -```{r camera} -load(url("http://bioinf.wehi.edu.au/software/MSigDB/mouse_c2_v5p1.rdata")) -``` -This will load `Mm.c2`, which is a list of gene sets, each a vector of Entrez Ids. -This can be converted to a list of index numbers: - -```{r id2indices.Mm.c2} -idx <- ids2indices(Mm.c2,id=rownames(y)) -``` - -First we compare virgin stem cells to virgin luminal cells: - -```{r cam.BVvsLV} -BvsL.v <- makeContrasts(B.virgin - L.virgin, levels=design) -cam <- camera(y, idx, design, contrast=BvsL.v, inter.gene.cor=0.01) -options(digits=2) -head(cam,14) -``` -With a large gene set collection, setting `inter.gene.cor = 0.01` gives a good compromise between biological interpretability and FDR control. -As expected, the mammary stem cell and mammary luminal cell signatures from Lim et al [@lim2010transcriptome] are top-ranked, and in the expected directions. - -We can visualize the top signature, combining the up and down mammary stem cell signatures to make a bi-directional signature set: - -```{r barcode2, fig.width=8, fig.height=6.4, fig.cap="Barcode plot showing strong enrichment of mammary stem cell signature in the stem cell vs luminal cell comparison. Red bars show up signature genes, blue bars show down genes. The worms show relative enrichment."} -res <- glmQLFTest(fit, contrast=BvsL.v) -barcodeplot(res$table$logFC, - index=idx[["LIM_MAMMARY_STEM_CELL_UP"]], - index2=idx[["LIM_MAMMARY_STEM_CELL_DN"]], - labels=c("L.virgin","B.virgin"), - main="LIM_MAMMARY_STEM_CELL", - alpha=1) -``` - - -Packages used {#packages-used .unnumbered} -============= - - -This workflow depends on various packages from version 3.3 of the Bioconductor project, running on R version 3.3.0 or higher. -The complete list of the packages used for this workflow are shown below: - -```{r info} -sessionInfo() -``` - - -Read alignment and quantification {#read-alignment-and-quantification .unnumbered} -================================= - -Download raw sequence files from the SRA {#download-raw-sequence-files-from-the-sra .unnumbered} ----------------------------------------- - -We now revisit the question of recreating the matrix of read counts from the raw sequence reads. -Unlike the above workflow, which works for any version of R, read alignment requires Unix or Mac OS and, in practice, a high performance Unix server is recommended. -Read alignment and read counting require only one Bioconductor package, *Rsubread*. -However the `fastq-dump` utility from the SRA Toolkit is also required to convert from SRA to FASTQ format. -This can be downloaded from the NCBI website () and installed on any Unix system. - -The first task is to download the raw sequence files, which are stored in SRA format on the SRA repository. -The SRA files need to be unpacked into FASTQ format using the `fastq-dump` utility. -The following R code makes a system call to `fastq-dump` to download each SRA file and convert it to FASTQ format: - - for (sra in targets$SRA) system(paste("fastq-dump", sra)) - -The `fastq-dump` utility automatically downloads the specified SRA data set from the internet. -The above code will produce 12 FASTQ files, in the current working directory, with file names given by the following vector: - - all.fastq <- paste0(targets$SRA, ".fastq") - - -Accuracy of base-calling {#accuracy-of-base-calling .unnumbered} ------------------------- - - -Sequencers typically store base-calling quality scores for each read in the FASTQ files. -*Rsubread*'s `qualityScores` function can be used to extract these scores from any particular file: - - QS <- qualityScores("SRR1552444.fastq") - -The boxplot function provides a compact way to view the quality scores by position across all reads: - - boxplot(QS, ylab="Quality score", xlab="Base position", - main="SRR1552444.fastq", cex=0.25, col="orange") - - -
-Boxplots of quality scores by base position for the first FASTQ file. -

-Figure 12: Boxplots of quality scores by base position for the first FASTQ file.
-

-
- - -The vertical axis shows the Phred quality score, equal to $-10\log_{10}(p)$ where $p$ is the probability of an erroneous call. -The maximum possible value is 40, and all values above 10 correspond to extremely small error probabilities. -The horizontal axis shows position within a read. -The file contains 100bp single-end reads, so the scale is from 1 to 100. -The plot displays a compact boxplot at each base position. -As is very commonly observed, the quality scores are best in the middle of the reads and decrease slightly towards the start and end of the reads. -However the quality remains generally good even near the ends of the reads: the scores would need to be very much lower than this before they would cause problems for the alignment. -Similar plots can be made for each of the FASTQ files. - - -Build a genome index {#build-a-genome-index .unnumbered} --------------------- - -Before the sequence reads can be aligned, we need to build an index for the GRCm38/mm10 (Dec 2011) build of the mouse genome. -Most laboratories that use *Rsubread* regularly will already have an index file prepared, as this is a once-off operation for each genome release. -If you are using *Rsubread* for mouse for the first time, then the latest mouse genome build can be downloaded from the NCBI webpage -(). -The full URL is: +完整网址为: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz -(Note that this link is for patch 6 of the RefSeq mm10 build, which is valid at the time of writing in May 2018. -The link will change as new patches are released periodically.) -An index can then be built by: +(注意，这个链接是针对RefSeq mm10构建的补丁6的，它在2018年5月撰写本文时是有效的。 +随着新补丁的定期发布，链接也会发生变化。) +然后可以通过以下方法建立索引: library(Rsubread) buildindex(basename = "mm10", reference = "GCF_000001635.26_GRCm38.p6_genomic.fna.gz") -Aligning reads {#aligning-reads .unnumbered} +比对 reads {#aligning-reads .unnumbered} -------------- -The sequence reads can now be aligned to the mouse genome using the `align` function: +现在使用“align”函数将序列read与小鼠基因组对齐: all.bam <- sub(".fastq", ".bam", all.fastq) align(index="mm10", readfile1=all.fastq, input_format="FASTQ", output_file=all.bam) -This produces a set of BAM files containing the read alignments for each RNA library. -The mapping proportions can be summarized by the `propmapped` function: +这将生成一组BAM文件，其中包含每个RNA库的读取对齐。 +映射比例可以通过‘propmapping’函数统计: > propmapped(all.bam) Samples NumTotal NumMapped PropMapped @@ -1781,31 +842,31 @@ The mapping proportions can be summarized by the `propmapped` function: 11 SRR1552448.bam 31425424 27326500 0.870 12 SRR1552449.bam 31276061 27204156 0.870 -Ideally, the proportion of mapped reads should be above 80%. -By default, only reads with unique mapping locations are reported by *Rsubread* as being successfully mapped. -Restricting to uniquely mapped reads is recommended, as it avoids spurious signal from non-uniquely mapped reads derived from, e.g., repeat regions. +理想情况下，映射读取的比例应该在80%以上。 +默认情况下，只有具有唯一映射位置的读取才会被“Rsubread”包报告为成功映射。 +建议将读取限制为惟一映射，因为这样可以避免来自重复区域等非惟一映射读取的伪信号。 -Quantifying read counts for each gene {#quantifying-read-counts-for-each-gene .unnumbered} +量化每个基因的read计数 {#quantifying-read-counts-for-each-gene .unnumbered} ------------------------------------- -The read counts for each gene can be quantified using the `featureCounts` function in *Rsubread*. -Conveniently, the *Rsubread* package includes inbuilt NCBI RefSeq annotation of the mouse and human genomes. -`featureCounts` generates a matrix of read counts for each gene in each sample: +可以使用“Rsubread”包的“featureCounts”函数量化每个基因的read计数。 +“Rsubread”包内含有内建的NCBI RefSeq注释小鼠和人类基因组。 +“featuts”为每个样本中的每个基因生成一个读取计数矩阵: fc <- featureCounts(all.bam, annot.inbuilt="mm10") -The output is a simple list, containing the matrix of counts (`counts`), a data frame of gene characteristics (`annotation`), a vector of file names (`targets`) and summary mapping statistics (`stat`): +输出是一个简单的列表，包含计数矩阵(' counts ')、基因特征数据框(' annotation ')、文件名向量(' targets ')和汇总映射统计信息(' stat '): > names(fc) [1] "counts" "annotation" "targets" "stat" -The row names of `fc$counts` are the Entrez gene identifiers for each gene. -The column names are the output file names from `align`, which we simplify here for brevity: +“fc$counts”的行名是每个基因的Entrez基因标识符。 +列名是来自“align”的输出文件名，这里我们将其简化： > colnames(fc$counts) <- rownames(targets) -The first six rows of the counts matrix are shown below. +count矩阵的前六行如下所示。 > head(fc$counts) MCL1.DG MCL1.DH MCL1.DI MCL1.DJ MCL1.DK MCL1.DL MCL1.LA @@ -1823,56 +884,56 @@ The first six rows of the counts matrix are shown below. 20671 25 18 8 3 10 27395 468 488 332 312 344 -Finally, a `DGEList` object can be assembled by: +最后，一个DGEList对象可以通过以下方式组装: y <- DGEList(fc$counts, group=group) y$genes <- fc$annotation[, "Length", drop=FALSE] -Data and software availability {#data-and-software-availability .unnumbered} +使用的数据和软件 {#data-and-software-availability .unnumbered} ============================== -This html document was generated as a vignette of the RnaSeqGeneEdgeRQL package, which can be installed by running +这个html文档是作为RnaSeqGeneEdgeRQL包的概要文件生成的，可以在R中通过运行安装RnaSeqGeneEdgeRQL包实现。 install.packages("BiocManager") BiocManager::install("RnaSeqGeneEdgeRQL") -at the R prompt. -Except for the targets file `targets.txt`, all data analyzed in the workflow is read automatically from public websites as part of the code. -All software used is publicly available as part of Bioconductor (Release 3.3 or later), except for the `fastq-dump` utility, which can be downloaded from NCBI website as described in the text. -The article includes the complete code necessary to reproduce the analyses shown. -The original Bioconductor 3.3 version of this workflow, including a pdf article version and a corresponding Rnw file, is available from and [@chen2016quasiworkflow]. +除了目标文件“targets.txt”，所有在工作流程中分析的数据将作为代码的一部分自动从公共网站读取。 +除“fastq-dump”程序之外，使用的所有软件都是作为Bioconductor(3.3版或更高版本)的一部分公开提供的，如文中所述，该程序可以从NCBI网站下载。 +本文包含了重现所示分析所必需的完整代码。 + +该工作流的原始Bioconductor3.3版本，包括pdf文章版本和相应的Rnw文件，可以从和[@chen2016quasiworkflow]获得。 -Author contributions {#author-contributions .unnumbered} +作者贡献 {#author-contributions .unnumbered} ==================== -All authors developed and tested the code workflow. -All authors wrote the article. +所有作者都开发并测试了代码工作流。 +所有作者都参与了这篇文章的写作。 -Grant information {#grant-information .unnumbered} +授权信息 {#grant-information .unnumbered} ================= -This work was supported by the National Health and Medical Research Council +这项工作得到了国家卫生和医学研究理事会的支持 (Fellowship 1058892 and Program 1054618 to G.K.S, Independent Research Institutes Infrastructure Support to the Walter and Eliza Hall Institute) -and by a Victorian State Government Operational Infrastructure Support Grant. +并由维多利亚州政府运营基础设施支持赠款。 -Acknowledgments {#acknowledgments .unnumbered} +致谢 {#acknowledgments .unnumbered} =============== -The authors thank Wei Shi and Yang Liao for advice with *Rsubread* and Yifang Hu for creating the mouse version of the MSigDB. -We also wish to acknowledge the early developers of the quasi-likelihood testing approach. -Davis McCarthy authored the first versions of the quasi-likelihood functions in *edgeR* in February 2011 and Steve Lund and Dan Nettleton worked independently on *QuasiSeq* around the same time. -We thank Steve Lund and Dan Nettleton for a valuable and enjoyable collaboration that lead to Lund et al [@lund2012quasiseq]. +感谢Wei Shi 和Yang Liao对“Rsubread”包提出的建议，感谢Yifang Hu 创建MSigDB的小鼠版本。 +我们也希望感谢准似然测试方法的早期开发人员。 +2011年2月，Davis McCarthy在“edgeR”包中撰写了准似似函数的第一个版本，Steve Lund和Dan Nettleton大约在同一时间独立研究了“QuasiSeq”包。 +我们感谢Steve Lund和Dan Nettleton的宝贵而愉快的合作，这促成了Lund等人[@lund2012quasiseq]。 -How to cite {#how-to-cite .unnumbered} +如何引用 {#how-to-cite .unnumbered} ==================== -Please cite this workflow article as +请将此工作流文章引用为 Chen, Yunshun, A. T. L. Lun, and G. K. Smyth. 2016. "From Reads to Genes to Pathways: Differential Expression Analysis of RNA-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline." *F1000Research* 5:1438. -References {#references .unnumbered} +参考 {#references .unnumbered} ========== From 92c8d39975fb576a4e01ffce5729c06cdf90c3a2 Mon Sep 17 00:00:00 2001 From: restore1997 <47097886+restore1997@users.noreply.github.com> Date: Thu, 1 Aug 2019 23:02:35 +0800 Subject: [PATCH 3/3] Update 10RnaSeqGeneEdgeRQL_en.Rmd --- BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd index f2ae507..ec88142 100644 --- a/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd +++ b/BiocWorkflow_todo/10RnaSeqGeneEdgeRQL_en.Rmd @@ -930,7 +930,7 @@ count矩阵的前六行如下所示。 如何引用 {#how-to-cite .unnumbered} ==================== -请将此工作流文章引用为 +请如下引用该篇文章 Chen, Yunshun, A. T. L. Lun, and G. K. Smyth. 2016. "From Reads to Genes to Pathways: Differential Expression Analysis of RNA-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline." *F1000Research* 5:1438.