\VignetteIndexEntry{BioC2014: RNA-Seq差分基因表达工作流}%\VignettePackage{DESeq2} %\VignetteEngine{knitr::knitr} %编译此文档% library('knitr');rm(列表= ls ());knit('RNA-Seq-Analysis-Lab.Rnw') \documentclass[12pt]{article} \newcommand{\deseqtwo}{\Biocpkg{DESeq2}} \newcommand{\lowtilde}{\raise. txt ')17ex\hbox{$\scriptstyle\mathtt{\sim}$}} \usepackage{xcolor} \usepackage[framemethod=default]{mdframed} \usepackage{amsmath,amsthm} \定理样式{定义}\ newmdtheormenv [linecolor=红色，backgroundcolor=粉色!]50]{exer}{Exercise} \widowpenalty10000 \clubpenalty10000 < >= library("knitr") opts_chunk$set(tidy=FALSE,dev="png"，图。show="hide"， fig.width=4,fig.height=4.5, message=FALSE) @ < >= BiocStyle::latex() @ \作者{Michael Love$^{1*}$， Simon Anders$^{2}$， Wolfgang Huber$^{2}$ \\[1em] \small{$^{1}$ Dana Farber癌症研究所生物统计系和}\\ \small{哈佛大学公共卫生学院，美国波士顿;}\\ \\ small{$^{2}$欧洲分子生物学实验室(EMBL)，德国海德堡}\\ \small{\texttt{$^*$michaelisaiahlove (at) gmail.com}}} \title{BioC2014:RNA-Seq差异基因表达工作流\begin{document} \maketitle \begin{abstract}本实验室将带您完成端到端RNA-Seq差异表达工作流。我们将从FASTQ文件开始，对准参考基因组，通过计数测序片段准备基因表达值作为计数表，进行差异基因表达分析，并可视化地探索结果。本实验涵盖了基本的介绍。要更深入地解释高级细节，我们建议您继续阅读\deseqtwo{}， \emph{计数数据的差分分析——\emph{DESeq2}包}的小插图。对于外显子水平差异表达的处理，我们参考\Biocpkg{DEXSeq}包的小插曲，\emph{分析差异外显子使用的RNA-seq数据\emph{DEXSeq}包}。%基本基因集富集\end{abstract} < >= options(数字=3，宽度=80)@ \tableofcontents \section{输入数据}\label{Haglund}作为本实验室的示例数据，我们将使用Felix Haglund等人的文章，\emph{甲状旁腺腺瘤中功能性雌激素受体的证据}，J Clin Endocrin Metab, 2012年9月\脚注{\url{http://www.ncbi.nlm.nih.gov/pubmed/23024189}}。本实验的目的是探讨雌激素受体在甲状旁腺肿瘤中的作用。研究人员从4名患者身上提取了甲状旁腺瘤细胞的原代培养。这些原代培养用二酰丙腈(DPN)、雌激素受体激动剂或4-羟他莫西芬(OHT)处理。分别于24小时和48小时从处理和对照培养物中提取RNA。实验的阻塞设计允许对治疗效果进行统计分析，同时控制患者之间的差异。本实验的部分数据在Bioconductor数据包\Biocexptpkg{parathyroidSE}中提供。\分段{准备计数矩阵}作为输入，\deseqtwo{}包期望获得的计数数据，e.\，g。，来自RNA-Seq或其他高通量测序实验，以整数值矩阵的形式。矩阵的$i$-第一行和$j$-第列的值表示样本$j$中有多少个读取被映射到基因$i$。类似地，对于其他类型的分析，矩阵的行可能对应e.\，g。\ to binding regions (with ChIP-Seq) or peptide sequences (with quantitative mass spectrometry). The count values must be raw counts of sequencing reads. This is important for \deseqtwo{}'s statistical model to hold, as only the actual counts allow assessing the measurement precision correctly. Hence, please do not supply other quantities, such as (rounded) normalized counts, or counts of covered base pairs -- this will only lead to nonsensical results. \subsection{Aligning reads to a reference} The computational analysis of an RNA-Seq experiment begins earlier: what we get from the sequencing machine is a set of FASTQ files that contain the nucleotide sequence of each read and a quality score at each position. These reads must first be aligned to a reference genome or transcriptome. It is important to know if the sequencing experiment was single-end or paired-end, as the alignment software will require the user to specify both FASTQ files for a paired-end experiment. The output of this alignment step is commonly stored in a file format called \href{http://samtools.github.io/hts-specs}{\texttt{BAM}}. A number of software programs exist to align reads to the reference genome, and the development is too rapid for this document to provide a current list. We recommend consulting benchmarking papers that discuss the advantages and disadvantages of each software, which include accuracy, ability to align reads over splice junctions, speed, memory footprint, and many other features. Here we use the TopHat2 spliced alignment software\footnote{\url{http://tophat.cbcb.umd.edu/}} \cite{Kim2013TopHat2} in combination with the Bowtie index available at the Illumina iGenomes page\footnote{\url{http://tophat.cbcb.umd.edu/igenomes.html}}. For full details on this software and on the iGenomes, users should follow the links to the manual and information provided in the links in the footnotes. For example, the paired-end RNA-Seq reads for the \Biocexptpkg{parathyroidSE} package were aligned using TopHat2 with 8 threads, with the call: \begin{verbatim} tophat2 -o file_tophat_out -p 8 genome file_1.fastq file_2.fastq \end{verbatim} The BAM files for a number of sequencing runs can then be used to generate count matrices, as described in the following section. \subsection{Example BAM files} Besides the main count table, which we will use later, the \Biocexptpkg{parathyroidSE} package also contains a small subset of the raw data from the Haglund et al.\ experiment, namely three BAM file each with a subset of the reads from three of the samples. We will use these files to demonstrate how a count table can be constructed from BAM files. Afterwards, we will load the full count table corresponding to all samples and all data, which is already provided in the same package, and will continue the analysis with that full table. We load the data package with the example data < >= library("parathyroidSE") @ R函数\Rfunction{系统。File}可用于查找软件包中的文件已安装在计算机上的哪个位置。这里我们请求\texttt{extdata}目录的完整路径，该目录是\Biocexptpkg{parathyroidSE}包的一部分:< >= extDataDir <- system. bat。在这个目录中，我们找到了三个BAM文件(以及一些其他文件):< > =列表。通常，我们有一个表，其中包含样本的实验元数据。对于这三个文件，它如下所示:\begin{tabular}{llll} sampleName & fileName & treatment & time \\ \hline Ctrl\_24h\_1 & SRR479052。bam & Control & 24h \\ Ctrl\_48h\_1 & SRR479053。bam & Control & 48h \\ DPN\_24h\_1 & SRR479054。为了避免错误，将这样一个示例表显式存储在一个文本文件中并加载它是很有帮助的。使用%电子表格应用程序(如Microsoft Excel)在磁盘上构造一个包含上述内容的表格，并以CSV(逗号分隔值)格式保存表格%(或简单地使用纯文本%编辑器)。% \end{exer}加载:< >= sampleTable <- read.csv("/path/to/your/sampleTable.csv"， header=TRUE) @ < >= sampleTable <- data.frame(sampleName = c("Ctrl_24h_1"， "Ctrl_48h_1"， "DPN_24h_1")， fileName = c("SRR479052. "bam”、“SRR479053。bam”、“SRR479054。bam")， treatment = c("Control"， "Control"， "DPN")， time = c("24h"， "48h"， "24h")) @这是样例表应该看起来的样子< 使用表中的fileName列，我们构造要对其执行计数操作的文件的完整路径:< >= bamFiles <- file。我们可以查看其中一个BAM文件，以查看序列(染色体)的命名风格。这里我们使用了\Biocpkg{Rsamtools}包中的\Rfunction{BamFile}函数。< >= library("Rsamtools") seqinfo(BamFile(bamFiles[1])) @我们要确保这些序列名与我们将在下一节中获得的基因模型的序列名是相同的风格。要计算每个基因有多少读映射，我们需要转录注释。从Ensembl下载带有人类基因注释的当前GTF文件。(如果网络太慢，可以使用这个文件的截短版本，名为\texttt{Homo\_sapiens. grch37.75 . set.gtf.gz}，我们已经把它放在课程服务器上了。)从这个文件中，\Biocpkg{GenomicFeatures}中的函数\Rfunction{makeTranscriptDbFromGFF}构造了一个包含所有注释转录本的数据库。< >= library("GenomicFeatures") hse <- makeTranscriptDbFromGFF("/path/to/your/genemodel_file. "exonsByGene <- exonsBy(hse, by="gene") @ < >= library("GenomicFeatures") hse <- makeTranscriptDbFromGFF(" homo_sapians . grch37.75 .subset.gtf.gz"， format="gtf") exonsByGene <- exonsBy(hse, by="gene") @在最后一步中，我们使用\Rfunction{exonsBy}函数将转录组数据库变成所有基因的列表形状，< >= exonsByGene @ \begin{exer}注意\Rfunction{makeTranscriptDbFromGFF}发出的警告。我们能安全地忽略它们吗?\end{exer} \begin{exer}在\Robject{exonsByGene}中，检查第一个基因的外显子的基因组间隔。注意它们不是不相交的(它们是重叠的)。为什么?这会影响下一步的结果吗?在这些准备之后，实际的计数就很容易了。\Biocpkg{GenomicAlignments}包中的函数\Rfunction{summarizeOverlaps}将完成此操作。< >= library("GenomicAlignments") se <- summarizeOverlaps(exonsByGene, BamFileList(bamFiles)， mode="Union"， singleEnd=FALSE，忽略。我们使用计数模式\Robject{"Union"}，这表明那些恰好重叠一个特征的任何部分的读取被计数。有关不同计数模式的更多信息，请参阅\Rfunction{summarizeOverlaps}的帮助页。由于这个实验产生了对端读取，我们指定\Robject{singleEnd = FALSE}。由于协议不是特定于字符串的，我们指定\Robject{ignore。strand = TRUE}。\Robject{fragments = TRUE}表示我们还想对未映射的对进行读取计数。最后一个参数只适用于配对实验。有关如何从BAM文件读取的详细信息可以使用\Rfunction{BamFileList}函数指定。例如，为了控制内存，我们可以指定一次读取$2 000 000$的批次:< >= BamFileList(bamFiles, yieldSize = 2000000) @请记住，我们只使用了原始实验中读取的一小部分:3个样本和100个基因。尽管如此，我们仍然可以通过查看\Robject{assay}插槽中的计数、\Robject{colData}插槽中样本的表型数据(在这种情况下是空的\Rclass{DataFrame})以及\Robject{rowData}插槽中的基因数据来研究结果\Rclass{summarizeexperiment}。图~\ref{Figure /beginner_plotSE}解释了\r类{summarizeexperiment}类的基本结构。< >= se head(化验(se)) colsum(化验(se)) colData(se) rowData(se) @注意\Robject{rowData}槽是一个\Rclass{GRangesList}，它包含了每个基因外显子的所有信息，即计数表的每一行。\Robject{colData}插槽，到目前为止是空的，应该包含所有的元数据。因此，我们将示例表赋值给它:< >= colData(se) <- DataFrame(sampleTable) @我们可以使用\Robject{\$}操作符从\Robject{colData}提取列，我们可以省略\Rfunction{colData}以避免额外的击键。< 我们还可以使用sampleName表来命名数据矩阵的列:< >= colnames(se) <- sampleTable$sampleName head(化验(se)) @这个\Rclass{summarizeexperiment}对象\Robject{se}是我们开始分析所需要的。< > = par (mar = c(0, 0, 0, 0))的阴谋(1,1,xlim = c (0100), ylim = c(0100),电池=“n”,类型=“n”,xlab = " ", ylab = " ", xaxt =“n”,yaxt =“n”)多边形(c(45、80、80、45),c(10、10、70、70),坳= rgb(1 0 0 5),边界= NA)多边形(c(45、80、80、45),c(68, 68, 70, 70),坳= rgb(1 0 0 5),边界= NA)文本(62.5,40岁的“试验(s)”)文本(如62.5,30。”计数”)多边形(c(20、40、40岁,20),c(10、10、70、70),坳= rgb(0, 0, 1, 0。5)边境= NA)多边形(c(20、40、40岁,20),c(68, 68, 70, 70),坳= rgb(0, 0, 1, 0。5)边境= NA)文本(30 40构成了rowData””)多边形(c(45、80、80、45),c(75, 75, 90, 90),坳= rgb(5 0。5。5)边境= NA)多边形(c(45岁,47岁,47岁,45岁),c(75, 75, 90, 90),坳= rgb(5 0。5。5)边境= NA)文本(62.5,82.5,“colData”)@ \ incfig图/ beginner_plotSE}{{。一个\Rclass{summarizeexperimental}对象的组成部分。\Robject{assay(s)}(红色块)包含汇总值的矩阵(或矩阵)，\Robject{rowData}(蓝色块)包含关于基因组范围的信息，\Robject{colData}(紫色块)包含关于样本或实验的信息。每个块中突出显示的行表示第一行(注意\Robject{colData}的第一行与\Robject{assay}的第一列对齐。{\分段{DESeqDataSet，列元数据和设计公式}Bioconductor软件包通常有一个特殊的数据对象类，其中包含特殊的插槽和需求。\deseqtwo{}中的数据对象类是\Rclass{DESeqDataSet}，它构建在\Rclass{summarizeexperiment}类之上。一个主要的区别是\Robject{assay}插槽是使用\Rfunction{count}访问器访问的，并且这个矩阵中的值必须是非负整数。第二个区别是\ r类{DESeqDataSet}有一个相关的“设计公式”。设计是在分析开始时指定的，因为它将通知许多\deseqtwo{}函数如何处理分析中的样本(一个例外是大小因子估计，即\，e。，不同库大小的调整，不依赖于设计公式)。设计公式告诉列元数据表(\Robject{colData})中的哪些变量指定实验设计，以及在分析中应该如何使用这些因素。 The simplest design formula for differential expression would be \Rfunction{\lowtilde{} condition}, where \Robject{condition} is a column in \Robject{colData(dds)} which specifies which of two (or more groups) the samples belong to. For the parathyroid experiment, we will specify \Rfunction{\lowtilde{} patient + treatment}, which means that we want to test for the effect of treatment (the last factor), controlling for the effect of patient (the first factor). You can use R's formula notation to express any experimental design that can be described within an ANOVA-like framework. Note that \deseqtwo{} uses the same formula notation as, for instance, the \Rfunction{lm} function of base R. If the question of interest is whether a fold change due to treatment is different across groups, for example across patients, ``interaction terms'' can be included using models such as \Rfunction{\lowtilde{} patient + treatment + patient:treatment}. More complex designs such as these are covered in the other \deseqtwo{} vignette. We now use R's \Rfunction{data} command to load a prepared \Rclass{SummarizedExperiment} that was generated from the publicly available sequencing data files associated with the Haglund et al.~paper, described on page~\pageref{Haglund}. The steps we used to produce this object were equivalent to those you worked through in Section~\ref{sec:count}, except that we used the complete set of samples and all reads. < >= data("parathyroidGenesSE") se <- parathyroidGenesSE @关于我们上面展示的工作流的一个额外好处是，我们使用的基因模型的信息不需要额外的努力就包括在内。甲状旁腺数据的\Robject{rowData}是通过\Rfunction{makeTranscriptDbFromBiomart}函数获得的，我们可以找到基因模型版本的所有元数据信息。\Rfunction{str} R函数用于紧凑地显示列表中数据的结构。< 假设我们已经使用前一节中描述的方法之一构造了一个\Rclass{summarizeexperimental}，我们现在需要确保对象包含关于样本的所有必要信息，即，一个在计数表的列上存储在\Robject{colData}插槽的元数据表:< >= colData(se)[1:5,1:4] @在这里，我们看到这个对象已经包含了一个信息性的\Robject{colData}插槽——因为我们已经为您准备好了它，正如\Biocexptpkg{parathyroidSE}小插图中所描述的那样。然而，当你使用你自己的数据时，你必须在这个阶段为实验添加相关的样本/表型信息。我们强烈建议将此信息保存在逗号分隔值(CSV)或制表符分隔值(TSV)文件中，该文件可以从Excel电子表格导出，并将其分配到\Robject{colData}插槽，如前一节所示。确保列数据表中的行顺序与\Robject{assay}数据槽中的列顺序匹配。在上一节中构建\Rclass{se}对象时，我们如何确保这一点?一旦我们有了完全注释的\Rclass{summarizeexperiment}对象，我们就可以从中构造一个\Rclass{DESeqDataSet}对象，它将形成实际\deseqtwo{}包的起始点，这将在下面几节中描述。在这里，我们使用从\Biocexptpkg{parathyroidSE}包中获得的\Rclass{summarizeexperiment}对象，并通过指定适当的设计公式来扩充它。< > =图书馆(DESeq2) ddsFull <——DESeqDataSet病人(se设计= ~ +治疗)@请注意,有两个可选功能,\ Rfunction {DESeqDataSetFromMatrix}和\ Rfunction {DESeqDataSetFromHTSeq},这允许您开始在你的数据而不是形式的\ Rclass {SummarizedExperiment}对象,但作为一个简单的矩阵计算值或输出文件从\ emph {htseq-count}的脚本\ emph {HTSeq} Python包。\分段{崩溃技术复制}有许多样品在多次运行中进行了测序。例如，样本\Robject{SRS308873}被测序了两次。为此，我们列出\Rfunction{colData}的各个列。(\Rfunction{as.data.frame}的使用迫使R向我们显示完整的列表，而不是像以前那样只显示开始和结束。)< >= as.data.frame(colData(ddsFull)[，c("sample"，"patient"，"treatment"，"time")] @我们建议首先将技术复制(即来自相同样本的库)添加在一起，这样每个样本都有一列。我们为此实现了一个方便的函数，它可以获取一个对象(\Rclass{summarizeexperiment}或\Rclass{DESeqDataSet})和一个分组因子(在本例中为样本名)，并返回每个唯一样本的求和计数对象。这还将重命名对象的列，使它们与分组因子中使用的惟一名称相匹配。可选地，我们可以提供第三个参数\Robject{run}，它可用于将折叠后创建新对象的运行的名称粘贴在一起。注意\Robject{dds\$variable}等价于\Robject{colData(dds)\$variable}。< >= ddscollapse <- collapsereplcopies (ddsFull, groupby = ddsFull$sample, run = ddsFull$run) head(as.data.frame(colData(ddscollapse)[，c("sample"，" runscollapse ")])， 12) @我们可以确认新对象的计数等于具有相同分组因子值的列的计数之和:< >= original <- rowsum (counts(ddsFull)[， ddsFull$sample == "SRS308873"]) all(original == counts(dds) [，"SRS308873"]) @ \section{运行DESeq2管道}\label{sec: rdiseq}在这里，我们将分析样本的一个子集，即48小时后的样本，考虑到控制、DPN或OHT处理，考虑到多因素设计。\分段{为感兴趣的分析准备数据对象}首先，我们从完整的数据集中子集相关列:< >= dds <- ddscollapse [， ddscollapse $time == "48h"] @有时有必要降低因子的级别，以防设计中一个或多个因子级别的所有样本都已被删除。如果在设计公式中包含时间，则可以使用以下代码来处理这一列中下降的级别。< >= dds$time <- droplevels(dds$time) @确保\Robject{Control}是治疗因子中的\emph{first}级别是很方便的，这样默认的log2倍变化被计算为治疗优于控制，而不是相反。函数\Rfunction{relevel}实现了以下功能:< >= dds$treatment <- relevel(dds$treatment， "Control") @一个快速检查我们现在是否有正确的样本:< 为了加快下面的一些注释步骤，删除对所有样本计数为零的基因是有意义的。这可以通过简单地索引\Robject{dds}对象来实现:< >= dds <- dds[rowsum (counts(dds)) > 0，] @ \分段{运行管道}最后，我们准备运行差分表达式管道。让我们回顾一下我们指定的设计:< \deseqtwo{}分析建模计算患者和治疗效果现在可以通过单个调用函数\Rfunction{DESeq}来运行:< >= dds <- DESeq(dds) @这个函数将为它执行的各个步骤打印一条消息。这些在\Rfunction{DESeq}的手册页中有更详细的描述，可以通过输入\Robject{?DESeq}来访问。简单地说，这些是:大小因素的估计(控制测序实验文库大小的差异)，每个基因离散度的估计，以及拟合广义线性模型。返回一个包含所有拟合信息的\Rclass{DESeqDataSet}，下面一节将描述如何从该对象中提取感兴趣的结果表。不带任何参数调用Rfunction{results}将提取设计公式中最后一个变量的估计log2倍变化和$p$值。如果这个变量有超过2个级别——就像本分析中的情况一样——\Rfunction{results}将提取结果表，以便比较最后一个级别与第一个级别。下一节描述如何提取其他比较。< >= res <- results(dds) res @由于\Robject{res}是一个\Robject{DataFrame}对象，它携带关于列含义的元数据:< >= mcols(res, use.names=TRUE) @第一列，\Robject{baseMean}，是所有样本归一化计数值的平均值，除以大小因子。其余四列表示特定的\emph{对比}，即因素变量\emph{处理}的\emph{DPN}水平与\emph{Control}水平的比较。请参阅\Rfunction{results}的帮助页面(输入\Rfunction{?results})，以了解如何获得其他对比。\Robject{log2FoldChange}列是效应大小估计值。它告诉我们，与对照组相比，由于DPN治疗，基因表达似乎发生了多少变化。这个值是以2为基数的对数刻度报告的:例如，log2倍的变化为1.5意味着基因的表达增加了2^{1.5}\约2.82美元的乘法因子。当然，这个估计有一个与之相关的不确定性，可以在\Robject{lfcSE}列中找到，log2倍变化估计的标准误差估计。我们还可以将特定效应量估计的不确定性表示为统计检验的结果。微分表达式检验的目的是检验数据是否提供了足够的证据来得出这个值确实不同于零的结论。\deseqtwo{}对每个基因执行一个\emph{假设检验}，看看是否有足够的证据来推翻\emph{原假设}，即处理对基因没有影响，并且处理与对照组之间观察到的差异仅仅是由实验变异性引起的(即，这种类型的变化，你也可以期待在同一处理组的不同样本之间)。 As usual in statistics, the result of this test is reported as a \emph{$p$ value}, and it is found in the column \Robject{pvalue}. (Remember that a $p$ value indicates the probability that a fold change as strong as the observed one, or even stronger, would be seen under the situation described by the null hypothesis.) We note that a subset of the $p$ values in \Robject{res} are \Robject{NA} (``not available''). This is \Rfunction{DESeq}'s way of reporting that all counts for this gene were zero, and hence not test was applied. In addition, $p$ values can be assigned \Robject{NA} if the gene was excluded from analysis because it contained an extreme count outlier. For more information, see the outlier detection section of the advanced vignette. We can examine the counts and normalized counts for the gene with the smallest p value: < >= idx <- that .min(res$pvalue) counts(dds)[idx，] counts(dds, normalized=TRUE)[idx，] @ \begin{exer}使用\Rfunction{plot}检查处理中p值最小的基因的归一化计数，即\Robject{that .min(res\$pvalue)}。尝试在处理列上使用\Rfunction{as.integer}来显示点而不是箱线图。尝试在patient列上使用\Rfunction{as.integer}，并更改点的绘图字符(pch)或颜色(col)。尝试\Robject{log="y"}，以便更清楚地看到差异。现在将其转换为一个函数，该函数可以将任何行下标作为参数。(绘制计数图是一种有用的诊断方法，因此，辅助函数\Rfunction{plotCounts}已被添加到2014年10月的下一个\deseqtwo{}版本中。)\end{exer} \分段{其他比较}一般来说，可以使用\Robject{contrast}参数对\Rfunction{results}提取变量的任意两个级别的比较结果。用户应该指定三个值:变量名、分子中的关卡名和分母中的关卡名。首先我们保存之前的结果表:< >= resOHT <- res @这里我们提取了DPN $/$ Control的折叠变化的log2的结果。< >= res <- results(dds, contrast = c("treatment"， "DPN"， "Control")) res @如果需要交互项的结果，应使用\Rfunction{results}的\Robject{name}参数。请参阅更高级的小插图了解更多细节。\分段{添加基因名}我们的结果表只使用Ensembl基因id，但基因名可能更有信息。Bioconductor的注释包有助于将各种ID方案相互映射。加载注释包\Rpackage{org.Hs.eg.db}: < 这是\emph{Homo sapiens} ('' Hs ")的生物体注释包('' org ")，组织为\Rpackage{AnnotationDbi}包('' db")，使用Entrez Gene IDs ('' eg ")作为主键。要获得所有可用键类型的列表，使用< 使用\Rpackage{AnnotationDbi}包中的本地函数转换id目前有点麻烦，所以我们提供了以下方便函数(没有解释它是如何工作的):< >= convertIDs <- function(ids, fromKey, toKey, db, ifMultiple=c("putNA"， "useFirst")) {stopifnot(继承(db， "AnnotationDb")) ifMultiple <- match。arg(ifMultiple) suppressWarnings(selRes <- AnnotationDbi::select(db, keys=ids, keytype=fromKey, columns=c(fromKey,toKey)) if(ifMultiple == "putNA") {duplicatedIds <- selRes[duplicate (selRes[，1])， 1] selRes <- selRes[!selRes[，1] %in% duplicatedIds，]} return(selRes[match(ids, selRes[，1])， 2])} @这个函数接受一个id列表作为第一个参数，它们的键类型作为第二个参数。第三个参数是我们要转换的键类型，第四个是要使用的\Rclass{AnnotationDb}对象。最后，最后一个参数指定了当一个源ID映射到多个目标ID时该怎么做:函数应该返回NA还是仅仅返回多个ID中的第一个?要将\Robject{res}行名中的Ensembl id转换为基因符号并将其添加为新列，我们使用:< >= res$hgnc_symbol <- convertIDs(row.names(res)， "ENSEMBL"， "SYMBOL"， org.Hs.eg.db) res$entrezid <- convertIDs(row.names(res)， "ENSEMBL"， " entrezid "， org.Hs.eg.db) @现在结果有所需的外部基因id: < >= head(res, 4) @ \begin{exer}转到Ensembl网站，选择\emph{BioMart}选项卡，并通过手动输入Ensembl ID并找到HGNC名称和Entrez ID(指定过滤器，属性，然后单击结果)来重做我们的BioMart查询。从这个集市还能得到什么其他数据?你可以修改上面的代码块，向\Robject{res}对象中添加一列\Robject{chrom}，告诉我们每个基因位于哪条染色体上?高通量生物学的新手通常认为，将这些$p$值设置在一个较低的值，比如0.01，就像在其他设置中经常做的那样，是合适的——但事实并非如此。我们简单地解释为什么:有\Sexpr{sum(res$pvalue < 0.01, na。\Sexpr{sum(!is.na(res$pvalue))}基因的$p$值低于0.01，测试成功报告了$p$值:< >= sum(res$pvalue < 0.01, na。现在，暂时假设原假设对所有基因都成立，即没有基因受到DPN治疗的影响。然后，根据$p$值的定义，我们期望有1%的基因的$p$值低于0.01。这相当于\Sexpr{round(sum(!is.na(res$pvalue)) * 0.01)}基因。如果我们只考虑$p$值低于0.01的基因列表作为差异表达，那么这个列表应该包含最多$\Sexpr{round(sum(!is.na(res$pvalue)) * 0.01)} / \Sexpr{sum(res$pvalue < 0.01, na。rm=TRUE)} = \Sexpr{round(sum(!is.na(res$pvalue)) * 0.01 / sum(res$pvalue < 0.01, na。rm=TRUE) * 100,0)}$\%误报!\deseqtwo{}使用所谓的Benjamini-Hochberg (BH)调整;简而言之，该方法为每个基因计算一个调整后的$p$值，回答以下问题:如果一个人称所有的$p$值小于或等于该基因的$p$值阈值的基因为显著性基因，那么其中的假阳性(错误发现率)的比例是多少(在上述计算的意义上)?这些值被称为BH-adjusted $p$ values，在\Robject{results}对象的\Robject{padj}列中给出。因此，如果我们认为10%的假阳性是可以接受的，我们可以认为所有调整后的p值低于10 %=0.1的基因都是显著的。有多少这样的基因? < >= sum(res$padj < 0.1, na。我们将结果表子集划分为这些基因，然后按log2倍变化估计值进行排序，以获得下调最强的显著基因< >= resSig <-子集(res, res$padj < 0.1) head(resSig[order(resSig$log2FoldChange)，]， 4) @ with with最强上调< >= head(resSig[order(-resSig$log2FoldChange)，]， 4) @ \分段{诊断图}所谓的MA图为两组比较的实验提供了有用的概述:< >= plotMA(res, ylim = c(- 3,3)) @图(图\ \ref{figure/beginner_MA})用一个点表示每个基因。x轴是所有样本的平均表达式，y轴是处理和对照组之间的log2倍变化。调整后的$p$值低于阈值(此处为0.1，默认值)的基因用红色表示。图/ beginner_MA \ incfig{}{。{MA-plot}{MA-plot显示了从处理到标准化计数的平均值的log2倍变化，即按大小因子标准化的计数的平均值。\deseqtwo{}包包含了log2倍的先验变化，导致来自低计数和高可变计数的基因的适度估计，可以从图左侧的点分布的缩小中看到。使用\Rfunction{identify}函数挑选出MA图中单个基因的结果。第一个参数应该是\Robject{res\$baseMean}，第二个参数应该是\Robject{res\$log2FoldChange}。按\ emh {Esc}完成。R会打印出你选择的基因的行号。 Use this to index the \Robject{res} object. \end{exer} This plot demonstrates that only genes with a large average normalized count contain sufficient information to yield a significant call. Also note \deseqtwo{}'s shrinkage estimation of log fold changes (LFCs): When count values are too low to allow an accurate estimate of the LFC, the value is ``shrunken'' towards zero to avoid that these values, which otherwise would frequently be unrealistically large, dominate the top-ranked log fold changes. Whether a gene is called significant depends not only on its LFC but also on its within-group variability, which \deseqtwo{} quantifies as the \emph{dispersion}. For strongly expressed genes, the dispersion can be understood as a squared coefficient of variation: a dispersion value of 0.01 means that the gene's expression tends to differ by typically $\sqrt{0.01}=10\%$ between samples of the same treatment group. For weak genes, the Poisson noise is an additional source of noise, which is added to the dispersion. The function \Rfunction{plotDispEsts} visualizes \deseqtwo's dispersion estimates: < >= plotdis虫害(dds, ylim = c(1e-6, 1e1)) @ \incfig{figure/ beginner_displot}{。5\textwidth}{离散度估计值图}{详情见文本}黑点是每个基因的离散度估计值，分别考虑来自每个基因的信息。除非有很多样本，否则这些值会在真实值附近剧烈波动。因此，我们拟合红色的趋势线，它显示了分散对均值的依赖，然后将每个基因的估计缩小到红线，以获得最终的估计(蓝点)，然后用于假设检验。在主要“云”点上方的蓝色圆圈是具有高基因分散估计的基因，标记为分散离群值。因此，这些估计值并没有向拟合的趋势线收缩。图/ beginner_pvalHist \ incfig{}{。另一个有用的诊断图是$p$值的直方图(图~\ref{figure/beginner_pvalHist})。< >= hist(res$pvalue, breaks=20, col="grey") @ \begin{exer} p值直方图有一些峰值，在1处有一个大峰值。如果你只绘制那些\Robject{res\$baseMean}大于$10$的基因的直方图会发生什么?MA图(图~\ref{图/beginner_MA})突出了RNA-Seq数据的一个重要属性。对于弱表达的基因，我们没有机会看到差异表达，因为低读计数受到很高的泊松噪声的影响，任何生物效应都淹没在读计数的不确定性中。我们还可以通过检查被平均归一化计数分类的基因的小$p$值(例如小于0.01)的比率来证明这一点:< >= #创建垃圾箱使用分位数函数qs <- c(0，分位数(res$baseMean[res$baseMean > 0]， 0:7/7)) #“切割”基因到垃圾箱垃圾箱<- cut(res$baseMean, qs) #重命名的级别的垃圾箱使用中间点水平(垃圾箱)<- paste0("~"，圆(。5*qs[-1] + .5*qs[-length(qs)])) #计算每个bin ratio <- tapply(res$pvalue, bins, function(p) mean(p < .01, na. 0)的p值小于.01的比值。rm=TRUE)) #绘制这些比率barplot(比率，xlab=“平均归一化计数”，ylab=“小p值的比率”)@ \incfig{图/beginner_ratioSmallP}{。9\textwidth}{小p值的比值}{对于被平均归一化计数分类的基因组}乍一看，过滤这些基因似乎没有什么好处。毕竟，测试发现它们并不显著。但是，这些基因对多重测试调整有影响，如果去除这些基因，多重测试调整的性能会提高。通过从FDR程序的输入中删除弱表达基因，我们可以在保留的基因中发现更多重要的基因，从而提高了测试的能力。这种方法被称为\emph{独立过滤}。\deseqtwo{}软件自动执行独立的过滤，最大限度地增加将调整$p$值小于临界值的基因数量(默认情况下，\Robject{alpha}被设置为0.1)。这种自动独立的过滤由\Rfunction{results}函数执行，并且可以由\Rfunction{results}函数控制。我们可以观察到基于平均归一化计数的各种截止点的拒绝数量是如何变化的。下面的最佳阈值和可能值的表存储为结果对象的\emph{属性}。 < >= attr(res，"filterThreshold") plot(attr(res，"filterNumRej")，type="b"， xlab="quantiles of 'baseMean'"， ylab="number of rejections") @ \incfig{figure/beginner_filtByMean}{。5 \ textwidth}{独立过滤。{\deseqtwo{}自动确定一个阈值，对平均规范化计数进行过滤，使调整后的$p$值小于临界值的基因数量最大化。“独立”一词强调了一个重要的警告。只有当过滤条件独立于实际测试统计量时，才允许这样的过滤~\cite{Bourgon:2010:PNAS}。否则，过滤将使测试无效，从而使BH过程的假设无效。这就是为什么我们对所有样本的平均值进行过滤:这个过滤器对处理组和对照组的样本分配是盲目的，因此是独立的。\分段{导出结果}最后，我们注意到，您可以轻松地将结果表保存在CSV文件中，然后您可以使用电子表格程序(如Excel)加载:< >= write.csv(as.data.frame(res)， file="results.csv") @ \分段{基因集富集分析}基因上调或下调的基因是否有共同之处?我们接下来进行基因集富集分析(GSEA)来检查这个问题。正如讲座中提到的，利用RNA-Seq数据进行基因集富集分析需要一些微妙之处。简单地说，已经提出了许多不同的方法，每一种方法都暗示着被检验的零假设略有不同，它们的生物学解释也不同。这是一个正在进行的研究课题。我们在这里提出了一个相对简单的方法，以证明基本思想，但注意，为了获得更明确的结果，需要更仔细的处理。我们使用Reactome数据库中的基因集。 >= library(" reacome .db") @这个数据库使用Entrez id，所以我们将需要之前添加到\Robject{res}对象中的\Robject{entrezid}列。首先，我们将结果表\Robject{res}子集化为Reactome数据库中有数据的基因(即，我们在\Robject{reacome .db}的相应键列中找到其Entrez ID)，并且\deseqtwo{}测试给出的调整后的$p$值不是\Robject{NA}。< >= res2 <- res[res$entrezid %in% keys(reactme .db， " entrezid ") & !使用\Rfunction{select}，一个来自\Rpackage{AnnotationDbi}的用于查询数据库对象的函数，我们得到一个从Entrez id到Reactome Path id的映射表。 >= reactomeTable <- AnnotationDbi::select(reactomeTable, keys=as.character(res2$entrezid)， keytype=" entrezid "， columns=c(" entrezid "，"REACTOMEID")) head(reactomeTable) @ %下一个代码块将这个表转换为\emph{关联矩阵}。这是一个布尔矩阵，其中一行表示每个反应组路径，一列表示\Robject{res2}中的每个唯一基因，它告诉我们哪些基因属于哪些反应组路径。(如果你想了解这个块是如何工作的，请阅读\Rfunction{tapply}函数。)%这段代码的问题是res2$entrez不是唯一的% < >= % ft2mat <- function(x, y) {% ux = unique(x) % uy = unique(y) % im = matrix(FALSE, nrow=length(ux)， ncol=length(uy)， dimnames=list(ux, uy)) % im[cbind(x, y)] = TRUE % return(im) %} % incm <- with(reactomeTable, ft2mat(REACTOMEID, ENTREZID)) % @ < >= incm <- do。call(rbind, with(reactomeTable, tapply(ENTREZID, factor(REACTOMEID)， function(x) res2$ ENTREZID %in% x)) colnames(incm) <- res2$entrez str(incm) @我们删除所有对应于Reactome路径的小于20或大于80个已分配基因的行。< >= within <- function(x, lower, upper) (x>=lower & x<=upper) incm <- incm[within(rowsum (incm)， lower=20, upper=80)，] @为了测试在我们的实验中，Reactome路径中的基因是否以特殊方式表现，我们计算了一些统计量，包括$t$-statistic，以查看基因集中基因的log2倍变化值的平均值是否不同于零。为了便于计算，我们定义了一个小的helper函数:< >= testCategory <- function(reactomeID) {isMember <- incm[reactomeID，] data.frame(reactomeID = reactomeID, numGenes = sum(isMember)， avgLFC = mean(res2$log2FoldChange[isMember])， sdLFC = sd(res2$log2FoldChange[isMember])， zValue = mean(res2$log2FoldChange[isMember]) / sd(res2$log2FoldChange[isMember])， strength = sum(res2$log2FoldChange[isMember]) / sqrt(sum(isMember))， pvalue = t.test(res2$log2FoldChange[isMember])$p。value, reactomeName = reactomePATHID2NAME[[reactomeID]]， stringsAsFactors = FALSE)} @该函数可以用一个reactompath ID来调用:< 正如你所看到的，这个函数不仅执行$t$测试并返回p值，而且还列出了其他有用的信息，如类别中的基因数量、平均对数折叠变化、“强度”测量(见下文)和Reactome描述路径的名称。我们为关联矩阵中的所有路径调用函数，并在数据帧中收集结果:% < >= reactomeResult <- do。由于我们执行了许多测试，我们应该再次使用多个测试调整。% < >= reactomeResult$padjust <- p.adjust(reactomeResult$pvalue， "BH") @这是一个Reactome路径列表，这些路径在我们的DPN处理与对照组的比较中有显著差异，根据信号的符号和强度排序:{\小< >= reactomeResultSignif <- reactomeresultesult $padjust < 0.05，] head(reactomeResultSignif[order(-reactomeResultSignif$strength)，]) @}然而，正如演讲中所讨论的，在这里使用t检验是非常值得怀疑的。毕竟，一组基因通常是相关的，这违反了t检验核心的独立性假设。因此，我们真的应该看t检验的p值吗?通过样本置换得到的p值可以解决这个问题，因为样本置换保留了基因-基因相关性。然而，由于只有四个受试者，我们没有足够的样本。因此，更谨慎的做法可能是完全忽略这些可疑的p值，而是查看更原始的统计数据，例如一条路径内的平均LFC，可能除以标准偏差。我们将其存储为\Robject{zValue}。{\ <小 如果这样的分析只被认为是探索性的，我们可以检查各种这样的表，看看路径的排名是否有助于我们理解数据。然而，这方面肯定有改进的空间。按$z$ value排名的前几个点击都有完全相同的值。为什么?用于多维数据探索性分析的许多常用统计方法，特别是用于聚类和排序(e.\ \，g.)的方法。，主成分分析等)，最适用于(至少近似)同方差数据;这意味着一个可观察量的方差(即，在这里，基因的表达强度)不依赖于平均值。然而，在RNA-Seq数据中，方差随着均值的增加而增加。例如，如果直接对归一化读取计数的矩阵进行主成分分析，结果通常只依赖于少数最强表达的基因，因为它们显示了样本之间最大的绝对差异。避免这种情况的一个简单且常用的策略是取规范化计数值的对数加上一个小的伪计数;然而，现在低计数的基因倾向于主导结果，因为由于小计数值固有的强泊松噪声，它们在样本之间表现出最强的相对差异。 As a solution, \deseqtwo{} offers the \emph{regularized-logarithm transformation}, or \emph{rlog} for short. For genes with high counts, the rlog transformation differs not much from an ordinary log2 transformation. For genes with lower counts, however, the values are shrunken towards the genes' averages across all samples. Using an empirical Bayesian prior in the form of a \emph{ridge penalty}, this is done such that the rlog-transformed data are approximately homoskedastic. Note that the rlog transformation is provided for applications other than differential testing. For differential testing we recommend the \Rfunction{DESeq} function applied to raw counts, as described earlier in this vignette, which also takes into account the dependence of the variance of counts on the mean value during the dispersion estimation step. The function \Rfunction{rlog} returns a \Rclass{SummarizedExperiment} object which contains the rlog-transformed values in its \Rclass{assay} slot: < >= rld <- rlog(dds) assay(rld)[1:3, 1:3] @为了显示转换的效果，我们首先简单地使用\Rfunction{log2}函数(在添加1之后，以避免取零的对数)绘制第一个样本与第二个样本的图，然后使用经过rlog转换的值。图/ beginner_rldPlot \ incfig{}{。9\textwidth}{样本2 vs样本1的散射图。}{左:使用普通的log2转换。右:使用rlog变换。} < >= par(mfrow = c(1,2)) plot(log2(1+counts(dds，归一化=TRUE)[， 1:2])， col="#00000020"， pch=20, cex=0.3) plot(assay(rld)[， 1:2]， col="#00000020"， pch=20, cex=0.3) @注意，为了更容易看到几个点在哪里相互重叠，我们将绘图颜色设置为半透明的黑色(编码为\Rfunction{\#00000020})，并将点更改为固体磁盘(\Rfunction{pch=20})，减小尺寸(\Rfunction{cex=0.3})\脚注{来自包\CRANpkg{LSD}的函数\Rfunction{heatscatter}提供了一个彩色的替代方案。在图~\ref{Figure /beginner_rldPlot}中，我们可以看到，在普通对数尺度上，低计数的基因似乎是过度可变的，而rlog变换压缩了数据无论如何都不能提供良好信息的基因的差异。\分段{样本距离}在RNA-Seq分析中有用的第一步通常是评估样本之间的总体相似性:哪些样本彼此相似，哪些不同?这符合实验设计的期望吗?图/ beginner_sampleDistHeatmap \ incfig{}{。5\textwidth}{rlog变换后的欧几里得样本距离热图。我们使用R函数\Rfunction{dist}来计算样本之间的欧氏距离。为了避免距离测量被少数高度可变的基因所主导，并且所有基因的贡献大致相等，我们在rlog转换数据上使用它:< >= sampledist <- dist(t(assay(rld))) as。矩阵(sampleDists)[1:3, 1:3] @注意使用函数\Rfunction{t}来转置数据矩阵。我们需要这样做是因为\Rfunction{dist}计算data \ emh {rows}和我们的样本构成列之间的距离。我们在热图中可视化距离，使用函数\Rfunction{热图。2}从\CRANpkg{gplots}包。< >= sampleDistMatrix <- as。matrix(sampleDists) rownames(sampleDistMatrix) <- paste(rld$treatment, rld$patient, sep="-") colnames(sampleDistMatrix) <- NULL library("gplots") library("RColorBrewer") colors = colorrampalalette (rev(brewer。pal(9，“Blues”))(255)热图。2(sampleDistMatrix, trace="none"， col= colors) @注意，我们已经改变了距离矩阵的行名，以包含治疗类型和患者编号，而不是样本ID，这样我们在查看热图时就可以看到所有这些信息(图.~\ref{figure/beginner_sampleDistHeatmap})。图/ beginner_samplePCA \ incfig{}{。6\textwidth}{主成分分析(PCA)}{经过rlog变换的样本。另一种可视化样本到样本距离的方法是主成分分析(PCA)。在这种排序方法中，数据点(即这里的样本)被投影到2D平面上，以便它们以最佳方式展开(图.\ \ref{figure/beginner_samplePCA})。< >= colors <- c(rgb(1:3/4,0,0)，rgb(0,1:3/4,0)，rgb(0,0,1:3/4)，rgb(1:3/4,0,1:3/4)) plotPCA(rld, intgroup =c ("patient"，"treatment")， col= colors) @这里，我们使用了\Rfunction{plotPCA}，它附带\deseqtwo{}。\Rfunction{intgroup}指定的两个项是样本数据中的列名;它们告诉函数使用它们来选择颜色。从这两个可视化图中，我们可以看到患者之间的差异远大于同一患者的治疗样本和对照样本之间的差异。这说明了为什么考虑这种配对设计是很重要的(“配对”，因为每个处理的样本都与来自同一患者的一个对照样本配对)。我们在一开始设置数据对象时使用设计公式\texttt{\lowtilde{} patient + treatment}来实现这一点。如果我们使用非配对分析，通过只指定\texttt{\lowtilde{}治疗，我们将不会发现很多匹配，因为这样，患者之间的差异将淹没任何治疗效果。在这里，我们在分析结束时进行了样本距离分析。然而，在实践中，这是一个适合于对数据进行初步概述的步骤。因此，人们通常会将此分析作为分析的第一步之一。 To this end, you may also find the function \Rfunction{arrayQualityMetrics}, from the package of the same name, useful. \subsection{Gene clustering} \incfig{figure/beginner_geneHeatmap}{.5\textwidth}{Heatmap with gene clustering.} In the heatmap of Fig.~\ref{figure/beginner_sampleDistHeatmap}, the dendrogram at the side shows us a hierarchical clustering of the samples. Such a clustering can also be performed for the genes. Since the clustering is only relevant for genes that actually carry signal, one usually carries it out only for a subset of most highly variable genes. Here, for demonstration, let us select the 35 genes with the highest variance across samples: < >= library("genefilter") topVarGenes <- head(order(rowVars(assay(rld))，递减=TRUE)， 35) @如果我们不看绝对表达强度，而是看特定样本中每个基因偏离所有样本中基因平均值的量，热图就会变得更有趣。因此，我们在样本中对每个基因的值进行集中和缩放，并绘制热图。< > =热图。2(assay(rld)[topVarGenes，]， scale="row"， trace="none"， dendrogram="column"， col = colorRampPalette(rev(brewer。pal(9， "RdBu")))(255)， ColSideColors = c(Control="灰色"，DPN="深绿色"，OHT="橙色")[colData(rld)$treatment]) @ % \begin{exer} % \Rfunction{热图的哪个选项。2}导致基因的定心和缩放?如果不缩放热图，%热图会发生什么变化?我们现在可以看到(图. \ref{figure/beginner_geneHeatmap})在患者之间共变的基因块。通常，这样的热图是有洞察力的，即使在这里，看到患者之间的这些变化是有限的价值，因为我们更感兴趣的是每个患者之间的治疗效果。这张热图显示了不同患者之间差异最大的基因。在我们的差异表达分析中，我们发现了一系列基因在治疗组和对照组之间有显著差异。在类似的热图中显示这些基因(或者可能只是那些在治疗中显示上调调节的基因)。能否确认热图目测的测试结果? \end{exer} \begin{exer} In the abstract of the publication for this dataset, 5 parathyroid-related genes are mentioned as differentially expressed due to OHT or DPN treatment in the 48 hour samples compared to control: \Robject{c("CASR","VDR","JUN","CALR","ORAI2")}. Using the results table for OHT vs Control and the results table for DPN vs Control, find the rank of these genes in terms of smallest p values using the \Rfunction{match} function. \end{exer} \clearpage \section{Going further} \begin{exer} Analyze more datasets: use the function defined in the following code chunk to download a processed count matrix from the ReCount website\footnote{\url{http://bowtie-bio.sourceforge.net/recount/}}. Suggested datasets: \begin{itemize} \item \emph{hammer} -- spinal nerve ligation of rats at 2 time points (note a misspelling in one of the time points in \Robject{colData}). Perform differential expression and describe the most significant gene sets. Do these have something to do with spinal nerve ligation? \item \emph{wang} -- human tissue comparison. Make a PCA plot of the samples \item \emph{bottomly} -- 2 inbred mouse strains. Make a PCA plot of the samples \end{itemize} Note that the annotation libraries for rat and mouse are \Biocannopkg{org.Rn.eg.db} and \Biocannopkg{org.Mm.eg.db}. \end{exer} The following function takes a name of the dataset from the ReCount website, e.g. \Robject{"hammer"}, and returns a \Rclass{SummarizedExperiment} object. < >= recret2se <- function(name) {filename <- paste0(name，" _set . rdata ") if (!file.exists(filename))下载。file(paste0("http://bowtie-bio.sourceforge.net/recount/ExpressionSets/"， filename)，filename) load(filename) e <- get(paste0(name，".eset")) se <- summarizeexperimental (SimpleList(counts=exprs(e))， colData=DataFrame(pData(e))) se} @ \clearpage \section{另见}\begin{itemize} \item \Biocpkg{DEXSeq}用于不同外显子的使用。参见附带的小插图，\emph{分析RNA-seq数据的差异外显子使用DEXSeq包}，这是类似于本教程的风格。\item其他用于RNA-Seq差异表达的Bioconductor软件包:\Biocpkg{edgeR}， \Biocpkg{limma}， \Biocpkg{DSS}， \Biocpkg{BitSeq}(转录水平)，\Biocpkg{EBSeq}， \Biocpkg{cummeRbund}(用于导入和显示袖扣结果)，\Biocpkg{monocle}(单细胞分析)。更多在\url{//www.andersvercelli.com/packages/release/BiocViews.html#___RNASeq} \item基因集测试方法:\Rfunction{romer}和\Rfunction{roast}在\Biocpkg{limma}，基于排列:\Biocpkg{safe} \item协变量归一化包(例如，GC内容):\Biocpkg{cqn}， \Biocpkg{EDASeq} \item检测批次包:\Biocpkg{sva}， \Biocpkg{RUVSeq} \item生成与外部资源链接的HTML结果表(基因描述):\Biocpkg{ReportingTools} \end{itemize} \bibliography{RNA-Seq-Analysis-Lab} \section{Session Info}作为本文档的最后一部分，我们调用函数\Rfunction{sessionInfo}，它报告R的版本号和在此会话中使用的所有包。始终保持这样的记录是一个很好的实践，因为它将有助于追踪发生了什么情况下，R脚本停止工作，因为一个包已经在一个新的版本中被更改。会话信息也应始终包含在任何电子邮件到Bioconductor邮件列表。< >= toLatex(sessionInfo()) @ %从摘要中查找基因的代码:% resOHT$hgnc_symbol <- convertIDs(line .names(resOHT)， "ENSEMBL"， "SYMBOL"， org.Hs.eg.db) % genes <- c("CASR"，"VDR"，"JUN"，"CALR"，"ORAI2") % match(genes, res$hgnc_symbol[order(res$pvalue)]) % match(genes, resOHT$hgnc_symbol[order(resOHT$pvalue)]) \end{document}