From db76942b61d8f1a651a83b71623d8ca47ef88afc Mon Sep 17 00:00:00 2001 From: Xiao W Date: Fri, 28 Jun 2019 00:02:10 -0400 Subject: [PATCH 01/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 158 ++++++++++++------------- 1 file changed, 79 insertions(+), 79 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index a443f6c..6337450 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -7,7 +7,7 @@ author: affiliation: EMBL Heidelberg, Meyerhofstrasse 1, 69117 Heidelberg, Germany, steffi.reisenauer@tum.de header-includes: - \usepackage{float} -abstract: In this article, we walk through an end-to-end Affymetrix microarray differential expression workflow using Bioconductor packages. This workflow is directly applicable to current "Gene" type arrays, e.g. the HuGene or MoGene arrays, but can easily be adapted to similar platforms. The data analyzed here is a typical clinical microarray data set that compares inflamed and non-inflamed colon tissue in two disease subtypes. For each disease, the differential gene expression between inflamed- and non-inflamed colon tissue was analyzed. We will start from the raw data CEL files, show how to import them into a Bioconductor ExpressionSet, perform quality control and normalization and finally differential gene expression (DE) analysis, followed by some enrichment analysis. +摘要: 在这篇文章里，我们将阐述用 Bioconductor 包 an end-to-end Affymetrix microarray differential expression 的工作流程。这个工作流程可以直接应用在目前 “基因” 类型的微陈列，例如 HuGene 或 MoGene 微陈列，同时也容易调整并适用于类似的平台。在这里进行分析的数据是来自典型的临床微陈列数据集，用来对比两种疾病亚型中发炎和非发炎的结肠组织。针对每种疾病，发炎和非发炎的结肠组织的差异基因表达被进行了分析。我们将从原始数据 CEL 文件开始，展示如何将其导入为 Bioconductor 的 ExpressionSet，实施质量控制和归一化和最终的差异基因表达分析，其次是一些显著性分析。 keywords: Microarray, Gene Expression bibliography: MAEndToEnd.bib output: @@ -94,39 +94,39 @@ opts_chunk$set(echo = TRUE,

-**R version**: `r R.version.string` +**R 版本**: `r R.version.string`
-**Bioconductor version**: `r BiocManager::version()` +**Bioconductor 版本**: `r BiocManager::version()`
-**Package version**: `r packageVersion("maEndToEnd")` +**包版本**: `r packageVersion("maEndToEnd")`

-# Introduction +# 介绍 -In this article we introduce a complete workflow for a typical (Affymetrix) -microarray analysis. Data import, preprocessing, differential expression and -enrichment analysis are discussed. +在这篇文章里我们将介绍典型的（Affymetrix）微阵列分析的完整 +工作流程。从数据导入，前期处理，差异基因表达和显著性分析都会 +被讨论。 -The data set used [@Palmieri_2015] is from a paper studying the differences in -gene expression in inflamed and non-inflamed tissue. -14 patients suffering from Ulcerative colitis (UC) and 15 patients with Crohn's -disease (CD) were tested, and from each patient inflamed and non-inflamed -colonic mucosa tissue was obtained via a biopsy. -This is a typical clinical data set consisting of 58 arrays in total. -Our aim is to analyze differential expression (DE) between the tissues. -Our results show a substantial overlap with the results of the original paper. +使用的数据 [@Palmieri_2015] 来自于一篇研究发炎和非发炎组 +织中差异基因表达的论文。 +测试了 14 个患有溃疡性结肠炎 (UC)的病人和 15 个患有克罗恩病 +（CD) 的病人，然后通过活检提取了每个病人发炎和未发炎的结肠粘 +膜组织。 +这个典型的临床数据集共由 58 个阵列构成。我们的目标是分析组织 +间的差异表达（DE）。 +我们的结果显示了与原始论文结果的大量重叠。 -# Workflow package installation +# 工作流程包的安装 -The workflow is wrapped in a package called `maEndToEnd`. +工作流程被包裹在叫做 `maEndToEnd` 的包里。 -The `maEndToEnd` package can currently be obtained from GitHub and is -available via the current development version of Bioconductor (3.8) (see here: ). +`maEndToEnd` 包目前可以从 GitHub 获得， +目前可用的开发版本为 Bioconductor (3.8) (看这里:). -## Workflow package installation from Bioconductor +## 从 Bioconductor 安装工作流程包 -You can install the package via the `r CRANpkg("BiocManager")`. +你可以通过 `r CRANpkg("BiocManager")`安装包。 ```{r installBioc, eval=FALSE} if (!require("BiocManager")) @@ -134,22 +134,22 @@ if (!require("BiocManager")) BiocManager::install("maEndToEnd", version = "devel") ``` -Currently, the workflow is available in the development version of Bioconductor -(3.8), which will become the release version in October 2018. +目前，工作流程可在 Bioconductor（3.8）的开发版本中获得， 这将 +成为 2018 年 10 月的发布版本。 -For details on how to use this version of Bioconductor see: +有关如何使用这个版本 Bioconductor 的详细信息见： -## Workflow package installation from Github +## 从 Github 安装工作流程包 ```{r, echo=TRUE, results="hide", warning=FALSE, eval=FALSE} -#In order to download the package from GitHub, we need the "install_github" -#function from the "remotes" package. We download the latest developer -#version of "remotes" from GitHub with the devtool::install_github -#function; note that this is necessary as the current "remotes" version on -#CRAN doesn't allow us to correctly download the "maEndToEnd" package: +#为了从 GitHub下载包，我们需要来自 “remotes” 的包里中 +#“install_github”的功能。我们用 +#devtool::install_github 从 GitHub下载了最新的开发 +#版本的 "remotes" 函数；注意：这是必要的。 +#目前在CRAN上面的 “remotes” 版本不允许我们正确的下载 "maEndToEnd"包： install.packages("devtools") library(devtools) @@ -162,56 +162,56 @@ remotes::install_github("b-klaus/maEndToEnd", ref="master") ``` -## Workflow package import +## 工作流程包的导入 -Once the workflow package has been successfully installed, we can use -a call to `library()` in order to load it. This will also load -all the other packages neccessary to run the workflow. +成功安装工作流程包后， +我们可以使用对library() 的调用来加载它。 +它还将加载运行工作流程所需的所有其他包。 ```{r maEndToEndImport} suppressPackageStartupMessages({library("maEndToEnd")}) ``` -## List of packages required for the workflow +## 工作流程所需要的包列表 -Below, you find a list of packages that are required by the workflow. -Some Helper/Styling packages have been commented here, as they are -not strictly neccesary to execute the workflow. +下面，您可以找到工作流程所需的软件包列表。 +这里已经注释了一些 Helper/Styling 包， +它们不是执行工作流所必需的。 ```{r pkgList, results="hide"} -#General Bioconductor packages +#通常的 Bioconductor 包 library(Biobase) library(oligoClasses) -#Annotation and data import packages +#注释和数据导入包 library(ArrayExpress) library(pd.hugene.1.0.st.v1) library(hugene10sttranscriptcluster.db) -#Quality control and pre-processing packages +#质量控制和前期处理包 library(oligo) library(arrayQualityMetrics) -#Analysis and statistics packages +#分析和统计包 library(limma) library(topGO) library(ReactomePA) library(clusterProfiler) -#Plotting and color options packages +#绘图和颜色选项包 library(gplots) library(ggplot2) library(geneplotter) library(RColorBrewer) library(pheatmap) -#Formatting/documentation packages +#格式化/文档包 #library(rmarkdown) #library(BiocStyle) library(dplyr) library(tidyr) -#Helpers: +#助手包 (Helpers) library(stringr) library(matrixStats) library(genefilter) @@ -221,14 +221,14 @@ not strictly neccesary to execute the workflow. ``` -# Downloading the raw data from ArrayExpress +# 从 ArrayExpress 下载原始数据 -The first step of the analysis is to download the raw data CEL files. These -files are produced by the array scanner software and contain the measured probe intensities. The data we use have been deposited at [ArrayExpress](https://www.ebi.ac.uk/arrayexpress/) -and have the accession code **E-MTAB-2967**. +分析的第一步是下载原始数据CEL文件。 +这些文件由阵列扫描仪软件生成，包含测量的探针强度。我们使用的数据已存放在 [ArrayExpress](https://www.ebi.ac.uk/arrayexpress/) 上， +其登记号为 **E-MTAB-2967**. -We will store these files in the directory **raw\_data\_dir** which defaults to -a temporary directory. +我们将这些文件存储在 **raw\_data\_dir** 目录中， +该目录默认为临时目录。 ```{r generateFolderForRawData} raw_data_dir <- tempdir() @@ -240,48 +240,48 @@ if (!dir.exists(raw_data_dir)) { ``` -Each ArrayExpress data set has a landing page summarizing the data set, -and we use the `getAE`function from the `r Biocpkg("ArrayExpress") ` Bioconductor package to obtain -the ftp links to the raw data files ([Data from Palmieri et. al. on ArrayEpress](https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2967/)). +每个 ArrayExpress 数据集都有一个总结数据集的登陆页面， +我们使用`r Biocpkg("ArrayExpress") ` Bioconductor 包中的 `getAE` 函数获取原始数据文件的 +ftp 链接 ([Data from Palmieri et. al. on ArrayEpress](https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2967/)). -With the code below, we download the raw data (also including annotation data) from -[ArrayExpress](https://www.ebi.ac.uk/arrayexpress/) [@Kolesnikov_2014] -by using the `getAE`-function. The data are saved in the `raw_data_dir` -created above. The names of the downloaded files are returned as a list. +使用下面的代码，我们使用`getAE`函数 +从 [ArrayExpress](https://www.ebi.ac.uk/arrayexpress/) [@Kolesnikov_2014] 下载原始数据（也包括注释数据）。 +数据保存在上面创建的 `raw_data_dir`中。 +下载文件的名称将作为列表返回。 ```{r getDataEBI, eval=TRUE, results='hide', message=FALSE} anno_AE <- getAE("E-MTAB-2967", path = raw_data_dir, type = "raw") ``` -We will now have a closer look at the data we downloaded from ArrayExpress +我们现在将仔细研究从 ArrayExpress 下载的数据。 -# Background information on the data +# 有关数据的背景信息 -## Information stored in ArrayExpress +## ArrayExpress 里储存的信息 -Each dataset at ArrayExpress is stored according to the MAGE-TAB -(MicroArray Gene Expression Tabular) specifications as a collection of -tables bundled with the raw data. The MAGE-TAB format specifies up to five -different types of files: +ArrayExpress 中的每个数据集都根据 +MAGE-TAB（微阵列基因表达式表格）规范存储为与 +原始数据捆绑在一起的表的集合。 +MAGE-TAB 格式指定最多五种不同类型的文件： -* **Investigation Description Format (IDF)** -* **Array Design Format (ADF)** -* **Sample and Data Relationship Format (SDRF)** -* **raw data files** -* **processed data files** +* **调查说明格式 (IDF)** +* **阵列设计格式 (ADF)** +* **样本和数据关系格式 (SDRF)** +* **原始数据文件** +* **处理过的数据文件** -Other than the raw data files, the IDF and the SDRF file are important for us. -The IDF file contains top level information -about the experiment including title, description, submitter contact details -and protocols. The SDRF file contains essential information on the experimental -samples, e.g. the experimental group(s) they belong to. +除原始数据文件外， +IDF 和 SDRF 文件对我们很重要。 +IDF 文件包含有关实验的顶级信息， +包括标题，说明，提交者联系详细信息和协议。 +SDRF 文件包含有关实验样本的基本信息，例如他们所属的实验组。 -Before we move on to the actual raw data import, we will briefly introduce the -`ExpressionSet` class contained in the `r Biocpkg("Biobase")` package. -It is commonly used to store microarray data in Bioconductor. +在我们继续进行实际的原始数据导入之前， +我们将简要介绍`r Biocpkg("Biobase")` 包中包含的 `ExpressionSet` 类。 +它通常用于在 Bioconductor 中存储微阵列数据。 ## Bioconductor ExpressionSets From 48628b656759706711b486a0e7cedae7fafc66d5 Mon Sep 17 00:00:00 2001 From: Xiao W Date: Fri, 28 Jun 2019 11:37:51 -0400 Subject: [PATCH 02/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 432 ++++++++++++------------- 1 file changed, 216 insertions(+), 216 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index 6337450..db6ec4f 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -285,37 +285,37 @@ SDRF 文件包含有关实验样本的基本信息，例如他们所属的实验 ## Bioconductor ExpressionSets -Genomic data can be very complex, -usually consisting of a number of different components, e.g. information -on the experimental samples, annotation of genomic features measured as well -as the experimental data itself. -In Bioconductor, the approach is taken that these components should be -stored in a single structure to easily manage the data. - -The package `r Biocpkg("Biobase")` contains standardized data structures -to represent genomic data. The `ExpressionSet` class is designed -to combine several different sources of information (i.e. as contained in the -various MAGE-TAB files) into a single convenient -structure. An ExpressionSet can be manipulated (e.g., subsetted, copied), -and is the input to or output of many Bioconductor functions. - -The data in an ExpressionSet consist of: - -* **assayData**: Expression data from microarray experiments -with microarray probes in rows and sample identifiers in columns +基因组数据可以非常复杂， +通常由许多不同的组分组成，例如， +有关实验样品的信息，测量的基因组特征的注释 +以及实验数据本身。 +在 Bioconductor 中， 采用的方法是将这些组件存储 +在一个单结构中以便于管理数据。 + +`r Biocpkg("Biobase")` 包含标准化数据结构以表示基因组数据。 +`ExpressionSet` 类旨在将几个不同的信息源 +（即包含在各种 MAGE-TAB 文件中）组合 +成一个方便的结构。 +可以操纵 ExpressionSet（例如，子集化，复制） +并且是许多 Bioconductor 函数的输入或输出。 + +ExpressionSet 中的数据包括： + +* **assayData**: 来自微阵列实验的表达数据，其中每行是微阵列探针，每列是样品标识符。 +其中每行是微阵列探针，每列是样品标识符。 * **metaData** - + **phenoData**: A description of the samples in the experiment with sample identifiers in rows and description elements in columns; holds the content of the SDRF file - + **featureData**: metadata about the features on the chip - or technology used for the experiment with same rows as assayData by default and freely assignable columns - + further annotations for the features, - for example gene annotations from biomedical databases (annotation). + + **phenoData**: 实验中样本的描述，每行是样本标识符，每列是描述元素; 储存的是 SDRF 文件的内容。 + + **featureData**: 有关芯片功能的元数据或者用于实验的技术， + 默认情况下使用与 assayData 相同的行和可自由分配的列。 + + 进一步注释特征，例如来自生物医学数据库（注释）的基因注释。 + 例如来自生物医学数据库（注释）的基因注释。 -* **experimentData**: A flexible structure to describe the experiment. +* **experimentData**: 一种灵活的结构来描述实验。 -The ExpressionSet class coordinates all of these data, so that one does not -have to worry about the details. However, one should keep in mind that the rownames of the `phenoData` have to match the column names of the assay data, while the row names of the assay data have to match the -row names of the `featureData`. -This is illustrated in Figure \@ref(fig:sumexp). +ExpressionSet 类协调所有这些数据，因此不必 +担心细节。 但是，应该记住，`phenoData` 的 rownames （行名）必须与 assay data 的列名匹配， +而 assayData 的行名必须与 `featureData` 的行名匹配。 +这在图1中说明。 \@ref(fig:sumexp). ```{r sumexp, fig.cap = "Structure of Bioconductor’s ExpressionSet class.", echo=FALSE, fig.show="asis", fig.width = 7, fig.height = 4.5} @@ -354,20 +354,20 @@ polygon(c(50,65, 65, 50),c(56,56,65,65),col="darkseagreen3",border=NA) text(57.5,60.5,"sample IDs", cex = 1) ``` -You can use the functions ` pData ` and ` fData ` to extract -the sample and feature annotation, respectively, from an ` ExpressionSet `. -The function ` exprs ` will return the expression data itself as a matrix. +您可以使用函数 ` pData ` 和 ` fData ` 分别从 +` ExpressionSet ` 中提取样本和要素注记。 +函数 ` exprs ` 将表达式数据本身作为矩阵返回。 -# Import of annotation data and microarray expression data as "ExpressionSet" +# 导入注释数据和微阵列表达数据为 "ExpressionSet" -We import the SDRF file with the `read.delim` function from the raw data -folder in order to obtain the sample annotation. +我们使用原始数据文件夹中的 `read.delim` 函数导入 +SDRF 文件，以获取样本注释。 -The sample names are given in the column Array.Data.File of the SDRF data table and will be used as rownames for the SDRF file. +样本名称在 SDRF 数据表的 Array.Data.File 列中给出，并将用作 SDRF 文件的 rownames 。 -We turn the SDRF table into an `AnnotatedDataFrame` -from the `r Biocpkg("Biobase")` package that we -will need later to create an `ExpressionSet` for our data [@Bioc]. +我们将 SDRF 表转换为 `r Biocpkg("Biobase")` 包中 +的 `AnnotatedDataFrame` ，稍后我们将需要为我们的数据 +创建一个 `ExpressionSet`[@Bioc]. ```{r getSDRF} sdrf_location <- file.path(raw_data_dir, "E-MTAB-2967.sdrf.txt") @@ -378,19 +378,19 @@ SDRF <- AnnotatedDataFrame(SDRF) ``` -We now create the Expression Set object -`raw_data`, which contains array data, pheno data -(from the SDRF file) as well as the information -of the chip annotation package used. +我们现在创建 Expression Set object +`raw_data`,它包含陈列数据，pheno data (表型数据) +（来自 SDRF 文件） 以及所使用 +的芯片注释包的信息。 -The analysis of Affymetrix arrays starts with CEL files. These are the result -of the processing of the raw image files using the Affymetrix software -and contain estimated probe intensity values. Each CEL file additionally -contains some metadata, such as a chip identifier. +Affymetrix 列阵的分析是从 CEL 文件开始。这些是使用 +Affymetrix 软件处理原始图像文件 +并包含估计的探针强度值的结果。每个 CEL 文件还包含一些元数据， +例如芯片标识符。 -We use the function `read.celfiles` from the `r Biocpkg("oligo") ` -package [@oligo] to import the files: +我们使用 `r Biocpkg("oligo") ` 包[@oligo]中的 `read.celfiles` +函数导入文件： ```{r importCelfiles, results="hide", eval=TRUE, dependson="getSDRF", warning = FALSE } raw_data <- oligo::read.celfiles(filenames = file.path(raw_data_dir, @@ -399,41 +399,41 @@ raw_data <- oligo::read.celfiles(filenames = file.path(raw_data_dir, stopifnot(validObject(raw_data)) ``` -This automatically creates an ExpressionSet, -fills the sections "array data" with the data from the CEL files -and uses the correct chip annotation package, in this case -`r Biocannopkg("pd.hugene.1.0.st.v1") ` (the -chip-type is also stored in the .CEL files). +这会自动创建一个 ExpressionSet, +用 CEL 文件中的数据填充部分 "array data" +并使用正确的芯片注释包，在本例中为 +`r Biocannopkg("pd.hugene.1.0.st.v1") ` ( +芯片类型也存储在 .CEL 的文件)。 -Furthermore, we specified our `AnnotatedDataFrame` "`SDRF`" -created earlier from the SDRF file as `phenoData`. Thus, -we had to make sure to import the CEL files -in the order that corresponds to the SDRF table --- to enforce this, -we used the column `Array.Data.File` of the `SDRF` table as the `filenames` -argument. +此外，我们将早先从 SDRF 文件创建的 `AnnotatedDataFrame` + "`SDRF`" 指定为 `phenoData`. 因此， +我们必须确保按照与 SDRF 表对应的顺序导入 CEL 文件 +- 为了强制执行此操作，我们使用 `SDRF` +表的 `Array.Data.File` 列作为 `filenames` +(文件名) 参数。 -Finally, we checked whether the object created is valid (e.g. sample names -match between the different tables). +最后，我们检查了创建的对象是否有效 +(例如，不同表之间的样本名称是否匹配）。 -We now have a first look on the raw data. +我们现在首先看一下原始数据。 -The `pData` function of the Biobase package directly accesses the phenoData in the ExpressionSet `raw_data`. With the `head()` function, we can view the -first six lines of the table. We have a look at the columns included -and retain only those columns that are related to the experimental -factors of interest. +Biobase 包的 `pData` 函数直接访问 ExpressionSet `raw_data` 中的 phenoData。使用 `head()` 函数，我们可以查看 +表的前六行。我们需要查看包含的列， +并仅保留与感兴趣的实验因子 +相关的列。 ```{r inspectPhenoData, eval=TRUE } head(Biobase::pData(raw_data)) ``` -The columns of interest for us are the following: +我们感兴趣的列如下： - * identifiers of the individuals, i.e. columns "Source.Name", "Characteristics.individual." - * disease of the individual, i.e. "Factor.Value.disease." - * mucosa type, i.e. "Factor.Value.phenotype." + * 个人的标识符，即列名 "Source.Name", "Characteristics.individual." + * 个体的疾病，即名为 "Factor.Value.disease." + * 粘膜类型，即名为 "Factor.Value.phenotype." -We now subselect the corresponding columns: +我们现在再次选择相应的列： ```{r reassignpData, eval=TRUE} Biobase::pData(raw_data) <- Biobase::pData(raw_data)[, c("Source.Name", @@ -442,29 +442,29 @@ Biobase::pData(raw_data) <- Biobase::pData(raw_data)[, c("Source.Name", "Factor.Value.phenotype.")] ``` -# Quality control of the raw data +# 原始数据的质量控制 -The first step after the initial data import is the quality control of the data. -Here we check for outliers and try to see whether the data clusters as expected, -e.g. by the experimental conditions. The expression intensity values are in the -assayData sub-object "exprs" and can be accessed by the `exprs(raw_data)` -function. The rows represent the microarray probes, i.e. -the single DNA locations on the chip, while the columns -represent one microarray, i.e. a sample of inflamed -and non-inflamed tissue of every patient, respectively. +初始数据导入后的第一步是数据的质量控制。 +在这里，我们检查异常值并尝试查看数据是否按预期聚类， +例如 由实验条件。 表达强度值在 +assayData 子对象 "exprs" 中，并且可以由 `exprs(raw_data)` +函数访问。 每行代表微阵列探针，即 +芯片上的单个 DNA 位置， +而每列代表一个微阵列， +即每个患者的发炎和非发炎组织的样品。 ```{r inspectAssayData, eval=TRUE} Biobase::exprs(raw_data)[1:5, 1:5] ``` -For quality control, we take the log2 of `Biobase::exprs(raw_data)`, -as expression data is commonly analyzed on a logarithmic scale. +对于质量控制，我们采用 `Biobase::exprs(raw_data) `的 log2， +因为表达数据通常以 log 分析。 -We then perform a principal component analysis (PCA) and plot it (Figure -\@ref(fig:qualityControlRawDataPCA)). Every point in the plot -represents one sample, with the colour indicating the mucosa type -(inflamed vs non-inflamed) and the shape indicating the disease (UC or CD). +然后，我们执行主成分分析（PCA）并绘制它 \@ref(fig:qualityControlRawDataPCA)) +图中的每个点代表一 +个样本，颜色表示粘膜类型 +(发炎与未发炎）和表示疾病的形状（ UC 或 CD ）。 @@ -491,36 +491,36 @@ ggplot(dataGG, aes(PC1, PC2)) + scale_color_manual(values = c("darkorange2", "dodgerblue4")) ``` -The PCA plot (Figure \@ref(fig:qualityControlRawDataPCA), performed on -the log-intensity scale) of the raw data shows that the -first principal component differentiates between the diseases. This means that -the disease type is a major driver of gene expression differences. -This might hinder our analysis, as we want to analyze the differential expression between inflamed and non-inflamed -tissues, independently of the disease a person suffers from. +原始数据的 PCA 图 (Figure \@ref(fig:qualityControlRawDataPCA), 在对数强度标度上进行 +) 显示 第一主成分 (first PC) 区分疾病。 +这意味着 +疾病类型是基因表达差异的主要驱动因素。 +这可能会妨碍我们的分析，因为我们想要分析发炎组织和非发炎组织之间的差异表达， +而不是人们患病的疾病的关系。 -We also represent the probe intensities via a boxplot graph with one box per -individual microarray. (Figure \@ref(fig:qualityControlRawDataBox)). Note that the -`oligo::boxplot` function, i.e. the boxplot function of the oligo package, can take -expression sets as argument. It accesses the expression data and performs a -log2-transformation by default. We therefore can use `raw_data` as argument here. +我们还通过箱形图图表示探针强度， +每个单独的微阵列是有一个箱。(Figure \@ref(fig:qualityControlRawDataBox)). 注意， +`oligo::boxplot` 函数，即 oligo 包的 boxplot() 函数， +可以将表达式集作为参数。 +它访问表达式数据并默认执行 log2 转换。 因此，我们可以在这里使用 `raw_data` 作为参数。 ```{r qualityControlRawDataBox, fig.cap="Intensity boxplots of the log2–transformed raw data."} oligo::boxplot(raw_data, target = "core", main = "Boxplot of log2-intensitites for the raw data") ``` -When looking at the boxplot (Figure -\@ref(fig:qualityControlRawDataBox)), we see that the intensity distributions of the -individual arrays are quite different, indicating the need for an appropriate -normalization, which we will discuss next. +当查看箱形图 (Figure \@ref(fig:qualityControlRawDataBox)) 时， +我们看到各个阵列的强度分布完全不同， +这表明需要进行适当的归一化， +我们将在下面讨论。 -Until now, we have only performed a very basic quality control; -more elaborate quality control plots are available in the package -`r Biocpkg("arrayQualityMetrics")` [@AQM]. The package produces an html report, -containing the quality control plots together with a description of their -aims and an identification of possible outliers. We do not discuss this tool in -detail here, but simply provide the code below, -which creates a report for our raw data. +到目前为止，我们只进行了非常基本的质量控制; +`r Biocpkg("arrayQualityMetrics")` [@AQM] 中提供了 +更精细的质量控制图。 该软件包可以生成一个 html 报告， +其中包含质量控制图以及对其目标的描述以及可能的异常值的识别。 +我们在此不详细讨论此工具， +只是提供下面的代码， +它为我们的原始数据创建报告。 ```{r arrayQualityMetricsRaw, eval = TRUE, warning=FALSE, message=FALSE} arrayQualityMetrics(expressionset = raw_data, @@ -529,73 +529,73 @@ arrayQualityMetrics(expressionset = raw_data, intgroup = c("Factor.Value.disease.", "Factor.Value.phenotype.")) ``` -# Background adjustment, calibration, summarization and annotation +# 背景调整，校准，摘要和注释 -## Background adjustment +## 背景调整 -After the initial import and quality assessment, the next step in processing of -microarray data is background adjustment. This is essential because -a proportion of the measured probe intensities are due to non-specific -hybridization and the noise in the optical detection system. -Therefore, observed intensities need to be adjusted to -give accurate measurements of specific hybridization. +在初始导入和质量评估之后， +处理微阵列数据的下一步是背景调整。 +这是必要的因为测量的探针强度的一部分 +是由于非特异性杂交和光学检测系统中的噪声。 +因此，需要调整观察到的强度 +以给出特定杂交的准确测量。 -## Across-array normalization (calibration) +## 跨阵列归一化（校准） -Normalization across arrays is needed in order to be able -to compare measurements from different array hybridizations -due to many obscuring sources of variation. -These include different efficiencies of -reverse transcription, labeling or hybridization -reactions, physical problems with the arrays, reagent batch effects, -and laboratory conditions. +因为许多差异的来源并不明确， +所以要标准化跨阵列以便能够比较来自不同阵列杂交的测量结果。 +这些包括逆转录， +标记或杂交反应的不同效率， +阵列的物理问题， +试剂批次效应 +和实验室条件。 -## Summarization +## 概要 -After normalization, summarization is necessary to be done because on -the Affymetrix platform, transcripts are represented -by multiple probes, that is multiple locations on the array. -For each gene, the background-adjusted and normalized intensities of all probes -need to be summarized into one quantity that estimates an amount proportional to -the amount of RNA transcript. +我们必须在归一化后进行汇总， +因为在Affymetrix平台上， +转录本由多个探针表示，即阵列上的多个位置。 +对于每个基因，需要将所有探针的背景调整 +和标准化强度总结为一个量， +让预计与 RNA 转录物的量成比例的量。 -After the summarization step, the summarized data can be annotated with various -information, e.g. gene symbols and ENSEMBL gene identifiers. There is an -annotation database available from Bioconductor -for our platform, namely the package -`r Biocannopkg("hugene10sttranscriptcluster.db") `. +在汇总步骤之后，汇总数据可以用各种信息注释， +例如，基因符号和 ENSEMBL 基因标识符。 +Bioconductor 为我们的平台提供了一个注释数据库， +即 +`r Biocannopkg("hugene10sttranscriptcluster.db") `包。 -You can view its content like this: +您可以像这样查看其内容： ```{r annotationDataBaseContent, eval = TRUE} head(ls("package:hugene10sttranscriptcluster.db")) ``` -Additional information is available from the reference manual of the package. -Essentially, the package provides a mapping from the transcript cluster -identifiers to the various annotation data. +更多信息可从包装的参考手册中获得。 +本质上该包提供从转录本簇标识 +符到各种注释数据的映射。 -## Old and new "probesets" of Affymetrix microarrays +## Affymetrix 微阵列的新旧 “探针组” -Traditionally, Affymetrix arrays (the so-called 3' IVT arrays) -were probeset based: a certain fixed group of probes were part of a probeset -which represented a certain gene or transcript (note however, that a -gene can be represented by multiple probesets). +传统来说，Affymetrix 阵列（所谓的3'IVT 阵列）基于探针组： +某一固定探针组是探针组的一部分， +该探针组代表某种基因或转录物 +但请注意，基因可由多个探针组表示）。 -The more recent "Gene" and "Exon" Affymetrix arrays are exon based and hence -there are two levels of summarization to get to the gene level. The "probeset" summarization leads to the exon level. The gene / transcript level is given by "transcript clusters". Hence, the appropriate annotation package for -our chip type is called +最近的“基因”和“外显子” Affymetrix 阵列是基于外显子的， +因此有两个级别的总结可以达到基因水平。探针组”总结导致外显子水平。基因/转录物水平由“转录物簇”给出。 因此，我们的芯片类型的适当注释包 +名为 `r Biocpkg("hugene10sttranscriptcluster.db") `. -"Gene" arrays were created as affordable versions of the "Exon" arrays, by only -taking the "good" probes from the Exon array. Initially on the Exon array, -at least four probes were part of one "Exon". -With the thinned out "Gene" array, many probesets were made up of three or fewer -probes. This is visualized in Figure -\@ref(fig:DifferenceBetweenExonAndGeneTypeArrays): Single probesets -are indicated by single colours; probesets representing one gene -are indicated by a colour shade: e.g., all yellow probes belong to one Exon, -and all yellow, orange and red probesets belong to one gene: +通过仅从外显子阵列中获取“好”的探针， +“基因”阵列被创建为“外显子”阵列的负担得起的版本。 +最初在外显子阵列上，至少有四个探针是一个“外显子”的一部分。 +通过稀释的“基因”阵列，许多探针组由三个或更少的探针组成。 +这在图4 中可视化：\@ref(fig:DifferenceBetweenExonAndGeneTypeArrays) +单个探针组由单色标注; +表示一个基因的探针组用色调表示： +例如，所有黄色探针属于一个外显子， +并且所有黄色，橙色和红色探针组属于一个基因： ```{r DifferenceBetweenExonAndGeneTypeArrays, fig.width=7, fig.height=6, echo=FALSE, fig.cap="Visualization of the difference between \"Exon\" type array (left) and \"Gene\" type array (right)."} library(grid) @@ -634,64 +634,64 @@ grid.rect(x=dat$x,y=dat$y, height=0.1, width=0.1, hjust=0,vjust=0,vp=vp2, ``` -On the left side, we see plenty of probes for each Exon / probeset -(i.e. each colour): therefore, a summarization on the probeset / exon level -makes sense. In the gene type array, however, only a small proportion of the -original probes per probeset is included. Thus, a summarization on the -probeset / exon level is not recommended for "Gene" arrays -but nonetheless possible by using the `r Biocannopkg("hugene10stprobeset.db") ` -annotation package. +在左侧，我们看到每个外显子/探针组（即每种颜色）有大量探针： +因此，探针组/外显子水平上的摘要是合理的。 +然而，在基因型阵列中， +每个探针组仅包括一小部分原始探针。 +因此，不建议对 “基因” 阵列使用探针组/外显子级别进行汇总， +但仍然可以使用 `r Biocannopkg("hugene10stprobeset.db") ` +注释包。 -Note that furthermore, there are also no longer designated match/mismatch probes -present on "Gene" and "Exon" type chips. The mismatch probe was initially -intended as base-level for background correction, but hasn't prevailed -due to more elaborate background correction techniques that -do not require a mismatch probe. +注意，此外，“基因”和“外显子”型芯片上 +也不再存在指定的匹配/错配探针。 +错配探针最初用作背景校正的基础水平， +但由于更精细的背景校正技术 +不需要不匹配探针而没有普及。 -## One-step preprocessing in oligo +## 一步预处理寡核苷酸 -The package `r Biocpkg("oligo") ` allows us to perform background correction, -normalization and summarization in one single step using a deconvolution -method for background correction, quantile normalization and -the RMA (robust multichip average) algorithm for summarization. +`r Biocpkg("oligo") ` 包允许我们在一个步骤里使用用于背景校正， +归一化，和汇总。 +这个包使用反卷积方法来背景校正，分位数归一化和 +RMA (robust multichip average) 算法来汇总。 -This series of steps as a whole is commonly referred to as RMA algorithm, -although strictly speaking RMA is merely a summarization method +整个系列步骤通常称为 RMA 算法， +但严格来说 RMA 仅仅是一种汇总的方法 [@Irizarry_2003; @Bolstad_2003; @Irizarry_2003a]. -# Relative Log Expression data quality analysis +# 相对日志表达式数据质量分析 -Before calibrating and evaluating the data, we want to perform another quality -control procedure, namely Relative Log Expression (RLE), as described in -the article by Gandolfo et al [@Gandolfo_2018]. To this end, we first -perform an RMA without prior normalization: +在校准和评估数据之前，我们要执行另一个质量控制程序， +即相对对数表达式（RLE）， +如 Gandolfo et al [@Gandolfo_2018]。 +我们首先执行没有事先规范化 RMA： ```{r RMAcalibrationForRLE, eval=TRUE} palmieri_eset <- oligo::rma(raw_data, target = "core", normalize = FALSE) ``` -Further details on the RMA algorithm will be provided after RLE analysis, when -the "full" RMA is carried out, including normalization. +在 RLE 分析之后，或当执行 “完整” RMA（包括归一化）时， +将提供关于 RMA 算法的进一步细节。 -## Computing the RLE +## 计算 RLE -The RLE is performed by calculating the median log2 intensity of every -transcript across all arrays. +通过计算所有阵列上每个转录物的 log2 +强度中值来进行 RLE。 -We do this by calculating the row medians of `exprs(palmieri_eset)`, -as the transcripts are represented by the rows and the single microarrays -by the columns. +我们通过计算 `exprs(palmieri_eset)`的行中值来做到这一点， +因为转录本由每行表示， +单个微阵列由每列表示。 -Note that we do not have to apply the log2 manually, -as the output data of the RMA function is in log2 scale by default. +注意，我们不必手动应用 log2， +因为默认情况下 RMA 函数的输出数据为 log2 比例。 -We then substract this transcript median intensity from every -transcript intensity via the `sweep` function. +然后我们通过 `sweep`函数从每个转录物强度中 +减去该转录物中值强度。 -## Plotting the RLE +## 绘制 RLE 图 -We finally reshape the data into a format in which we can use to create a boxplot -for each array, as before: +我们最终将数据重新整形为一种格式，我们可以使用这种格式 +为每个数组创建一个箱形图 (boxplot) ，如前所述： ```{r boxplotDataForRLETidy, fig.cap="Boxplot for the RLE values", warning=FALSE} @@ -713,34 +713,34 @@ ggplot2::ggplot(RLE_data_gathered, aes(patient_array, face = "bold")) ``` -Note that the y-axis now displays for each microarray the deviation -of expression intensity from the median expression of the respective single -transcripts across arrays. +注意，现在 y 轴显示的是针对每个微阵列表达强度 +与跨阵列的相应单个转录物 +的表达中值的偏差。 -Boxes with a larger extension therefore indicate an unusually high deviation -from the median in a lot of transcripts, suggesting that these arrays -are different from most of the others in some way. +因此，具有较大延伸的盒子表示在许多转录物中 +与中值的异常高偏差， +表明这些阵列在某些方面与大多数其他阵列不同。 -Boxes that are shifted in y-direction indicate a systematically higher or lower -expression of the majority of transcripts in comparison to most of the other -arrays. This could be caused by quality issues or batch effects. +与大多数其他阵列相比，在 y 方向上移位的框表示 +大多数转录物的系统性更高或更低的表达。 +这可能是由质量问题或批量效应引起的。 -Therefore, if shape and median of a given box varies too much from the -bulk, they should be inspected and potentially removed. +因此，如果给定盒子的形状和中位数与体积相差太大， +则应检查并可能将其移除。 -By inspecting the boxplot in Figure \@ref(fig:boxplotDataForRLETidy), -five arrays could be considered as outliers: -2826_I, 2826_II, 3262_II, 3302_II and 3332_II are negatively y-shifted. +通过检查图 5 中的箱线图 \@ref(fig:boxplotDataForRLETidy), +可以将五个阵列视为异常值: +2826_I，2826_II，3262_II，3302_II 和 3332_II 被负 y 移位。 -We will keep these samples in mind for heatmap cluster analysis later on in -the workflow. Arrays that are confirmed to be outliers by heatmap -analysis could be removed for subsequent analysis. +我们将记住这些样本， +以便稍后在工作流程中进行热图聚类分析。 +通过热图分析确认为异常值的阵列可以在后续分析前移除。 -# RMA calibration of the data +# RMA 校准数据 -Now, we can apply the full RMA algorithm to our data in order to -background-correct, normalize and summarize: +现在，我们可以将完整的 RMA 算法应用于我们的数据， +以便进行背景校正，规范化和总结： ```{r RMAcalibrationWITHnormalization, eval=TRUE} palmieri_eset_norm <- oligo::rma(raw_data, target = "core") From 23232235c6780cbe63801324b224c8ad974642da Mon Sep 17 00:00:00 2001 From: Xiao W Date: Mon, 1 Jul 2019 12:12:58 -0400 Subject: [PATCH 03/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 208 ++++++++++++------------- 1 file changed, 104 insertions(+), 104 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index db6ec4f..c763491 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -746,63 +746,63 @@ ggplot2::ggplot(RLE_data_gathered, aes(patient_array, palmieri_eset_norm <- oligo::rma(raw_data, target = "core") ``` -The parameter `target` defines the degree of summarization, the -default option of which is "core", using transcript clusters containing -"safely" annotated genes. For summarization on the exon level (not recommended for Gene -arrays), one can use "probeset" as the target option. -Although other methods for background correction and normalization exist, -RMA is usually a good default choice. -RMA shares information across arrays and -uses the versatile quantile normalization method that -will make the array intensity distributions match. However, it is preferable -to apply it only after outliers have been removed. -The quantile normalization algorithm used by RMA -works by replacing values by the average of identically -ranked (within a single chip) values across arrays. A more detailed -description can be found on the [Wikipedia page](https://en.wikipedia.org/wiki/Quantile_normalization). - -An alternative to quantile normalization would be the `r Biocpkg("vsn") ` -algorithm,that performs background correction and normalization by robustly -shifting and scaling intensity values within arrays before log-transforming -them. This is less "severe" than quantile normalization [@vsn]. - -## Some mathematical background on normalization (calibration) and background correction - -A generic model for the value of the intensity \(Y\) of a single probe on -a microarray is given by +参数 `target` 用来定义摘要的程度， +默认选项是“核心”， +使用含有“安全”注释基因的转录物簇。对于外显子级别的汇总（不建议用于基因阵列)， +可以使用 "probeset" 作为 target 选项。 +虽然存在其他背景校正和归一化方法， +但 RMA 通常是一个很好的默认选择。 +RMA 跨阵列共享信息， +并使用通用的分位数归一化方法， +使阵列强度分布匹配。但是，最好只在 +删除异常值后才应用它。 +RMA 使用的分位数归一化算法 +通过将数值替换为数组中相同排名（单个芯片内）值的平均值来工作。 +更详细的描述可以在 +维基百科页面上找到 [Wikipedia page](https://en.wikipedia.org/wiki/Quantile_normalization). + +分位数归一化的替代方案是 `r Biocpkg("vsn") ` 算法， +其通过在对数变换之前对阵列内的强度值 +进行鲁棒地移位和缩放来执行背景校正和归一化。 +这跟分位数归一化 [@vsn] 比不那么“严重”。 + +## 关于归一化（校准）和背景校正的一些数学背景知识 + +微阵列上单个探针的强度 \(Y\) 值的通用模型 +由下式给出： \[ Y = B + \alpha \cdot S \] -where B is a random quantity due to background noise, usually composed of -optical effects and non-specific binding, \(\alpha\) is a gain factor, -and \(S\) is the amount of measured specific binding. The signal \(S\) -is considered a random variable as well and accounts for measurement -error and probe effects. -The measurement error is typically assumed to be multiplicative so we can write: +其中 B 是由于背景噪音引起的随机数量， +通常由光学效应和非特异性结合组成， \(\alpha\) 是增益因子， + \(S\) 是测量的特异性结合量。 信号 \(S\) +也被认为是随机变量， +并且考虑了测量误差和探测效应。 +通常假设测量误差是乘法的，因此我们可以写： \[ \log(S) = \theta + \varphi + \varepsilon \] -Here \(\theta\) represents the logarithm of the true abundance, -\(\varphi\) is a probe-specific effect, and $\varepsilon$ accounts for the -nonspecific error. -This is the additive-multiplicative-error model for microarray data used -by RMA and also the `r Biocpkg("vsn") ` algorithm [@vsn]. The algorithms differ -in the way that \(B\) is removed and an estimate of \(\theta\) is obtained. +这里 \(\theta\) 表示真实丰度的对数， +\(\varphi\) 是探针特异性效应 $\varepsilon$ 表示 +非特异性误差。 +这是 RMA 使用的微阵列数据的加法 - 乘法 - 误差模型 +以及 `r Biocpkg("vsn") ` 算法 [@vsn]. 算法的不同之处 +在于去除 \(B\) 并获得 \(\theta\) 的估计。 -## Quality assessment of the calibrated data +## 校准数据的质量评估 -We now produce a clustering heatmap and another PCA plot using the -calibrated data. +我们现在使用校准数据生成聚类热图 +和另一个 PCA 图。 -**PCA analysis** +**PCA 分析** -First, we perform a PCA analysis of the calibrated data analogously -to the one with the raw data: +首先，我们对校准数据进行 PCA 分析， +类似于原始数据： ```{r PCAMetricsCalibrated, fig.cap = "PCA plot of the calibrated, summarized data.", eval = TRUE } exp_palmieri <- Biobase::exprs(palmieri_eset_norm) @@ -829,18 +829,18 @@ ggplot(dataGG, aes(PC1, PC2)) + scale_color_manual(values = c("darkorange2", "dodgerblue4")) ``` -In comparison to the first PCA analysis before RMA (Figure \@ref(fig:qualityControlRawDataPCA)), we see that now the first principal -component separates between the tissues types (Figure \@ref(fig:PCAMetricsCalibrated)). -This indicates that now differential expression between the tissue types is the -dominant source of variation. Note that the second principal component -separates the diseases. +与 RMA 之前的第一次 PCA 分析相比 (Figure \@ref(fig:qualityControlRawDataPCA)), 我们看到现在第一个主要成分 PC +在组织类型之间分离 (Figure \@ref(fig:PCAMetricsCalibrated))。 +这表明现在组织类型之间的差异表达 +是变异的主要来源。 注意，第二 PC +分离疾病。 -**Heatmap clustering analysis** +**热图聚类分析** -We want to plot a heatmap with the sample-to-sample distances with the -sample names as row-names. Also, we want to see how well the samples -cluster for phenotype (inflamed and non-inflamed tissue) and disease -(UC and CD), respectively. We use row annotation for that: it means that these features get a colour code and will be displayed to the left of each row. +我们想要绘制具有样本到样本距离的热图， +样本名称为行名。 此外，我们希望看到样本分别表现为表型 +（发炎和非发炎组织）和疾病（UC和CD）的程度。 +我们使用行注释：这意味着这些功能获得颜色代码并将显示在每行的左侧。 ```{r rownamesForHeatmap, fig.height = 8.5, fig.width = 7, eval = TRUE, echo=TRUE} phenotype_names <- ifelse(str_detect(pData @@ -858,17 +858,17 @@ row.names(annotation_for_heatmap) <- row.names(pData(palmieri_eset_norm)) ``` -In order to map the sample-to-sample distances, we first compute the distances -using the `dist` function. We need to transpose the expression values since -the function computes the distances between the rows (i.e. genes in our case) by -default. The default distance is the Euclidean one. However this can be -changed and we choose the Manhattan distance here (it uses absolute -distances along rectangular paths instead of squared distances -of the direct path), as it is more robust. -We set the diagonal of the distance matrix to ` NA ` in order -to increase the contrast of the color coding. Those diagonal entries -do not contain information since the distance of a sample -to itself is always equal to zero. +为了映射样本到样本的距离， +我们首先使用 `dist` 函数计算距离。我们需要转换表达式值， +因为函数默认计算行之间的距离（例如我们的基因)。 +默认距离是欧几里德距离。然而，这可以改变， +我们在这里选择曼哈顿距离 +（它使用沿矩形路径的绝对距离而不是直接路径的平方距离）， +因为它更稳健。 +我们将距离矩阵的对角线设置为 ` NA ` ， +以增加颜色编码的对比度。 +那些对角线条目不包含信息， +因为样本到自身的距离总是等于零。 ```{r HeatmapWithAnnotation, fig.height = 8.5, fig.width = 7, fig.cap="Heatmap of the sample-to-sample distances"} @@ -893,39 +893,39 @@ pheatmap(dists, col = (hmcol), legend_labels = (c("small distance", "large distance")), main = "Clustering heatmap for the calibrated samples") ``` -On the heatmap plot (Figure \@ref(fig:HeatmapWithAnnotation)) we also see that -the samples do not cluster strongly by tissue, confirming the impression -from the PCA plot (Figure \@ref(fig:PCAMetricsCalibrated)) that the separation -between the tissues is not perfect. The yellow stripes in the heatmap might correspond -to outliers that could potentially be removed: the ones that could be flagged -here are 2826_II, 3262_II, 3271_I, 2978_II and 3332_II. -2826_II, 3262_II and 3332_II were found to be outliers in both RLE and heatmap -analysis and could therefore potentially be removed; however, in order -to stay as close as possible to the -original paper, we continue with the complete set of samples. -Note again that more elaborate metrics to identify and remove outliers -are provided by the `r Biocpkg("arrayQualityMetrics")` package. - - -# Filtering based on intensity - -We now filter out lowly expressed genes. Microarray data commonly show -a large numberof probes in the background intensity range. These probes -also do not change much across arrays. Hence they combine a low variance -with a low intensity. Thus, they could end up being detected as differentially -expressed although they are barely above the "detection" limit and are not very -informative in general. - -We will perform a "soft" intensity based filtering here, since this -is recommended by the `r Biocpkg("limma")` [@limma; @Smyth_2004] user guide -(a package we will use below for the differential -expression analysis). - -However, note that a variance based filter might exclude a similar set of -probes in practice. -For intensity-based filtering, we calculate the row-wise medians from the -expression data, as they represent the transcript medians, and assign them to -`palmieri_medians`. From this we create a histogram: +在热图 (Figure \@ref(fig:HeatmapWithAnnotation)) 中， +我们还发现样品不会被组织强烈聚集， +证实了 PCA 图 (Figure \@ref(fig:PCAMetricsCalibrated)) 的印象， +即组织之间的分离并不完美。热图中的黄色条纹 +可能对应于可能被删除的异常值：这里可以标记的异常值 +是 2826_II, 3262_II, 3271_I, 2978_II 和 3332_II。 +因为发现 2826_II, 3262_II 和 3332_II 在 RLE 和热图分析中都是异常值， +因此可能被删除; 但是， +为了尽可能接近原始论文， +我们继续提供完整的样本集。 +注意 `r Biocpkg("arrayQualityMetrics")`包提供了更精细的标识 +来识别和删除异常值。 + + +# 根据强度过滤 + +我们现在过滤掉低表达的基因。微阵列数据通常 +在背景强度范围内显示大量探针。 这些探针 +在数组之间也没有太大变化。因此， +它们结合了低方差和低强度。因此，它们最终可能被检测为差异表达， +尽管它们几乎不高于 “detection” 限制并且 +通常不是非常有用的信息。 + +我们将在这里执行基于“软”强度的过滤，因为这是 + `r Biocpkg("limma")` [@limma; @Smyth_2004] 用户指南推荐的 +(我们将在下面用于差异表达式分析的包 +差异表达式分析的包)。 + +但请注意，基于方差的过滤器在实践中可能会排除 +一组类似的探测。 +对于基于强度的过滤，我们从表达式数据计算行式中位数， +因为它们代表转录中位数，并将它们分配给 +`palmieri_medians`。从这里我们创建一个直方图： ```{r IntensityBasedManualFiltering, fig.cap="Histogram of the median intensities per gene"} palmieri_medians <- rowMedians(Biobase::exprs(palmieri_eset_norm)) @@ -936,14 +936,14 @@ hist_res <- hist(palmieri_medians, 100, col = "cornsilk1", freq = FALSE, xlab = "Median intensities") ``` -In the histogram of the gene-wise medians (Figure -\@ref(fig:IntensityBasedManualFiltering)), we -can clearly see an enrichment of low medians on the left hand side. -These represent the genes we want to filter. In order to infer a cutoff from the -data, we inspect the histogram: We visually set a cutoff line `man_threshold` -to the left of the histogram peak in order not to exclude too many genes. In our -example, we choose a threshold of 4. We plot the same histogram as before -and add the threshold line with the `abline()` function: +在基因中位数的直方图中 (Figure +\@ref(fig:IntensityBasedManualFiltering))， +我们可以清楚地看到左侧低中值的富集。 +这些代表了我们想要过滤的基因。 为了推断数据的截止， +我们检查直方图：我们在直方图峰的左侧可视地设置截止线 `man_threshold` +以便不排除太多基因。 在我们的示例中， +我们选择阈值 4 。我们绘制与之前相同的直方图， +并使用 `abline()` 函数添加阈值线： ```{r setManualThreshold, fig.cap="Histogram of the median intensities per gene with manual intensity filtering threshold (red line)."} man_threshold <- 4 From ae0ca6552f6e8954432b98c1bf5e064aef0c9703 Mon Sep 17 00:00:00 2001 From: Xiao W Date: Mon, 1 Jul 2019 19:42:46 -0400 Subject: [PATCH 04/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 362 ++++++++++++------------- 1 file changed, 181 insertions(+), 181 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index c763491..df06fd9 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -956,11 +956,11 @@ hist_res <- hist(palmieri_medians, 100, col = "cornsilk", freq = FALSE, abline(v = man_threshold, col = "coral4", lwd = 2) ``` -Transcripts that do not have intensities larger than the threshold in -at least as many arrays as the smallest experimental group are excluded. +排除在至少与最小实验组一样多的阵列中 +不具有大于阈值的强度的转录物。 -In order to do so, we first have to get a list with the number -of samples (=arrays) (`no_of_samples`) in the experimental groups: +为了做到这一点，我们首先必须得到一个列表， +其中包含实验组中的样本数（=数组）(`no_of_samples`): ```{r expGroups, dependson="PCAMetricsCalibrated"} no_of_samples <- @@ -970,21 +970,21 @@ no_of_samples ``` -We now filter out all transcripts -that do not have intensities greater than the threshold in at least -as many arrays as the smallest experimental group (`r min(no_of_samples)`) -which we define as `samples_cutoff`. +我们现在过滤掉所有 +不具有大于阈值强度的转录本， +至少与我们定义为 (`r min(no_of_samples)`) + `samples_cutoff` 一样多。 -A function `idx_man_threshold` is then applied to each row, i.e. -to each transcript across all arrays. -It evaluates whether the number of arrays where the median intensity passes the -threshold (`sum(x > man_threshold)`) is greater than the `samples_cutoff` and -returns TRUE or FALSE for each row, i.e. each transcript. +然后将函数 `idx_man_threshold` 应用于每一行， +即应用于所有数组的每个转录本。 +它评估中值强度超过阈值的数组的数量 +(`sum(x > man_threshold)`) 是否大于 `samples_cutoff` 并且 +对于每一行（即每个转录本）返回 TRUE 或 FALSE。 -We then create a table of `idx_man_threshold` to summarize the results and -get an overview over how many genes are filtered out. -In the last step, we subset our expression set to `palmieri_manfiltered` -and keep the TRUE elements of `idx_man_threshold`. +然后，我们创建一个 `idx_man_threshold` 表来汇总结果， +并概述过滤掉多少基因。 +在最后一步中，我们将表达式集合为 `palmieri_manfiltered`， +并保留 `idx_man_threshold`的 TRUE 元素。 ```{r filteringOfLowIntensity_transcripts} samples_cutoff <- min(no_of_samples) @@ -998,12 +998,12 @@ palmieri_manfiltered <- subset(palmieri_eset_norm, idx_man_threshold) ``` -# Annotation of the transcript clusters +# 抄本集群的注释 -Before we continue with the linear models for microarrays and differential -expression, we first add "feature data", i.e. annotation -information to the transcript cluster identifiers stored in the featureData of -our ExpressionSet: +在继续使用微阵列和差异表达的线性模型之前， +我们首先将“特征数据”， +即注释信息添加到存储在 ExpressionSet 的 featureData +中的转录集群标识符： ```{r annotateData, eval=TRUE, dependson="intensityBasedFiltering", message = FALSE} anno_palmieri <- AnnotationDbi::select(hugene10sttranscriptcluster.db, @@ -1014,20 +1014,20 @@ anno_palmieri <- AnnotationDbi::select(hugene10sttranscriptcluster.db, anno_palmieri <- subset(anno_palmieri, !is.na(SYMBOL)) ``` -We used the function `select` from `r Biocpkg("AnnotationDbi") ` to query -the gene symbols and associated short descriptions for the transcript clusters. -For each cluster, we added the gene symbol (`SYMBOL`) and a short description -of the gene the cluster represents (`GENENAME`). +我们使用 `r Biocpkg("AnnotationDbi") ` 中的 `select` 来查询 +转录本簇的基因符号和相关的简短描述。 +对于每个簇，我们添加了基因符号 (`SYMBOL`) 和簇代表的基因 +的简短描述 (`GENENAME`)。 -In a second step, we filtered out the probes that do not map to a gene, -i.e. that do not have a gene symbol assigned. +在第二步中，我们过滤掉了未映射到基因的探针， +即没有分配基因符号的探针。 -## Removing multiple mappings +## 删除多个映射 -Many transcript-cluster identifiers will map to multiple gene symbols, i.e. -they can't be unambigously assigned. +许多转录物簇标识符将映射到多个基因符号， +即它们不能被明确地分配。 -We compute a summary table in the code below to see how many there are: +我们在下面的代码中计算一个汇总表，看看有多少： ```{r multipleMappings, dependson="annotateData"} @@ -1046,25 +1046,25 @@ probe_stats <- anno_filtered nrow(probe_stats) ``` -First, we grouped `anno_palmieri` by their PROBEID; that way, -the subsequent operations are not carried through for each single row, -but for each group, i.e. each PROBEID. -We then summarized the groups and indicate the number of different genes -assigned to a transcript cluster in the column `no_of_matches`. -Finally, we filtered for PROBEIDs with multiple matches, -i.e. `no_of_matches > 1`. +首先，我们用他们的 PROBEID 对 `anno_palmieri` 进行了分组; +这样，后续操作不会针对每一行进行， +而是针对每个组，即每个 PROBEID 。 +然后，我们对这些组进行了总结， +并在 `no_of_matches` 列中指出了分配给转录物簇的不同基因的数量。 +最后，我们使用多个匹配过滤了 PROBEIDs ， +即 `no_of_matches > 1`。 -With `dim(probe_stats)`, we could see how many -probes have been mapped to multiple genes. +使用 `dim(probe_stats)`， 我们可以看到 +有多少探针被映射到多个基因。 -We have close to 2000 transcript clusters that map to multiple gene symbols. -It is difficult to decide which mapping is "correct". Therefore, -we exclude these transcript clusters. +我们有接近 2000 个转录物簇，映射到多个基因符号。 +很难确定哪个映射是“正确的”。 +因此，我们排除这些转录本群集。 -We want to remove those probe IDs that match the ones in `probe_stats`, -as those are the probes with multiple mappings. We assign these IDs -to the variable `ids_to_exclude`. Then, we generate `palmieri_final`, -an expression set without the `ids_to_exclude`. +我们想要删除那些与 `probe_stats` 中的 probe IDs 匹配的 ID。 +因为它们是具有多个映射的探测器。我们将这些 ID +分配给变量 `ids_to_exclude`。然后，我们生成 `palmieri_final`, +一个没有 `ids_to_exclude`的表达式集。 ```{r excludeMultipleMappingsFromAssayData, dependson="multipleMappings"} ids_to_exlude <- (featureNames(palmieri_manfiltered) %in% probe_stats$PROBEID) @@ -1076,29 +1076,29 @@ palmieri_final <- subset(palmieri_manfiltered, !ids_to_exlude) validObject(palmieri_final) ``` -As we have just excluded probe IDs from the assay data, we now have -to also exclude them from the feature data `anno_palmieri`: +由于我们刚刚从化验数据中排除了 probe IDs， +我们现在还必须将它们从特征数据 `anno_palmieri`中排除： ```{r recallAnnoPalmieri} head(anno_palmieri) ``` -Recall that `fData` enables us to access the feature data of an expression set. -Until now, no feature data whatsoever is stored in the `fData(palmieri_final)`. -Only the row names are the row names of the assay data by default, which are the -PROBEIDs of the transcripts. +回想一下， `fData` 使我们能够访问表达式集的特征数据。 +到目前为止，没有任何特征数据存储在 `fData(palmieri_final)` 中。 +默认情况下，只有行名称是化验数据的行名称，它们是转录本的 +PROBEIDs。 -Therefore, we generate a column PROBEID in `fData(palmieri_final)` and -assign the row names of `fData(palmieri_final)` to it: +因此，我们在 `fData(palmieri_final)`中生成一列 PROBEID， +并为其分配 `fData(palmieri_final)` 的行名： ```{r excludeMultipleMappingsFromFeatureData} fData(palmieri_final)$PROBEID <- rownames(fData(palmieri_final)) ``` -Then, we left-join `fData(palmieri_final)`with `anno_palmieri`, which already contains -the columns "SYMBOL" and "GENENAME". -A left-join keeps the rows and columns of the first argument -and adds the corresponding column entries of the second argument: +然后，我们左手连接数据 `fData(palmieri_final)` 和 `anno_palmieri`， +它已经包含 columns "SYMBOL" 和 "GENENAME" 列。 +left-join 保留第一个参数的行和列， +并添加第二个参数的相应列条目： ```{r excludeMultipleMappingsFromFeatureData2} fData(palmieri_final) <- left_join(fData(palmieri_final), anno_palmieri) @@ -1110,75 +1110,75 @@ rownames(fData(palmieri_final)) <- fData(palmieri_final)$PROBEID validObject(palmieri_final) ``` -By left-joining with `anno_palmieri`, we thus add the "SYMBOL" and "GENENAME" -columns from `anno_palmieri` for only the PROBEIDs that are in -`fData(palmieri_final)` and thus get the feature Data for the filtered probes. +通过与 `anno_palmieri` 的左连接，我们因此只为 `anno_palmieri` 中的 +"SYMBOL" 和 "GENENAME" 列添加了 `fData(palmieri_final)` 中的 +PROBEID ，从而获得过滤探测的特征数据。 -## Building custom annotations +## 构建自定义注释 -Alternatively, one can re-map the probes of the array -to a current annotation. A workflow to do this for Illumina arrays is given in -@Arloth_2015. -Essentially, the individual probe sequences are re-aligned to an in-silico -"exome" that consists of all annotated transcript exons. +或者，可以将阵列的探针重新映射到当前注释。 +在 @Arloth_2015 中给出了为 Illumina 阵列执行此操作的工作流程。 +基本上， +将各个探针序列重新比对由包含所有注释的 +转录物外显子的计算机“外体组”。 -In any case, the package `r Biocpkg("pdInfoBuilder") ` can be used to build -customannotation packages for use with `r Biocpkg("oligo") `. -In order to do this, PGF / CLF files (called "Library files" on the Affymetrix -website) as well as the probeset annotations are required. -The probesets typically represent small stretches -of the genome (such as a single exon) and multiple probesets -are then used to form a transcript cluster. +在任何情况下，`r Biocpkg("pdInfoBuilder") ` 包都可用于构建用于 +`r Biocpkg("oligo") ` 的 customannotation 包。 +为此，需要 PGF / CLF 文件（在 Affymetrix 网站上称为“库文件”） +以及探针集注释。 +探针组通常代表基因组的小片段 +（例如单个外显子）， +然后使用多个探针组形成转录物簇。 -The CLF file contains information about the location of -individual probes on the array. The PGF file then contains the individual probe -sequences and shows the probeset they belong to. Finally, the probeset -annotation .csv then contains information about which probesets are used -in which transcript cluster. Commonly, multiple probesets are used in one -transcript cluster and some probesets are contained in multiple transcript -clusters. +CLF 文件包含有关阵列上各个探针位置的信息。 +然后，PGF 文件包含各个探测序列， +并显示它们所属的探针组。最后，probeset 注释 .csv +然后包含有关在哪个 transcript 集群中 +使用哪些 probeset 的信息。 +通常，在一个转录物簇中使用多个探针组， +并且一些探针组包含在多个转录物簇中。 -# Linear models +# 线性模型 -In order to analyse which genes are differentially expressed between -inflamed and non-inflamed tissue, we have to fit a linear model to our -expression data. Linear models are the "workhorse" for the analysis of -experimental data. They can be used to analyse almost arbitrarily -complex designs, however they also take a while to learn and understand -and a thorough description is beyond the scope of this workflow. +为了分析哪些基因在发炎组织和非发炎组织之间差异表达， +我们必须将线性模型拟合到我们的表达数据中。 +线性模型是分析实验数据的“主力”。 +它们可用于分析几乎任意复杂的设计， +但是它们也需要一段时间来学习和理解， +并且详尽的描述超出了本工作流程的范围。 -Mike Love's and Rafael Irizzary's ["Biomedical Data Science"](http://genomicsclass.github.io/book/) [@Irizarry_2015] is a very good resource, especially the section on [interactions and contrasts](http://genomicsclass.github.io/book/pages/interactions_and_contrasts.html). -It might also be helpful to learn some linear algebra to better understand -the concepts here. The Khan Academy offers helpful (and free) [online courses](https://www.khanacademy.org/math/linear-algebra). +Mike Love's 和 Rafael Irizzary's ["Biomedical Data Science"](http://genomicsclass.github.io/book/) [@Irizarry_2015] 是一个非常好的资源，尤其是关于[interactions and contrasts](http://genomicsclass.github.io/book/pages/interactions_and_contrasts.html).的部分。 +学习一些线性代数以更好地理解这里的概念也可能有所帮助。 +The Khan Academy 提供有用（和免费）的在线课程。 [online courses](https://www.khanacademy.org/math/linear-algebra). -## Linear models for microarrays +## 微阵列的线性模型 -We now apply linear models to microarrays. Specifically, we discuss -how to use the `r Biocpkg("limma") ` package for differential expression -analysis. The package is designed to analyze complex experiments involving -comparisons between many experimental groups simultaneously while remaining -reasonably easy to use for simple experiments. The main idea is to fit -a linear model to the expression data for each gene. Empirical Bayes and other -methods are used to borrow information across genes for the residual variance -estimation leading to "moderated" \(t\)-statistics, and stabilizing -the analysis for experiments with just a small number of arrays [@Smyth_2004]. -Conceptually, the final per gene variance is a mix of a prior variance and -the per gene variance. +我们现在将线性模型应用于微阵列。 具体来说，我们讨论如何使用 +`r Biocpkg("limma") ` 包进行差异表达式分析。 +该包用于分析复杂实验， +包括同时在许多实验组之间进行比较， +同时保持相当容易用于简单实验。 +主要思想是将线性模型拟合到每个基因的表达数据。 +经验贝叶斯和其他方法被用于借助基因信息来进行剩余方差估计， +从而导致“适度” \(t\)-statistics, 并稳定 +仅使用少量阵列的实验分析 [@Smyth_2004]. +从概念上讲， +每个基因的最终方差是先验方差和每个基因方差的混合。 -Typical experimental designs are disussed in chapter 9 of `r Biocpkg("limma") ` -"User Guide", which can be found on the Bioconductor landing page -of `r Biocpkg("limma") `. +典型的实验设计在 `r Biocpkg("limma") ` +"用户指南"的第9章中讨论，可以在`r Biocpkg("limma") ` +登陆页面上找到。 -In the following, we use appropriate design and contrast matrices for -our linear models and fit a linear model to each gene separately. +在下文中，我们为线性模型使用适当的设计和对比矩阵， +并分别为每个基因拟合线性模型。 -## A linear model for the data +## 数据的线性模型 -For the subsequent linear modelling of the data, we introduce the abbreviations -"UC" and "CD" for the disease types, and "non_infl." and "infl." -for the phenotypes, respectively: +对于随后的数据线性建模， +我们分别为疾病类型引入缩写“UC”和“CD”， +并为表型引入 “non_infl” 和 “infl”： ```{r introductionOfAbbreviations} individual <- @@ -1198,38 +1198,38 @@ disease <- "Crohn"), "CD", "UC") ``` -The original paper is interested in changes in transcription that occur -between inflamed and adjacent non-inflamed mucosal areas of the colon. -This is studied in both inflammatory bowel disease types. +原始论文对发炎和邻近的非发炎的结肠粘膜区域之间 +发生的转录变化感兴趣。 +这在两种炎症性肠病中都有研究。 -For building our linear model, we have to think about -which experimental variables we want to consider. -As we want to find differential expression -between non-inflamed and inflamed tissue, in principle, -those are the only two variables -we would have to consider. +为了构建我们的线性模型， +我们必须考虑我们想要考虑哪些实验变量。 +由于我们希望找到非发炎和发炎组织之间的 +差异表达， +原则上， +这些是我们必须考虑的唯一两个变量。 -However, since we have two arrays per individual patient, we have a -"Paired Samples" design (see section 9.4 of the `r Biocpkg("limma")` -user guide). This means that the samples might be biased by the person they come -from. Whenever a feature in an experimental setup -is expected to have a systematic influence on the result, blocking factors -on these features should be introduced. +但是，由于每个患者有两个阵列， +我们有一个“配对样本”设计 (参见 第 9.4 节 of the `r Biocpkg("limma")` +用户指南). 这意味着样本可能会受到来自的人的偏见。 +每当预期实验装置中的特征对结果 +具有系统性影响时， +应引入这些特征的阻塞因子。 -Thus, the first factor we need is a blocking factor for the individuals -that will absorb differences in expression between -them. Therefore, we block on patients, which means that the -patient IDs become variables of the linear model. +因此，我们需要的第一个因素是个体的阻滞因子， +它将吸收它们之间表达的差异。 +因此，我们阻止患者， +这意味着患者 ID 成为线性模型的变量。 -Then we create factors that give us the grouping for -the tissue types (non-inflamed and inflamed). +然后我们创建了一些因子， +这些因子可以为我们提供组织类型的分组（非发炎和发炎）。 -Finally, we create two design matrices, -one for each of the two diseases as we will analyze them separately -in order to follow the analysis strategy of the original paper -closely (one could also fit a joint model to -the complete data set; however, the two diseases might show very different -characteristics so that a joint fit might not be appropriate). +最后，我们创建两个设计矩阵， +分别针对两种疾病中的每一种， +因为我们将分别对它们进行分析， +以便密切关注原始论文的分析策略 +也可以将一个联合模型拟合到完整的数据集; +但是，两种疾病可能表现出非常不同的特征，因此联合适合可能不合适）。 ```{r createDesign, eval=TRUE, dependson="excludeMultipleMappings" } i_CD <- individual[disease == "CD"] @@ -1243,7 +1243,7 @@ colnames(design_palmieri_UC)[1:2] <- c("I", "nI") rownames(design_palmieri_UC) <- i_UC ``` -We can inspect the design matrices: +我们可以检查设计矩阵： ```{r inspectDesignMatrix, eval = TRUE, dependson="createDesign"} head(design_palmieri_CD[, 1:6]) @@ -1251,55 +1251,55 @@ head(design_palmieri_CD[, 1:6]) head(design_palmieri_UC[, 1:6]) ``` -In the design matrix, the rows represent the patient array, -and the columns are the variables we include in our linear model. The variables -correspond to our blocking factors: there are two for non-inflamed and inflamed -tissue, respectively, and one for each patient. "`i_UC2424`" for example is the -blocking variable of patient 2424; UC stands for the disease the patient is -suffering from. -For example, the first two rows of the design matrix `design_palmieri_CD` -correspond to the two arrays for individual "164". +在设计矩阵中，行表示患者阵列， +列是我们在线性模型中包含的变量。 +变量对应于我们的阻塞因素: 对于非发炎和发炎的组织分别有两个， +每个患者一个。例如， "`i_UC2424`" 是患者 2424 的阻断变量; + "`i_UC2424`" 是患者 2424 的阻断变量; +UC 代表患者患的疾病。 +例如，设计矩阵 `design_palmieri_CD` 的前两行对应于单个 +的两个数组。 -The design matrix entries are 0 or 1 and thus tell us -which variables are "active" for which sample: +设计矩阵条目为 0 或 1 +因此告诉我们哪个变量对于哪个样本是“活动的”： -It can be seen as a switch that turns the single variables "on" (with a 1 at the -corresponding position) or "off" (with a 0 at the corresponding position) for -each row, i. e. each patient array. If we take a closer look at the single rows, -we see that for each sample, there are two "ones" assigned: one to one of the -variables "nI" or "I" corresponding to the tissue the array came from, -and one to the corresponding patient-specific blocking variable. +它可以看作是一个开关， +用于将每个行中的单个变量“打开”（在相应位置处为 1）或“关闭”（在相应位置处为 0 ）， +即每个患者阵列。如果我们仔细查看单行， +我们会看到每个样本都分配了两个 “1”： +一个到一个变量 “nI” 或 “I” 对应于阵列来自的组织， +一个对应于相应的患者特异性阻断变量。 -Note that in the linear model, individual 164 serves as the baseline for all -other individuals and thus isn't included in the sample variables. +注意，在线性模型中，个体 164 用作所有其他个体的基线， +因此不包括在样本变量中。 -Note that instead of blocking on individuals, it would also be possible to use a -"mixed model" approach with the `duplicateCorrelation()`function from the -`r Biocpkg("limma")` package. -It has advantages over the "fixed patient effect model" presented here -in terms of applicability to more complicated experimental designs, -where we want to perform comparisons both within and between the patients (e.g. -comparing between the two diseases; "split-plot-designs"). +请注意，不是阻止个人，也可以使用 `r Biocpkg("limma")` 包中的 +`duplicateCorrelation()` 函数的 +“混合模型”方法。 +它比这里提出的“固定患者效应模型”具有优势， +适用于更复杂的实验设计， +我们希望在患者内部和患者之间进行比较 +（例如，比较两种疾病;“分裂情节设计”））。 -More information on it can be found in the [limma User's Guide](https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf) (section 17.3.6). -However, the above explained is more intuitive and is therefore used here. +有关它的更多信息可以在 [limma User's Guide](https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf) (section 17.3.6). + 然而，以上说明更直观，因此在此使用。 -Before heading on to find all differentially expressed genes for both diseases, -we will first have a look at how this works in principle for one gene. -We will fit the linear model for one gene and run a t-test in order -to see whether the gene is differentially expressed or not. +在继续找到这两种疾病的所有差异表达基因之前， +我们首先要看看这对于一种基因原理是如何起作用的。 +我们将拟合一个基因的线性模型并进行t-test， +以查看该基因是否差异表达。 -## Analysis of differential expression based on a single gene +## 基于单个基因的差异表达分析 -For linear model fitting and subsequent testing for differential expression -by t-test, we will pick the gene with the PROBEID 8164535. It has the gene -symbol *CRAT* and will be named as such in the following code. +对于线性模型拟合和随后通过 t-test 进行差异表达的测试， +我们将使用 PROBEID 8164535 挑选基因。 +它具有基因符号 *CRAT* 并将在以下代码中命名。 -**Illustration of the fitted linear model on the *CRAT* gene** +***CRAT* 基因拟合线性模型的插图** -Before fitting the linear model, we have a look at the expression intensities -of this gene for each patient in non-inflamed and inflamed tissue, -respectively: +在拟合线性模型之前， +我们分别研究了该基因在非发炎和发炎组织中 +的表达强度： ```{r visualizeExpressionChanges, fig.cap="Visualization of expression changes"} From b4730227286aba73199d1742a2a78a0e2fc44fbe Mon Sep 17 00:00:00 2001 From: Xiao W Date: Tue, 2 Jul 2019 13:14:05 -0400 Subject: [PATCH 05/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 238 ++++++++++++------------- 1 file changed, 119 insertions(+), 119 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index df06fd9..b20b3e3 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1317,22 +1317,22 @@ ggplot(data = crat_data, aes(x = tissue_CD, y = org_value, ggtitle("Expression changes for the CRAT gene") ``` -We see that overall, this gene is expressed less in inflamed tissue -(Figure \@ref(fig:visualizeExpressionChanges)). We also see that the absolute -expression intensities vary greatly between patients. -However, we have already taken care of this -problem by introducing blocking factors based on the individuals, which allows -us to compare the tissues within each individual as represented by -the single lines. - -If we had not blocked for individuals, the linear model would -treat them interchangably and a graphical depiction would only -include a single line. Since the individuals haver very different -baseline expression levels, this would lead to a very high variance of -the estimated fold changes. - -We now compute the variable coefficients by fitting a linear model. -We get a vector `crat_coef` with one entry for each variable. +总体而言，我们看到这种基因在炎症组织中的表达较少 +(Figure \@ref(fig:visualizeExpressionChanges))。我们还发现 +患者之间的绝对表达强度差别很大。 +然而，我们已经 +通过引入基于个体的阻断因子来解决这个问题， +这允许我们比较每个个体内的组织， +如单线所示。 + +如果我们没有阻止个体， +线性模型将可互换地对待它们， +并且图形描述将仅包括单行。 +由于个体具有非常不同的基线表达水平， +这将导致估计的倍数变化的非常高的方差。 + +我们现在通过拟合线性模型来计算变量系数。 +我们得到一个向量 `crat_coef` ，每个变量有一个条目。 ```{r obtainFitCRAT} crat_coef <- lmFit(palmieri_final[,disease == "CD"], @@ -1342,8 +1342,8 @@ crat_coef ``` -In order to now obtain the expression values fitted by the model, -we multiply the design matrix with this vector of coefficients `crat_coef`: +为了现在获得模型拟合的表达式值， +我们将设计矩阵乘以系数 `crat_coef`的向量： ```{r 2obtainFitCRAT} crat_fitted <- design_palmieri_CD %*% crat_coef @@ -1353,21 +1353,21 @@ colnames(crat_fitted) <- "fitted_value" crat_fitted ``` -Recall that for every row in the design matrix (i.e. every patient sample) -only the variables with a 1 in the design matrix are taken into account for -calculating the fitted expression value. +回想一下，对于设计矩阵中的每一行（即每个患者样本）， +仅考虑设计矩阵中具有1的变量 +来计算拟合表达值。 -This means that as output of the multiplication, we get -a vector `crat_fitted` whose entries are the sum of relevant -variable coefficients for each sample, respectively. +这意味着作为乘法的输出， +我们得到一个矢量 `crat_fitted`， +其条目分别是每个样本的相关变量系数的总和。 -For example, the fitted value for patient sample `2114_I.CEL` -is `r crat_coef["nI"] + crat_coef["i_CD2114"]`: it is the sum of the respective activated variable coefficients "nI" (`r crat_coef["nI"]`) -and "i_CD2114" (`r crat_coef["i_CD2114"]`). +例如，患者样本 `2114_I.CEL` 的拟合值是 +`r crat_coef["nI"] + crat_coef["i_CD2114"]`: 它是各个激活变量系数 "nI" (`r crat_coef["nI"]`) +和 "i_CD2114" (`r crat_coef["i_CD2114"]`) 的总和。 -Let's visualize the difference between non-inflamed and inflamed tissue -again after fitting: +让我们在拟合后 +再次可视化非发炎和发炎组织之间的差异： ```{r 3obtainFitCRAT, fig.cap="Expression changes for the CRAT gene"} crat_data$fitted_value <- crat_fitted @@ -1379,51 +1379,51 @@ ggplot(data = crat_data, aes(x = tissue_CD, y = fitted_value, ``` -Note that the difference of the fitted expression values between inflamed and -non-inflamed samples of one patient is the same for all patients and -is determined by the difference between the variable coefficients of -`I` (`r crat_coef["I"]`) and `nI` (`r crat_coef["nI"]`), -which is `r crat_coef["I"] - crat_coef["nI"]` +注意，一名患者的发炎和非发炎样本之间拟合表达值的差异对于所有患者是相同的， +并且由 +`I` (`r crat_coef["I"]`) 和 `nI` (`r crat_coef["nI"]`) +的变量系数之间的差异确定 +即 `r crat_coef["I"] - crat_coef["nI"]` (Figure \@ref(fig:3obtainFitCRAT)). -This is the case because the same blocking variable is activated by the design -matrix for both samples from a single patient, leading to a comparison within -patients only. These blocking variables correct the fitted tissue specific -expression values towards the expression levels of the individual patients. -Therefore the final estimate is like an average of all the within-individual -differences. +情况就是这样， +因为来自单个患者的两个样本的设计矩阵激活相同的阻塞变量， +导致仅在患者内进行比较。 +这些阻断变量使拟合的组织特异性表达值朝向个体患者的表达水平进行校正。 +因此， +最终估计值就像所有个体内部差异的平均值。 -The "difference" between non-inflamed and inflamed tissue of -`r crat_coef["I"] - crat_coef["nI"]` is actually a log2 fold change, as our -expression data is on the log2 scale. `r crat_coef["I"] - crat_coef["nI"]` -therefore is our log2 fold change for the CRAT gene. +非发炎组织和发炎组织之间的 “差异” +`r crat_coef["I"] - crat_coef["nI"]` 实际上是 log2 倍变化， +因为我们的表达数据是 log2 标度。因此，`r crat_coef["I"] - crat_coef["nI"]` +是 CRAT 基因的 log2 倍变化。 -**Differential expression analysis of the CRAT gene** +**CRAT 基因的差异表达分析** -In order to test whether the gene is differentially expressed or not, -a \(t\)-test with the null hypothesis that there is no difference -in the expression between non-inflamed and inflamed tissue is carried out. -Our blocking design is conceptually similar to a paired t-test for -which the statistic is given by: +为了测试基因是否差异表达， +进行具有零假设的 \(t\)-test， +即非发炎组织和发炎组织之间的表达没有差异。 +我们的阻塞设计在概念上类似于配对t检验， +其统计数据由下式给出： \[ t = \frac{\bar{d} }{s/\sqrt{n}} \] -Where, \(\bar{d}\) is the mean difference in expression values -between the individuals. -The paired t-test computes the variance \(s^2\) from the paired differences. -This is lower than the variance of a standard t-test and thus the paired -t-test has higher power as long as the expression values -for the same individual are correlated ([see e.g. the article on Wikipedia](https://en.wikipedia.org/wiki/Paired_difference_test)). +其中, \(\bar{d}\) 是个体之间表达值 +的平均差异。 +paired t-test 从配对差异计算方差 \(s^2\)。 +这低于标准 t-test 的方差， +因此只要同一个体的表达值相关， +配对t检验就具有更高的功效（参见例如维基百科上的文章）。 ([see e.g. the article on Wikipedia](https://en.wikipedia.org/wiki/Paired_difference_test)). -We thus have improved the power of the ordinary \(t\)-test -by reducing the variance via blocking on individuals. +因此，我们通过阻止个体来减少方差， +从而提高了普通 \(t\)-test 的能力。 -We now conduct the \(t\)-test on the linear model in order to find -out whether the difference between non-inflamed and inflamed -tissue differs significantly from 0: +我们现在对线性模型进行 \(t\)-test， +以确定非发炎组织和发炎组织之间的差异 +是否与 0 显着不同: ```{r tTest} crat_noninflamed <- na.exclude(crat_data$org_value[tissue == "nI"]) @@ -1432,28 +1432,28 @@ res_t <- t.test(crat_noninflamed ,crat_inflamed , paired = TRUE) res_t ``` -We get a low p-value close to `r round(res_t[["p.value"]], 4) ` -and thus can conclude that the *CRAT* gene is differentially -expressed between non-inflamed and inflamed tissue. +我们得到接近 `r round(res_t[["p.value"]], 4) ` 的低 p-value 值， +因此可以得出结论 *CRAT* 基因 +在非发炎和发炎组织之间差异表达。 -Note that the p-value isn't exactly the same one as below when analyzing the -differential expression of all genes. This is due to the variance moderation -performed by `r Biocpkg("limma")` . +注意，在分析所有基因的差异表达时 p-value p值与下面的p值不完全相同。 +这是由于由 `r Biocpkg("limma")` +执行的方差调节。 -## Contrasts and hypotheses tests +## 对比和假设测试 -We now fit the linear model for all genes and -define appropriate contrasts to test hypotheses of interest. +我们现在适合所有基因的线性模型， +并定义适当的对比来检验感兴趣的假设。 -We want to compare the inflamed to the non-inflamed tissue. -Thus, we create a contrast matrix consisting of only one contrast "I-nI": -`r Biocpkg("limma")`'s function -`makeContrasts` creates this matrix from -a symbolic description of the contrast of -interest. +我们想比较发炎的和非发炎的组织。 +因此，我们创建一个仅由一个对比度 "I-nI" 组成的对比矩阵: +`r Biocpkg("limma")` 的函数 +`makeContrasts` 根据感兴趣的 +对比度的符号 +描述创建该矩阵。 -We now fit a linear model to our data and apply the `contrasts.fit()` function -to it in order to find genes with significant differential expression between non-inflamed and inflamed tissue: +我们现在将线性模型拟合到我们的数据中并对其应用对比度 `contrasts.fit()` 函数， +以便找到在非发炎组织和发炎组织之间具有显着差异表达的基因： @@ -1471,44 +1471,44 @@ palmieri_fit_UC <- eBayes(contrasts.fit(lmFit(palmieri_final[,disease == "UC"], contrast_matrix_UC)) ``` -We applied the empirical Bayes variance moderation method to the model -via the `eBayes()` function, which -computes moderated \(t\)-statistics. In microarray analysis, -the number of arrays often is quite small, -and thus variance estimation is difficult. -Using a combination of the per-gene-variance -and a prior variance we can improve the -variance estimate, hence the term "moderation". "Empirical Bayes" -means that the prior is estimated from the data. - -The result of the `eBayes()` step is that the individual variances are -shrunken towards the prior value. - -## Extracting results - -Finally, we extract the number of differentially expressed genes. -Results can be extracted by use of the `topTable` function. -We extract the results for both Crohn's disease and ulcerative colitis, and the -results are sorted by their absolute \(t\)-statistics. -As a diagnostic check, we also plot the p-value histogram: We expect a uniform -distribution for the p-values that correspond to true null hypotheses, while -a peak near zero shows an enrichment for low p-values corresponding to differentially -expressed (DE) genes. +我们通过 `eBayes()` 函数将 +经验贝叶斯方差调节方法应用于模型， +该函数计算调节的t统计量。 +在微阵列分析中，阵列的数量通常非常小， +因此难以估计方差。 +使用每基因方差和先验方差的组合， +我们可以改进方差估计， +因此术语“适度”， +“经验贝叶斯”意味着先验是根据数据估计的。 + +`eBayes()` 步骤的结果是 +个体差异朝着先前值收缩。 + +## 提取结果 + +最后，我们提取差异表达基因的数量。 +可以使用 `topTable` 函数提取结果。 +我们提取克罗恩病和溃疡性结肠炎的结果， +结果按其绝对t统计量排序。 +作为诊断检查， 我们还绘制了 p 值直方图： +我们期望对应于真零假设的 p 值的均匀分布， +而接近零的峰显示对应于差异表达的低 p 值 +的富集（DE）基因。 -Note that if the p-value distribution for a dataset is very different -from the ones in the histograms below, this might lead to quality loss in the -subsequent analysis. Reasons for a divergent p-value-distribution might be -batch effects or a lack of consideration of other blocking factors in the design -matrix. Thus, if the p-value is not as expected, -try to include possible blocking factors and batches and rerun the analysis. -If this does not help, empirical Bayes / null estimation methods -for multiple testing are useful. - -A good starting point to learn about these methods is the article on -false discovery rate estimation by Korbininan Strimmer [@Strimmer_2008] -and chapter 1--6 of Efron's book on Large-Scale Inference [@Efron_2010], -as well as the blog-post on "How to interpret a p-value histogram" by David -Robinson [@Robinson_2014]. +请注意，如果数据集的 p 值分布与下面的直方图中的 p 值分布非常不同， +则可能会导致后续分析中的质量损失。 +p 值分布不同的原因可能是 +批量效应或设计矩阵中缺乏对其他阻塞因素的考虑。 +因此，如果 p 值不符合预期，请尝试包括可能的阻塞因子和批次， +然后重新运行分析。 +如果这没有帮助， +用于多次测试的经验贝叶斯/零估计方法是有用的。 + +了解这些方法的一个很好的起点是关于 +Korbininan Strimmer [@Strimmer_2008] 的错误发现率估计和 +Efron's 的大规模推理书 [@Efron_2010], +以及关于“如何”解释 p值直方图的博客文章 by David +Robinson [@Robinson_2014]。 ```{r extractResultsCD, eval = TRUE, dependson="createContrastMatrixAndFitModel", message=FALSE, fig.cap="Histogram of the p–values for Crohn’s disease."} table_CD <- topTable(palmieri_fit_CD, number = Inf) @@ -1529,11 +1529,11 @@ hist(table_UC$P.Value, col = brewer.pal(3, name = "Set2")[2], ``` -## Multiple testing FDR, and comparison with results from the original paper +## 多次测试 FDR，并与原始论文的结果进行比较 -In the original paper, a p-value of 0.001 was used as a significance -cutoff. Using this we get `r nrow(subset(table_UC, P.Value < 0.001))` -genes identified as differentially expressed for UC: +在原始论文中，p 值为 0.001 用作显着性截止值。 +使用这个，我们得到 `r nrow(subset(table_UC, P.Value < 0.001))` +个被鉴定为 UC 差异表达的基因： ```{r pCutForCD} From 52c55a19541a29001fdbfdb9e3bddb9574bb945f Mon Sep 17 00:00:00 2001 From: Xiao W Date: Tue, 2 Jul 2019 13:16:10 -0400 Subject: [PATCH 06/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index b20b3e3..179f395 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1295,7 +1295,7 @@ UC 代表患者患的疾病。 我们将使用 PROBEID 8164535 挑选基因。 它具有基因符号 *CRAT* 并将在以下代码中命名。 -***CRAT* 基因拟合线性模型的插图** +** *CRAT* 基因拟合线性模型的插图** 在拟合线性模型之前， 我们分别研究了该基因在非发炎和发炎组织中 From 8131b881934f2c261c7bdf2067f945e611f2435f Mon Sep 17 00:00:00 2001 From: Xiao W Date: Tue, 2 Jul 2019 13:18:30 -0400 Subject: [PATCH 07/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index 179f395..348e999 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1295,7 +1295,7 @@ UC 代表患者患的疾病。 我们将使用 PROBEID 8164535 挑选基因。 它具有基因符号 *CRAT* 并将在以下代码中命名。 -** *CRAT* 基因拟合线性模型的插图** +**插图：*CRAT* 基因拟合线性模型** 在拟合线性模型之前， 我们分别研究了该基因在非发炎和发炎组织中 From f494347ab76605d524ae2c1a0851fb8dadf0f695 Mon Sep 17 00:00:00 2001 From: Xiao W Date: Wed, 3 Jul 2019 11:05:27 -0400 Subject: [PATCH 08/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 198 ++++++++++++------------- 1 file changed, 99 insertions(+), 99 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index 348e999..dd0a333 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1541,62 +1541,62 @@ nrow(subset(table_UC, P.Value < 0.001)) ``` -However, it is impossible -to determine a precise bound on the number of false positive genes -in this list. All that we can say using p-values is that we have -at most `r nrow(table_UC)` (total number of tests) * 0.001 -= `r nrow(table_UC)* 0.001` false positive genes in our list. Therefore, -by choosing a p-value cutoff of 0.001, as much as +然而， +不可能确定该列表中假阳性基因数量的精确界限。 +我们所说的使用 p 值的是我们列表中最多 +`r nrow(table_UC)` (测试总数) * 0.001 += `r nrow(table_UC)* 0.001` 个假阳性基因。 +因此，通过选择 0.001 的 p 值截止值， +鉴定为差异表达的基因中多达 `r round(nrow(table_UC)* 0.001 / nrow(subset(table_UC, P.Value < 0.001))*100, 2)`% -of our genes -identified as differentially expressed might be false positives. +可能是假阳性。 -Thus, we can see that the "raw" p-values are very "liberal" when -looking at many tests simultaneously. We therefore need error -rates adapted to the multiple testing situation. By far the -most popular one in molecular biology is the **false discovery -rate** or FDR for short. It is the percentage of false positives -among all positives. As we have seen, the FDR of our genes list -using a simple p-value cutoff might be quite high. +因此，我们可以看到，当同时观察许多测试时，“原始” p 值非常“自由”。 +因此，我们需要适应多种测试情况的错误率。 +到目前为止，分子生物学中最受欢迎的是错误发现率 **false discovery rate** 或简称 FDR。 +这是所有积极因素中误报的百分比。 +正如我们所看到的， +使用简单的 p 值截止值列出的基因列表 +的 FDR 可能非常高。 -On the other hand, we can see a clear peak in the p-value histogram -(Figure \@ref(fig:extractResultsUC)), caused by the differentially -expressed genes. There we expect the actual FDR of our list -to be lower. +另一方面，我们可以看到由差异表达的基因引起的 p 值直方图中的清晰峰值 +(Figure \@ref(fig:extractResultsUC)), +在那里， +我们预计列表的实际 FDR 会更低。 -The FDR at a given cutoff is given by the "adjusted" p-value in -the results table. +给定截止值下的 FDR +由结果表中的“调整后的” p 值给出。 ```{r FDRforUC} tail(subset(table_UC, P.Value < 0.001)) ``` -The adjusted p-value for a raw p-value of 0.001 in -the table is `r max(subset(table_UC, P.Value < 0.001)$adj.P.Val)`, -which is an order of magnitude lower than the FDR we can infer -from p-values alone. +表中原始 p 值 0.001 +的调整后的 p 值为 `r max(subset(table_UC, P.Value < 0.001)$adj.P.Val)`, +这比我们仅从 p 值推断 +的 FDR 低一个数量级。 -So although this is not recommended in general, -we also use a p-value cutoff at 0.001 in the following -in order to be able to compare our workflow results to the paper results. +因此，尽管一般不推荐这样做， +但我们还在下面使用 0.001 的 p 值截止值， +以便能够将我们的工作流程结果与纸张结果进行比较。 -The paper results can be downloaded as excel files from -[http://links.lww.com/IBD/A795](http://links.lww.com/IBD/A795) and -should be saved as an .xlsx file named `palmieri_DE_res.xlsx` in -your working directory. +论文结果可以从 +[http://links.lww.com/IBD/A795](http://links.lww.com/IBD/A795) +下载为 .xlsx 文件， +并应保存为工作目录中名为 `palmieri_DE_res.xlsx` -Note that while the paper uses p-value cutoffs, it also reports the -corresponding FDRs (just as we did for the UC data here). +请注意，虽然本文使用 p 值截止值， +但它也会报告相应的 FDR（就像我们在此处对 UC 数据所做的那样）。 -For a p-value cutoff of 0.001, the corresponding FDRs are 0.05 in -Crohn's disease and 0.02 in ulcerative colitis. There are four tables in total, -giving the list of up and downregulated genes in CD and UC, respectively. -We calculate the overlap between our results and the ones from the -paper as the ratio of the genes that were found in both analyses -and the genes that were only found in the paper. +对于 p 值截止值为 0.001，相应的 FDR 在克罗恩病中为 0.05， +在溃疡性结肠炎中为 0.02。 +总共有四个表， +分别给出 CD 和 UC 中上调和下调基因的列表。 +我们计算了我们的结果与论文中的结果之间的重叠， +即两种分析中发现的基因与仅在论文中发现的基因的比率。 -We also calculate the total number of diffentially expressed -genes that we find in our workflow analysis. +我们还计算了在工作流分析中 +发现的差异表达基因的总数。 ```{r compareDEgenes, dependson=c("extractResultsUC", "extractResultsCD")} fpath <- system.file("extdata", "palmieri_DE_res.xlsx", package = "maEndToEnd") @@ -1624,27 +1624,27 @@ total_genenumber_CD total_genenumber_UC ``` -We find `r total_genenumber_CD` (CD) and `r total_genenumber_UC` (UC) -differentially expressed genes ("DE-genes"). +我们发现 `r total_genenumber_CD` (CD) 和 `r total_genenumber_UC` (UC) +差异表达的基因（“DE-基因”）。 -In the paper, 298 (CD) and 520 (UC) DE-genes were found for the two diseases -at the same cutoff. This higher number of DE-genes identified is probably -due to the increased power of the blocking according to the individuals -and the moderated variance estimation that `r Biocpkg("limma") ` performs. +在该论文中，在相同的截止值下发现了两种疾病的 298 (CD) and 520 (UC) DE基因。 +鉴定的这种更高数量的 DE-基因可能是由于 +根据个体的阻断功率增加以及 +`r Biocpkg("limma") ` 执行的调节方差估计。 -We see that we get a moderate overlap of `r overlap_CD` for CD and -`r overlap_UC` for UC, showing that both analyses lead to -somewhat comparable results. +我们看到我们得到 CD 的 `r overlap_CD` 和 +UC 的 `r overlap_UC` 的适度重叠， +表明两种分析都会导致一些可比较的结果。 -## Visualization of DE analysis results - volcano plot +## DE 分析结果的可视化 - 火山图 -For a visualization of the differentially expressed genes, we create a volcano -plot, which is commonly used to summarize the results of a differential expression analysis in a single figure. +为了显示差异表达的基因，我们创建了一个火山图， +通常用于总结单个图中差异表达分析的结果。 -For a better overview, we only show gene symbols of genes with -a fold change greater than 1, which we define in the `volcano_names` object. The -`highlight` option in the `volcanoplot` function is set to 100 and thus only -labels the 100 genes with the lowest p-values. +为了更好地概述，我们只显示倍数大于 1 的基因的基因符号， +我们在 `volcano_names` 对象中定义。 +`volcanoplot` 函数中 `highlight` 高亮选项设置为 100， +因此仅标记具有最低 p 值的 100 个基因。 ```{r VolcanoPlot, fig.cap="Volcano plot of the DE-genes", fig.height=8, fig.width=7} @@ -1657,51 +1657,51 @@ volcanoplot(palmieri_fit_CD, coef = 1L, style = "p-value", highlight = 100, xlab = "Log2 Fold Change", ylab = NULL, pch=16, cex=0.35) ``` -We can now do a little research on the biological function of genes that show a -high foldchange, for example the gene with the symbol S100A8 on the right side -of the plot (Figure \@ref(fig:VolcanoPlot)). If we search for this gene symbol on -[genecards.org](https://www.genecards.org), we find that it encodes for a -protein that builds a pro-inflammatory complex in association with another -protein. +我们现在可以对显示高倍数变化的基因的生物学功能进行一些研究， +例如在图的右侧具有符号 S100A8 的基因 +(Figure \@ref(fig:VolcanoPlot)). 如果我们在 +[genecards.org](https://www.genecards.org) 上搜索这个基因符号， +我们发现它编码的蛋白质构建了 +与另一种蛋白质相关的促炎复合物。 -# Gene ontology (GO) based enrichment analysis +# 基因本体论（GO）为基础的浓缩分析 -As discussed above, it is recommended to use an FDR cutoff in differential -expression analysis rather than a p-value cutoff, -since this way you control an explicitly defined error rate -and the results are easier to interpret and -to compare. For the following enrichment analysis, we create tables with -differentially expressed genes for CD and UC, respectively, -and choose an FDR cutoff of 10\%. Here, we focus on the CD subset of the data. +如上所述，建议在差分表达式分析中使用 FDR 截止而不是 p 值截止， +因为这样您可以控制明确定义的错误率， +并且结果更易于解释和比较。 +对于以下富集分析， +我们分别创建具有 CD 和 UC 的差异表达基因的表， +并选择 10％ 的 FDR 截止值。 +在这里，我们关注数据的 CD 子集。 ```{r FDRcontrolledDEgenes, dependson=c("extractResultsUC", "extractResultsCD"), eval=TRUE} DE_genes_CD <- subset(table_CD, adj.P.Val < 0.1)$PROBEID ``` -We can now try to characterize the identified differentially expressed genes -more in detail by performing a GO enrichment analysis. Essentially the -gene ontology ([http://www.geneontology.org/](http://www.geneontology.org/)) is -a hierarchically organized collection of functional gene sets +我们现在可以通过进行 GO 富集分析来更详细地表征鉴定的差异表达基因。 +本质上， +基因本体论 ([http://www.geneontology.org/](http://www.geneontology.org/)) +是一个分层组织的功能基因集合 [@Ashburner_2000; @GO_2015; @du_Plessis_2011]. -## Matching the background set of genes +## 匹配背景基因组 -The function ` genefinder ` from the `r Biocpkg("genefilter") ` -package [@Bourgon_2010] -will be used to find a background set of genes that are similar in expression -to the differentially expressed genes. We then check whether -the background has roughly the same distribution -of average expression strength as the foreground. +来自 `r Biocpkg("genefilter") `包 [@Bourgon_2010] 的 +` genefinder ` 函数 +将用于找到与差异表达基因的表达相似的基因背景组。 +然后我们检查背景是否具有 +与前景大致相同的 +平均表达强度分布。 -We do this in order not to select a biased background since the gene set testing -is performed by a simple Fisher test on a 2x2 table. Note that this approach -is very similar to commonly used web tools like GOrilla [@Eden_2009]. +我们这样做是为了不选择有偏见的背景， +因为基因集测试是在 2x2 表上通过简单的 Fisher 测试进行的。 +请注意，此方法与常用的 Web 工具 GOrilla [@Eden_2009] 非常相似。 -For every differentially expressed gene, we try to find genes with similar -expression with `genefinder`. The `genefinder` function returns a list with two -elements for each gene: one with the indices of the -background genes found and one with the distances to the DE-genes: +对于每个差异表达的基因，我们尝试用 +`genefinder` 找到具有相似表达的基因。`genefinder` 函数返回一个列表， +其中包含每个基因的两个元素： +一个具有发现的背景基因的索引，另一个具有与 DE基因的距离： ```{r GOAnalysisCreateBackgrounds, eval=TRUE, warning=FALSE, message=FALSE} back_genes_idx <- genefilter::genefinder(palmieri_final, @@ -1709,23 +1709,23 @@ back_genes_idx <- genefilter::genefinder(palmieri_final, method = "manhattan", scale = "none") ``` -We have to extract the PROBEIDs, which correspond to the indices. -We do that by using the `sapply` -function, which gives us a single matrix with the DE-genes as column names -and the PROBEIDs of the corresponding background genes in the cells below: +我们必须提取与索引相对应的 PROBEIDs, 我们通过使用 +`sapply` 函数来实现这一点， +它为我们提供了一个单独的矩阵， +其中 DE-genes 作为列名，以及下面单元格中相应背景基因的 PROBEID： ```{r GOAnalysisCreateBackgrounds2, eval=TRUE, warning=FALSE, message=FALSE} back_genes_idx <- sapply(back_genes_idx, function(x)x$indices) ``` -We then create a vector `back_genes` containing all background gene PROBEIDs: +然后我们创建一个包含所有背景基因 PROBEID 的载体 `back_genes`： -In order to eliminate foreground genes, i.e. DE-genes, from the `back_genes` -set, we use the `setdiff` function. It returns all elements from the first -argument (`back_genes`) that are not part of the second argument -(`DE_genes_CD`). +为了从 `back_genes` 集中消除前景基因，即 DE 基因， +我们使用 `setdiff` 函数。 它返回第一个参数 +(`back_genes`) 中不属于第二个参数 +(`DE_genes_CD`) 的所有元素。 With the `intersect` function, we verify that we were successful: it should return 0, as there shouldn't be any intersect anymore between `back_genes`and `DE_genes_CD`: From c3e690548f76cec45fd82621561824008ff54930 Mon Sep 17 00:00:00 2001 From: Xiao W Date: Wed, 3 Jul 2019 11:44:48 -0400 Subject: [PATCH 09/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 190 ++++++++++++------------- 1 file changed, 95 insertions(+), 95 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index dd0a333..8b9383d 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1726,9 +1726,9 @@ back_genes_idx <- sapply(back_genes_idx, function(x)x$indices) 我们使用 `setdiff` 函数。 它返回第一个参数 (`back_genes`) 中不属于第二个参数 (`DE_genes_CD`) 的所有元素。 -With the `intersect` function, we verify that we were successful: -it should return 0, as there shouldn't be any intersect anymore between -`back_genes`and `DE_genes_CD`: +使用 `intersect` 函数，我们验证我们是否成功： +它应该返回 0，因为 +`back_genes`和 `DE_genes_CD` 之间不应该有任何交叉： ```{r GOAnalysisCreateBackgrounds3, eval=TRUE, warning=FALSE, message=FALSE} back_genes <- featureNames(palmieri_final)[back_genes_idx] @@ -1739,11 +1739,11 @@ intersect(back_genes, DE_genes_CD) length(back_genes) ``` -We create a multidensity plot with mean expression on the x-axis -and curves for all genes, foreground genes and background genes, respectively -(Figure \@ref(fig:multidensityPlot)). We want to see whether the background -genes show a plot similar to the foreground genes so that the background -is not biased for the gene enrichment analysis: +我们创建了一个多密度图， +其中 x 轴分别为平均表达，所有基因，前景基因和背景基因分别为曲线 +(Figure \@ref(fig:multidensityPlot))。 +我们想看看背景基因是否显示出与前景基因类似的图， +这样背景基因浓缩分析就没有偏差： ```{r multidensityPlot, fig.cap="Selecting a background set of genes for the gene ontology analysis."} multidensity(list( @@ -1755,35 +1755,35 @@ multidensity(list( main = "DE genes for CD-background-matching") ``` -When comparing the "background gene" curve to the "foreground gene" curve, -we see a similar curve shape, indicating a sensible background matching -(Figure \@ref(fig:multidensityPlot)). Note that the right-shift of the -"foreground-gene" curve in comparison to the "background-gene" curve -indicates that DE-genes are generally very -highly expressed, so that it wasn't possible to find background-genes -with exactly equal overall expression distribution. +当将“背景基因”曲线与“前景基因”曲线进行比较时， +我们看到类似的曲线形状，表明合理的背景匹配 +(Figure \@ref(fig:multidensityPlot))。 +请注意，与“背景基因”曲线相比，“前景 - 基因”曲线的右移 +表明DE基因通常表达非常高， +因此无法找到完全相同的背景基因。 +整体表达分布。 -The "all gene" curve has the leftmost curve maximum; -this can be explained by a high number of lowly expressed genes -in all samples and shows that a background matching is sensible -in order to avoid biases. +“全基因” 曲线最左边的曲线最大； +这可以通过所有样品中的 +大量低表达基因来解释， +并且表明背景匹配是合理的以避免偏差。 -For the actual testing of which GO gene sets are enriched in inflamed tissue, -we use the `r Biocpkg("topGO")` package which implements a nice interface -to Fisher testing and also has additional algorithms -taking the GO structure into account, by e.g. only reporting the most specific -gene set in the hierarchy [@Alexa_2006]. +对于 GO 基因组在炎症组织中富集的实际测试， +我们使用 `r Biocpkg("topGO")` 包，其实现了与 +Fisher 测试的良好界面， +并且还具有考虑 GO 结构的其他算法， +例如，仅报告层次结构中最具体的基因集 [@Alexa_2006]. -The GO has three top ontologies: Cellular component (CC), biological processes -(BP), and molecular function (MF). For illustrative purposes we limit ourselves -to the BP category here. +GO 有三个顶级本体：细胞成分（CC），生物过程（BP）和分子功能（MF）。 +出于说明目的， +我们将此限制在 BP 类别中。 -## Running topGO +## 运行 topGO -topGO requires a topGOdata object containing the necessary information for the -analysis. We follow the steps described in the topGO vignettes: First, -we will create a named vector `all_genes` with all genes to be analyzed, i.e. -DE-genes and background genes: +topGO 需要一个 topGOdata 对象，其中包含分析所需的信息。 +我们按照 topGO 小插图中描述的步骤进行操作： +首先，我们将创建一个命名向量 `all_genes` 其中包含所有要分析的基因， +即 DE 基因和背景基因： ```{r createFactorOfInterestingGenes, dependson="GOAnalysisCreateBackgrounds", eval=TRUE} gene_IDs <- rownames(table_CD) @@ -1795,56 +1795,56 @@ all_genes <- factor(as.integer(in_selection[in_universe])) names(all_genes) <- gene_IDs[in_universe] ``` -The following steps were carried through: +执行以下步骤： - 1. we created an `in_universe` vector by using the `%in%` matching function. - We want to know which elements from `gene_IDs` are also contained in - `DE_genes_CD` and `back_genes`, as the latter two are our gene universe - we use for enrichment analysis. - We got a vector `in_universe` with the length of `gene_IDs` that has the - entry `TRUE` when the corresponding gene in `gene_IDs` could be also found in - `DE_genes_CD` or `back_genes`, and `FALSE` otherwise. + 1. 我们使用`%in%` matching 函数创建了一个 `in_universe`向量。 + 我们想知道来自 `gene_IDs` 的哪些元素也包含在 + `DE_genes_CD` 和 `back_genes` 中， + 因为后两者是我们用于浓缩分析的基因世界。 + 我们得到一个带有 `gene_IDs` 长度的的载体 `in_universe`， + 当 `gene_IDs` 中的相应基因也可以在 `DE_genes_CD` 或 `back_genes` 中中找到时， + 其条目为 `TRUE`，否则为 `FALSE`。 - 2. We did the same for our DE-genes and call this vector `in_selection`. + 2. 我们为我们的 DE 基因做了同样的事情，并将这个载体称为 `in_selection`. - 3. We created the `all_genes` vector: - a) First, we selected all the elements from `in_selection` - that are `TRUE` in `in_universe` by applying - `all_genes <- in_selection[in_universe]`. - b) Then, we converted the elements in `all_genes` from `TRUE` and `FALSE` - to 0 and 1 by converting the vector to an integer vector. - This way, each element in the vector is a 0 if the corresponding gene is a - background gene and a 1 if the corresponding gene is a DE-gene. Also, we - converted the vector to a factor. - c) We named the vector elements with the corresponding `gene_IDs`. - - -We now initialize the `r Biocpkg("topGO") ` data set, using the GO annotations -contained in the annotation data base for the chip we are using. The `nodeSize` -parameter specifies a minimum size of a GO category we want to use: i.e. here, -categories with less than 10 genes are not included in the testing. + 3. 我们创建了 `all_genes` 向量： + a）首先，我们通过应用 `all_genes <- in_selection[in_universe]` + 选择 `in_selection` 中 + `in_universe` 中 `TRUE`的所有元素。 + b) 然后，我们通过将向量转换为整数向量，将 `all_genes` 中的元素从 `TRUE` 和 `FALSE` + 转换为 0 和 1。 + 这样，如果相应的基因是背景基因， + 载体中的每个元件是 0，如果相应的基因是 DE 基因则是 1。 + 此外，我们将矢量转换为因子。 + c) 我们用相应的 `gene_IDs` 命名了 vetor 元素。 + + +我们现在使用我们正在使用的芯片的注释数据库中包含的 GO 注释初始化 `r Biocpkg("topGO") ` 数据集。 +`nodeSize` 参数指定我们想要使用的 GO 类别的最小大小； +即，此处测试中不包括具有少于 +10 个基因的类别。 ```{r createTopGODataSet, dependson="createFactorOfInterestingGenes", eval=TRUE, message = FALSE, warning=FALSE} top_GO_data <- new("topGOdata", ontology = "BP", allGenes = all_genes, nodeSize = 10, annot = annFUN.db, affyLib = "hugene10sttranscriptcluster.db") ``` -Now the tests can be run. `r Biocpkg("topGO") ` offers a wide range of options, -for details see the paper [@Alexa_2006] or the package vignette. +现在可以运行测试了。`r Biocpkg("topGO") ` 包提供多种选择，有关详细信息， +请参阅 [@Alexa_2006] 或包插图。 -We run two common tests: an ordinary Fisher test for every GO category, and the -"elim" algorithm, which tries to incorporate the hierarchical structure of the -GO and tries to "decorrelate" it in order to report the most specific -significant term in the hierarchy. +我们运行两个常见的测试： 每个 GO 类别的普通 Fisher 测试，以及 +"elim" 算法，它试图结合 +的层次结构并 +尝试“去相关”它以报告最具体的重要术语。 -The algorithm starts processing the nodes / GO categories -from the highest (bottommost) level and then iteratively -moves to nodes from a lower level. If a node is scored as significant, -all of its genes are marked as removed in all ancestor nodes. -This way, the "elim" algorithm aims at finding the most specific node -for every gene. +该算法从最高（最低）级别开始处理节点/ GO 类别， +然后从较低级别迭代地移动到节点。 +如果节点评分为显着， +则其所有基因在所有祖先节点中被标记为已移除。 +这样， "elim" 算法旨在为 +每个基因找到最特定的节点。 -The test uses a 0.01 p-value cutoff by default. +该测试默认使用 0.01 p-value p值截止值。 ```{r runtopGOTests, results='hide', eval=TRUE, dependson = "createTopGODataSet", message = FALSE} result_top_GO_elim <- @@ -1853,16 +1853,16 @@ result_top_GO_classic <- runTest(top_GO_data, algorithm = "classic", statistic = "Fisher") ``` -We can now inspect the results. We look at the top 100 GO categories according -to the "Fisher elim" algorithm. The function `GenTable` produces -a table of significant GO categories, the function `printGenes` -gives genes annotated to them; the significant ones are denoted with a "2" -in the "raw p-value" column, the non-significant ones with a "1". -We therefore select `raw p-value == 2`. +我们现在可以检查结果。 我们根据 “Fisher elim” 算法查看前 100 个 GO 类别。 +`GenTable` 函数生成一个重要类别的表格， +`printGenes` 函数给出了注释它们的基因; +重要的那些在 “原始p值” 列中用 “2” 表示， +非重要的用 “1” 表示。 +因此，我们选择原始p值 `raw p-value == 2`. -Note that we do not get the actual p-values here because our `all_genes` vector -doesn't contain this information; it only tells us -whether a gene is differentially expressed or not. +请注意，我们在这里没有得到实际的p值， +因为我们的 `all_genes` 向量不包含这些信息; +它只告诉我们基因是否差异表达。 ```{r processtopGOResults, eval=TRUE, dependson="runtopGOTests"} res_top_GO <- GenTable(top_GO_data, Fisher.elim = result_top_GO_elim, @@ -1880,11 +1880,11 @@ res_top_GO$sig_genes <- sapply(genes_top_GO, function(x){ head(res_top_GO[,1:8], 20) ``` -## Visualization of the GO-analysis results +## GO 分析结果的可视化 -A graph of the results can also be produced. Here we visualize the three most -significant nodes according to the Fisher elim algorithm in the context of -the GO hierarchy. +还可以产生结果图。 +在这里， +我们根据 GO 层次结构中的 Fisher elim 算法可视化三个最重要的节点。 ```{r graphOfResults, fig.height = 6, eval=TRUE, results='hide', dpi=600, fig.cap="Significantly enriched GO nodes in the GO hierarchy"} @@ -1892,15 +1892,15 @@ showSigOfNodes(top_GO_data, score(result_top_GO_elim), firstSigNodes = 3, useInfo = 'def') ``` -We can see that indeed GO categories related to inflammation, signalling and -immune response come up as significant (Figure \@ref(fig:graphOfResults)) -Gene set enrichment analysis has been a field of very extensive -research in bioinformatics. For additional approaches -see the `r Biocpkg("topGO") ` vignette and the references therein -and also in the [GeneSetEnrichment view](http://bioconductor.org/packages/release/BiocViews.html#___GeneSetEnrichment). +我们可以看到确实与炎症， +信号传导和免疫反应相关的 GO 类别显着 (Figure \@ref(fig:graphOfResults)) +基因集富集分析一直是生物信息学中广泛研究的领域。 +有关其他方法， +请参阅 `r Biocpkg("topGO") ` 小插图及其中的参考文献以及 +[GeneSetEnrichment view](http://bioconductor.org/packages/release/BiocViews.html#___GeneSetEnrichment). -# A pathway enrichment analysis using reactome +# 使用 reactome 的途径富集分析 The package `r Biocpkg("ReactomePA") ` offers the possibility to test enrichment of specific pathways using the free, open-source, curated and peer reviewed @@ -1943,7 +1943,7 @@ head(as.data.frame(reactome_enrich))[1:6] Note that we trimmed pathway names to 20 characters. -## Visualizing the reactome based analysis results +## 可视化基于反应组的分析结果 The top pathways can be displayed as a bar chart that displays all categories with a p-value below the specified cutoff (Figure \@ref(fig:reactomeBar)). @@ -1973,7 +1973,7 @@ pathways. This allows the mapping of the expression data at hand to known regulatory interactions. -# Session information +# 会话信息 As the last part of this document, we call the function *sessionInfo*, which reports the version numbers of R and all the packages used in @@ -1997,7 +1997,7 @@ sessionInfo() ``` -# Acknowledgements +# 致谢 The authors would like to thank Vladislava Milchevskaya, Julian Gehring and Mike Smith for helpful comments on and small contributions @@ -2015,7 +2015,7 @@ for microarrays in Bioconductor. -# Author Contributions +# 作者贡献 BK produced the first version of this article. BK implemented the main scaffold of the workflow and designed @@ -2023,5 +2023,5 @@ the workflow steps. SR implemented additional workflow steps and improved existing ones from version 1. BK and SR wrote the article. -# References +# 参考 From c8195a724f35ad0954ed0a75d2563e2fee4eb6bb Mon Sep 17 00:00:00 2001 From: Xiao W Date: Wed, 3 Jul 2019 13:56:40 -0400 Subject: [PATCH 10/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 82 +++++++++++++------------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index 8b9383d..a9a44b8 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -1902,14 +1902,14 @@ showSigOfNodes(top_GO_data, score(result_top_GO_elim), firstSigNodes = 3, # 使用 reactome 的途径富集分析 -The package `r Biocpkg("ReactomePA") ` offers the possibility to test enrichment -of specific pathways using the free, open-source, curated and peer reviewed -[Reactome](http://www.reactome.org/) pathway database [@Croft_2013; @Fabregat_2017]. -The package requires entrez identifiers, so we convert our PROBEIDs -(transcript cluster identifiers) to entrez identifiers -using the function `mapIDs` from the package `r Biocpkg("AnnotationDbi")`. -This will create a named vector that maps the PROBEIDs to the entrez ones, -with the PROBEIDs as names and the entrez ids as vector elements. +`r Biocpkg("ReactomePA") ` 包提供了 +使用免费，开源，策划和同行评审的 +[Reactome](http://www.reactome.org/) 途径数据库测试特定途径富集的可能性 [@Croft_2013; @Fabregat_2017]. +该软件包需要 entrez 标识符，因此我们使用 +`r Biocpkg("AnnotationDbi")` 包中的 +`mapIDs` 函数将我们的 PROBEID（转录集群标识符）转换为 entrez 标识符。 +这将创建一个命名向量，将 PROBEID 映射到 entrez，将 PROBEID 作为名称， +将 entrez ID 作为向量元素。 ```{r mapIDsToEntrez, dependson="createFactorOfInterestingGenes", message = FALSE} @@ -1919,11 +1919,11 @@ entrez_ids <- mapIds(hugene10sttranscriptcluster.db, column = "ENTREZID") ``` -We can now run the enrichment analysis that performs a statistical test -based on the hypergeoemtric distribution that is the same as a one sided -Fisher-test, which `r Biocpkg("topGO")` calls "Fisher-classic". -Details can be found in the vignette of the `r Biocpkg("DOSE")` -package [@Yu_2014]. +我们现在可以运行富集分析， +该分析基于超对称分布执行统计测试，该分布与单侧 Fisher 测试相同， +`r Biocpkg("topGO")` 称为 "Fisher-classic"。 +详细信息可以在 `r Biocpkg("DOSE")` 包的插图中找到 +[@Yu_2014]。 ```{r runReactomeEnrichment, fig.cap="Enriched Reactome pathways and their p–values as a bar chart.", eval = TRUE, warning=FALSE} reactome_enrich <- enrichPathway(gene = entrez_ids[DE_genes_CD], @@ -1941,52 +1941,52 @@ reactome_enrich@result$Description <- paste0(str_sub( head(as.data.frame(reactome_enrich))[1:6] ``` -Note that we trimmed pathway names to 20 characters. +请注意，我们将通道名称修剪为20个字符。 ## 可视化基于反应组的分析结果 -The top pathways can be displayed as a bar chart that displays all categories -with a p-value below the specified cutoff (Figure \@ref(fig:reactomeBar)). +顶部路径可以显示为条形图， +显示 p 值低于指定截止值的所有类别 (Figure \@ref(fig:reactomeBar)). ```{r reactomeBar, dependson="runReactomeEnrichment", fig.cap="Enriched Reactome pathways and their p–values as a bar chart."} barplot(reactome_enrich) ``` -The "enrichment map" from the package `r Biocpkg("enrichplot")` -displays the results of the enrichment analysis as -a graph, where the color represents the p-value of the pathway and the -edge-thickness (that is the line connecting two pathways) is proportional -to the number of overlapping genes between two pathways. +`r Biocpkg("enrichplot")` 包中的“富集图”将 +富集分析的结果显示为图形， +其中颜色表示路径的 p 值， +边缘厚度（即连接两个路径的线）与数量成比例 +两个途径之间的重叠基因 ```{r emapplot, dependson="runReactomeEnrichment", fig.width=6, fig.height = 7, fig.cap="Enriched Reactome pathways enrichment results as a graph."} emapplot(reactome_enrich, showCategory = 10) ``` -Again, the graph in Figure \@ref(fig:emapplot) shows pathways related -to signalling and immune response. +同样，图 \@ref(fig:emapplot) 中的图表显示了 +与信号传导和免疫应答相关的途径。 -The package `r Biocpkg("clusterProfiler") ` [@Yu_2012] can also perform these -analyses using downloaded KEGG data. Furthermore, -the package `r Biocpkg("EnrichmentBrowser")` -[@Geistlinger_2016] additionally offers network-based enrichment analysis of individual -pathways. This allows the mapping of the expression data at hand to known -regulatory interactions. +`r Biocpkg("clusterProfiler") ` [@Yu_2012] 包 +也可以使用下载的 KEGG 数据执行这些分析。 +此外，`r Biocpkg("EnrichmentBrowser")` +[@Geistlinger_2016] 包另外提供了对个体途径的基于网络的富集分析。 +这允许将手头的表达数据映射 +到已知的调节相互作用。 # 会话信息 -As the last part of this document, we call the function *sessionInfo*, -which reports the version numbers of R and all the packages used in -this session. It is good practice to always keep such a record of this -as it will help to track down what has happened in case an R script -ceases to work or gives different results because the functions have -been changed in a newer version of one of your packages. By including -it at the bottom of a script, your reports will become more reproducible. - -The session information should also *always* -be included in any emails to the -[Bioconductor support site](https://support.bioconductor.org) along -with all code used in the analysis. +作为本文档的最后一部分，我们调用函数 *sessionInfo*, +它报告 R 的版本号以及此会话中使用的所有包。 +优良作法是始终保留这样的记录， +因为它将有助于跟踪 R 脚本停止工作或给出不同结果时发生的情况， +因为功能已在您的某个软件包的较新版本中更改。 +通过将其包含在脚本的底部， +您的报告将变得更加可重复。 + +会话信息也应始终包含在任何发送到 +[Bioconductor support site](https://support.bioconductor.org) +支持网站的电子邮件中， +以及分析中使用的所有代码。 ```{r} gc() From fd067257b55401619e37a960b7d8c9ede1e8b11d Mon Sep 17 00:00:00 2001 From: Xiao Wang Date: Fri, 6 Mar 2020 19:14:31 -0800 Subject: [PATCH 11/11] Update 09rnaEndToEnd_en.Rmd --- BiocWorkflow_todo/09rnaEndToEnd_en.Rmd | 77 +++++++++++++------------- 1 file changed, 38 insertions(+), 39 deletions(-) diff --git a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd index a9a44b8..e7491bf 100644 --- a/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd +++ b/BiocWorkflow_todo/09rnaEndToEnd_en.Rmd @@ -7,7 +7,7 @@ author: affiliation: EMBL Heidelberg, Meyerhofstrasse 1, 69117 Heidelberg, Germany, steffi.reisenauer@tum.de header-includes: - \usepackage{float} -摘要: 在这篇文章里，我们将阐述用 Bioconductor 包 an end-to-end Affymetrix microarray differential expression 的工作流程。这个工作流程可以直接应用在目前 “基因” 类型的微陈列，例如 HuGene 或 MoGene 微陈列，同时也容易调整并适用于类似的平台。在这里进行分析的数据是来自典型的临床微陈列数据集，用来对比两种疾病亚型中发炎和非发炎的结肠组织。针对每种疾病，发炎和非发炎的结肠组织的差异基因表达被进行了分析。我们将从原始数据 CEL 文件开始，展示如何将其导入为 Bioconductor 的 ExpressionSet，实施质量控制和归一化和最终的差异基因表达分析，其次是一些显著性分析。 +在本文中，我们使用Bioconductor包探索了end-to-end Affymetrix microarray差异表达流程。该工作流程直接适用于当前的“基因”类型阵列，例如， HuGene或MoGene阵列，但也很容易地应用于类似的平台。本文分析的数据是典型的临床微阵列数据集，它比较两种疾病亚型中的发炎和未发炎的结肠组织。每种疾病分析发炎和未发炎的结肠组织之间的差异基因表达。我们将从原始数据CEL文件开始，展示如何将它们导入Bioconductor ExpressionSet，进行质量控制和归一化处理，最后进行差异基因表达（DE）分析和一些富集分析。 keywords: Microarray, Gene Expression bibliography: MAEndToEnd.bib output: @@ -102,27 +102,26 @@ opts_chunk$set(echo = TRUE,

-# 介绍 +# 引言 在这篇文章里我们将介绍典型的（Affymetrix）微阵列分析的完整 -工作流程。从数据导入，前期处理，差异基因表达和显著性分析都会 -被讨论。 +工作流程。下面将从数据导入，前期处理，差异基因表达和显著性分析几个方面进行讨论。 使用的数据 [@Palmieri_2015] 来自于一篇研究发炎和非发炎组 -织中差异基因表达的论文。 +织基因表达差异的文章。 测试了 14 个患有溃疡性结肠炎 (UC)的病人和 15 个患有克罗恩病 -（CD) 的病人，然后通过活检提取了每个病人发炎和未发炎的结肠粘 +（CD) 的患者，然后通过活检提取了每个病人发炎和未发炎的结肠粘 膜组织。 这个典型的临床数据集共由 58 个阵列构成。我们的目标是分析组织 间的差异表达（DE）。 -我们的结果显示了与原始论文结果的大量重叠。 +我们的结果与原始论文的结果有大量重合。 # 工作流程包的安装 -工作流程被包裹在叫做 `maEndToEnd` 的包里。 +工作流程被包含在在叫做 `maEndToEnd` 的包里。 `maEndToEnd` 包目前可以从 GitHub 获得， -目前可用的开发版本为 Bioconductor (3.8) (看这里:). +目前可用的开发版本为 Bioconductor (3.8) (详情见:). ## 从 Bioconductor 安装工作流程包 @@ -179,7 +178,7 @@ suppressPackageStartupMessages({library("maEndToEnd")}) 它们不是执行工作流所必需的。 ```{r pkgList, results="hide"} -#通常的 Bioconductor 包 +#通用的 Bioconductor 包 library(Biobase) library(oligoClasses) @@ -240,7 +239,7 @@ if (!dir.exists(raw_data_dir)) { ``` -每个 ArrayExpress 数据集都有一个总结数据集的登陆页面， +每个 ArrayExpress 数据集都有一个总结数据集的登录页面， 我们使用`r Biocpkg("ArrayExpress") ` Bioconductor 包中的 `getAE` 函数获取原始数据文件的 ftp 链接 ([Data from Palmieri et. al. on ArrayEpress](https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2967/)). @@ -287,7 +286,7 @@ SDRF 文件包含有关实验样本的基本信息，例如他们所属的实验 基因组数据可以非常复杂， 通常由许多不同的组分组成，例如， -有关实验样品的信息，测量的基因组特征的注释 +有关实验样本的信息，测量的基因组特征的注释 以及实验数据本身。 在 Bioconductor 中， 采用的方法是将这些组件存储 在一个单结构中以便于管理数据。 @@ -302,7 +301,7 @@ SDRF 文件包含有关实验样本的基本信息，例如他们所属的实验 ExpressionSet 中的数据包括： * **assayData**: 来自微阵列实验的表达数据，其中每行是微阵列探针，每列是样品标识符。 -其中每行是微阵列探针，每列是样品标识符。 +其中每行是微阵列探针，每列是样本标识符。 * **metaData** + **phenoData**: 实验中样本的描述，每行是样本标识符，每列是描述元素; 储存的是 SDRF 文件的内容。 + **featureData**: 有关芯片功能的元数据或者用于实验的技术， @@ -313,7 +312,7 @@ ExpressionSet 中的数据包括： * **experimentData**: 一种灵活的结构来描述实验。 ExpressionSet 类协调所有这些数据，因此不必 -担心细节。 但是，应该记住，`phenoData` 的 rownames （行名）必须与 assay data 的列名匹配， +担心细节。 但是，应该记住，`phenoData` 的行名必须与 assay data 的列名匹配， 而 assayData 的行名必须与 `featureData` 的行名匹配。 这在图1中说明。 \@ref(fig:sumexp). @@ -354,7 +353,7 @@ polygon(c(50,65, 65, 50),c(56,56,65,65),col="darkseagreen3",border=NA) text(57.5,60.5,"sample IDs", cex = 1) ``` -您可以使用函数 ` pData ` 和 ` fData ` 分别从 +您可以分别使用函数 ` pData ` 和 ` fData ` 从 ` ExpressionSet ` 中提取样本和要素注记。 函数 ` exprs ` 将表达式数据本身作为矩阵返回。 @@ -363,7 +362,7 @@ text(57.5,60.5,"sample IDs", cex = 1) 我们使用原始数据文件夹中的 `read.delim` 函数导入 SDRF 文件，以获取样本注释。 -样本名称在 SDRF 数据表的 Array.Data.File 列中给出，并将用作 SDRF 文件的 rownames 。 +样本名称在 SDRF 数据表的 Array.Data.File 列中给出，并将用作 SDRF 文件的行名。 我们将 SDRF 表转换为 `r Biocpkg("Biobase")` 包中 的 `AnnotatedDataFrame` ，稍后我们将需要为我们的数据 @@ -379,7 +378,7 @@ SDRF <- AnnotatedDataFrame(SDRF) 我们现在创建 Expression Set object -`raw_data`,它包含陈列数据，pheno data (表型数据) +`raw_data`,它包含 array data (陈列数据)，pheno data (表型数据) （来自 SDRF 文件） 以及所使用 的芯片注释包的信息。 @@ -429,7 +428,7 @@ head(Biobase::pData(raw_data)) 我们感兴趣的列如下： - * 个人的标识符，即列名 "Source.Name", "Characteristics.individual." + * 个体的标识符，即列名 "Source.Name", "Characteristics.individual." * 个体的疾病，即名为 "Factor.Value.disease." * 粘膜类型，即名为 "Factor.Value.phenotype." @@ -446,7 +445,7 @@ Biobase::pData(raw_data) <- Biobase::pData(raw_data)[, c("Source.Name", 初始数据导入后的第一步是数据的质量控制。 在这里，我们检查异常值并尝试查看数据是否按预期聚类， -例如 由实验条件。 表达强度值在 +例如 按照实验条件聚类。 表达强度值在 assayData 子对象 "exprs" 中，并且可以由 `exprs(raw_data)` 函数访问。 每行代表微阵列探针，即 芯片上的单个 DNA 位置， @@ -459,7 +458,7 @@ Biobase::exprs(raw_data)[1:5, 1:5] 对于质量控制，我们采用 `Biobase::exprs(raw_data) `的 log2， -因为表达数据通常以 log 分析。 +因为表达数据通常以log标度分析。 然后，我们执行主成分分析（PCA）并绘制它 \@ref(fig:qualityControlRawDataPCA)) 图中的每个点代表一 @@ -496,7 +495,7 @@ ggplot(dataGG, aes(PC1, PC2)) + 这意味着 疾病类型是基因表达差异的主要驱动因素。 这可能会妨碍我们的分析，因为我们想要分析发炎组织和非发炎组织之间的差异表达， -而不是人们患病的疾病的关系。 +而与人们患病的疾病类型无关。 我们还通过箱形图图表示探针强度， 每个单独的微阵列是有一个箱。(Figure \@ref(fig:qualityControlRawDataBox)). 注意， @@ -529,35 +528,35 @@ arrayQualityMetrics(expressionset = raw_data, intgroup = c("Factor.Value.disease.", "Factor.Value.phenotype.")) ``` -# 背景调整，校准，摘要和注释 +# 背景调整，校准，概要和注释 ## 背景调整 在初始导入和质量评估之后， 处理微阵列数据的下一步是背景调整。 这是必要的因为测量的探针强度的一部分 -是由于非特异性杂交和光学检测系统中的噪声。 +是由于非特异性杂交和光学检测系统中的噪声造成的。 因此，需要调整观察到的强度 以给出特定杂交的准确测量。 ## 跨阵列归一化（校准） 因为许多差异的来源并不明确， -所以要标准化跨阵列以便能够比较来自不同阵列杂交的测量结果。 -这些包括逆转录， -标记或杂交反应的不同效率， +所以要归一化跨阵列以便能够比较来自不同阵列杂交的测量结果。 +这些来源逆转录， +标记或杂交反应的效率不同， 阵列的物理问题， 试剂批次效应 和实验室条件。 ## 概要 -我们必须在归一化后进行汇总， +归一化后，必须进行汇总， 因为在Affymetrix平台上， 转录本由多个探针表示，即阵列上的多个位置。 对于每个基因，需要将所有探针的背景调整 和标准化强度总结为一个量， -让预计与 RNA 转录物的量成比例的量。 +它估计与 RNA 转录物的量成比例的量。 在汇总步骤之后，汇总数据可以用各种信息注释， 例如，基因符号和 ENSEMBL 基因标识符。 @@ -572,13 +571,13 @@ head(ls("package:hugene10sttranscriptcluster.db")) ``` 更多信息可从包装的参考手册中获得。 -本质上该包提供从转录本簇标识 +实际上该包提供从转录本簇标识 符到各种注释数据的映射。 ## Affymetrix 微阵列的新旧 “探针组” 传统来说，Affymetrix 阵列（所谓的3'IVT 阵列）基于探针组： -某一固定探针组是探针组的一部分， +某一组固定的探针是探针组的一部分 该探针组代表某种基因或转录物 但请注意，基因可由多个探针组表示）。 @@ -587,8 +586,8 @@ head(ls("package:hugene10sttranscriptcluster.db")) 名为 `r Biocpkg("hugene10sttranscriptcluster.db") `. -通过仅从外显子阵列中获取“好”的探针， -“基因”阵列被创建为“外显子”阵列的负担得起的版本。 +“基因”阵列比“外显子”阵列更便宜， +它仅从外显子阵列中获取“好”的探针。 最初在外显子阵列上，至少有四个探针是一个“外显子”的一部分。 通过稀释的“基因”阵列，许多探针组由三个或更少的探针组成。 这在图4 中可视化：\@ref(fig:DifferenceBetweenExonAndGeneTypeArrays) @@ -635,7 +634,7 @@ grid.rect(x=dat$x,y=dat$y, height=0.1, width=0.1, hjust=0,vjust=0,vp=vp2, 在左侧，我们看到每个外显子/探针组（即每种颜色）有大量探针： -因此，探针组/外显子水平上的摘要是合理的。 +因此，探针组/外显子水平上的总结是合理的。 然而，在基因型阵列中， 每个探针组仅包括一小部分原始探针。 因此，不建议对 “基因” 阵列使用探针组/外显子级别进行汇总， @@ -644,14 +643,14 @@ grid.rect(x=dat$x,y=dat$y, height=0.1, width=0.1, hjust=0,vjust=0,vp=vp2, 注意，此外，“基因”和“外显子”型芯片上 也不再存在指定的匹配/错配探针。 -错配探针最初用作背景校正的基础水平， -但由于更精细的背景校正技术 -不需要不匹配探针而没有普及。 +错配探针最初用作基础级别的背景校正， +它没有普及， +因为更精确的背景校正技术不需要错配探针。 ## 一步预处理寡核苷酸 -`r Biocpkg("oligo") ` 包允许我们在一个步骤里使用用于背景校正， -归一化，和汇总。 +在单个步骤里，我们用 `r Biocpkg("oligo") ` +包进行背景校正，归一化，和汇总。 这个包使用反卷积方法来背景校正，分位数归一化和 RMA (robust multichip average) 算法来汇总。 @@ -664,7 +663,7 @@ RMA (robust multichip average) 算法来汇总。 在校准和评估数据之前，我们要执行另一个质量控制程序， 即相对对数表达式（RLE）， 如 Gandolfo et al [@Gandolfo_2018]。 -我们首先执行没有事先规范化 RMA： +为此，我们首先执行 RMA, 这个 RMA 之前没有进行归一化处理。 ```{r RMAcalibrationForRLE, eval=TRUE} palmieri_eset <- oligo::rma(raw_data, target = "core", normalize = FALSE)