RNA-seq工作流程：基因水平的探索性分析和差异表达
RNA-seq工作流程：基因水平的探索性分析和差异表达
迈克尔·爱1，西蒙·安德斯3，弗拉迪斯拉夫·金4和沃尔夫冈·胡贝尔4
1美国北卡罗莱纳州教堂山UNC-Chapel Hill生物统计学系
2美国北卡罗莱纳州教堂山UNC-Chapel Hill遗传学系
3德国海德堡海德堡大学分子生物学研究所
4欧洲分子生物学实验室（ EMBL），德国海德堡
2019年10月16日
抽象
在这里，我们介绍了使用Bioconductor软件包进行的端到端基因水平RNA-seq差异表达工作流程。我们将从FASTQ文件开始，展示如何将它们量化为参考转录本，并准备基因水平计数数据集以进行下游分析。我们将进行探索性数据分析（EDA）以进行质量评估，并探索样品之间的关系，进行差异基因表达分析，并直观地观察结果。
内容
- 1引言
- 1.1实验数据
- 2准备量化输入到DESeq2
- 2.1笔录量化和tximport / tximeta
- 2.2用鲑鱼定量
- 2.3用tximeta读数据
- 2.4 DESeq2导入功能
- 2.5总结实验
- 2.6分支点
- 3的DESeqDataSet对象，样品信息和设计配方
- 3.1从总结实验开始
- 3.2从计数矩阵开始
- 4探索性分析和可视化
- 4.1预过滤数据集
- 4.2方差稳定化变换和对数
- 4.3采样距离
- 4.4 PCA图
- 4.5使用广义PCA的PCA图
- 4.6 MDS图
- 5差异表达分析
- 5.1运行差异表达式管道
- 5.2建立结果表
- 5.3其他比较
- 5.4多重测试
- 6绘图结果
- 6.1计数图
- 6.2 MA图
- 6.3基因聚类
- 6.4独立过滤
- 6.5独立假设加权
- 7注释和导出结果
- 7.1导出结果
- 7.2绘制基因组空间的倍数变化
- 8删除隐藏的批处理效果
- 8.1将SVA与DESeq2一起使用
- 8.2将RUV与DESeq2一起使用
- 9时程实验
- 10会议信息
- 参考文献
R版本：R版本4.0.3（2020-10-10）
生物导体版本：3.12
包装：1.14.0
1个介绍
Bioconductor具有许多支持高通量序列数据分析的软件包，包括RNA测序（RNA-seq）。我们将在此工作流程中使用的软件包包括由Bioconductor核心团队维护的核心软件包，用于处理基因注释（基因组中的基因和转录本位置以及基因ID查找）。我们还将使用贡献包进行统计分析和序列数据的可视化。通过每6个月的预定发行，Bioconductor项目可确保发行中的所有程序包都可以和谐地协同工作（因此称为“导体”隐喻）。此工作流程中使用的软件包已加载了 库功能，可以按照 Bioconductor软件包安装说明进行安装。
- F1000Research上提供了此工作流程的已发布版本，包括审阅者报告和评论。您正在阅读的版本与此版本有所不同，主要是因为我们现在提供用于执行快速转录定量的代码， 然后导入R / Bioconductor以执行基因水平分析。
- 另一个涉及差异转录本使用 （DTU）的生物导体工作流程是 rnaseqDTU工作流程，其发布的版本也可以在 F1000Research获得。
- 如果您对此工作流程或任何Bioconductor软件有疑问，请将其发布到 Bioconductor支持站点。如果问题与特定软件包有关，则可以用软件包的名称标记帖子，或者对于有关工作流程的一般问题，请使用标记帖子
rnaseqgene。请注意 发布指南 ，以为支持站点制定最佳问题。
1.1实验数据
该工作流程中使用的数据存储在气道包装中，该数据 总结了RNA-seq实验，其中用地塞米松（一种具有抗炎作用的合成糖皮质激素类固醇）处理了气道平滑肌细胞 （Himes等，2014）。糖皮质激素被例如哮喘患者用来减轻气道炎症。在该实验中，用1微摩尔地塞米松处理了四个主要的人气道平滑肌细胞系18小时。对于这四个细胞系中的每一个，我们都有一个处理过的样品和一个未处理过的样品。有关实验的更多说明，请参见 PubMed条目24926665 ，有关原始数据，请参见 GEO条目GSE52778。
2准备量化输入到DESeq2
作为输入，以计数为基础的统计方法，如 DESeq2 （爱，胡贝尔和安德斯2014） ， 磨边机 （罗宾逊，麦卡锡和2009年史密斯）， LIMMA与VOOM法（Law等人。2014） ， DSS （吴，Wang和Wu 2013）， EBSeq （Leng等人，2013）和 baySeq （Hardcastle和Kelly，2010）期望输入数据例如是从RNA-seq或其他高通量测序实验中获得的，其形式为非标准化计数矩阵。第i行和第j行中的值矩阵的第4列告诉您可以将样本j中的基因i分配给基因i多少个读数（或片段，用于配对末端RNA-seq）。类似地，对于其他类型的测定，基质的行可能对应于例如结合区（具有ChIP-Seq）或肽序列（具有定量质谱）。
矩阵中的值应为测序读数/片段的计数或估计计数。这对于保持DESeq2的统计模型非常重要，因为只有计数才能正确评估测量精度。重要的是，永远不要提供针对测序深度/库大小进行了预归一化的计数，因为当将统计模型应用于非归一化计数时，该统计模型最为强大，并且该统计模型旨在内部解决库大小的差异。
2.1笔录量化和tximport / tximeta
该工作流程的先前版本（包括已发布的版本）演示了如何将读段与基因组对齐，然后计算与基因模型一致的读段数量。现在，我们建议使用更快的替代途径进行基因组比对和读数计数。该工作流程将演示如何导入转录本水平的定量数据，并通过tximport或 tximeta聚合到基因水平。转录物定量方法，例如 Salmon （Patro等人2017）， kallisto （Bray等人2016）或 RSEM （Li and Dewey 2011）对参考转录本进行读图的映射或比对，输出每个转录本的估计计数以及有效的转录本长度（总结了偏倚效应）。运行这些工具之一后，可以使用 tximport （Soneson，Love和Robinson 2015）或 tximeta （Love等人2020） 软件包来组装估计的计数和补偿矩阵，以与Bioconductor差异基因表达软件包一起使用。如下所示。
可在此处找到有关如何使用Salmon软件定量转录本丰度的教程。我们建议使用标记 来估计RNA-seq数据中常见的系统偏差的校正因子（Love，Hogenesch和Irizarry 2016; Patro等人2017），除非您确定您的数据不包含此类偏差。--gcBias
将转录本丰度定量词与tximport / tximeta结合使用以产生基因水平的计数矩阵并归一化偏移量的优点是：（1）这种方法可以校正样本中基因长度的任何潜在变化（例如，来自不同的同工型）（ Trapnell et al.2013）；（2）与基于对齐的方法相比，这些方法中的某些方法实质上更快，并且所需的内存和磁盘使用量更少；（3）可以避免丢弃那些可以与具有同源序列的多个基因比对的片段（Robert and Watson 2015）。请注意，转录本丰度量词会跳过存储读取比对的大文件的生成，而不是生成存储每个转录本的估计丰度，计数和有效长度的较小文件。有关更多详细信息，请参见描述此方法的手稿 （Soneson，Love和Robinson 2015），以及有关软件详细信息的tximport软件包插图。
tximeta （Love等人，2020年） 扩展了tximport，提供了相同的功能，以及为常用转录组（GENCODE，Ensembl，人和小鼠的RefSeq）自动添加注释元数据的额外好处。有关 更多详细信息，请参见 tximeta小插图包装小插图。tximeta产生一个 SummarizedExperiment，可以使用该函数轻松将其加载到DESeq2中DESeqDataSet，这将在下面进行演示。我们还将讨论是否可以使用tximport， tximeta， htseq来 输入DESeq2的各种可能的输入（Anders，Pyl和Huber 2015） ，或预先计算的计数矩阵。
2.2用鲑鱼定量
如上所述，可以在此处找到有关如何使用Salmon的简短教程，因此，我们将提供用于量化此工作流中使用的文件的代码。使用Snakemake工作流管理系统，可以在集群上方便地运行Salmon （Köster和Rahmann，2012年）。
以下Snakemake文件用于量化从SRA下载的八个样本（SRR标识符是运行标识符，对于这八个样本，每个样本只有一个运行）。
DATASETS = ["SRR1039508","SRR1039509","SRR1039512","SRR1039513","SRR1039516","SRR1039517","SRR1039520","SRR1039521"]SALMON = "/path/to/salmon_0.14.1/bin/salmon"rule all:input: expand("quants/{dataset}/quant.sf", dataset=DATASETS)rule salmon_quant:input:r1 = "fastq/{sample}_1.fastq.gz",r2 = "fastq/{sample}_2.fastq.gz",index = "/path/to/gencode.v29_salmon_0.14.1"output:"quants/{sample}/quant.sf"params:dir = "quants/{sample}"shell:"{SALMON} quant -i {input.index} -l A -p 6 --validateMappings \--gcBias --numGibbsSamples 20 -o {params.dir} \-1 {input.r1} -2 {input.r2}"最后一行是运行Salmon的关键行。它表示要使用特定索引进行量化，并使用6个线程通过自动映射类型检测，并使用验证映射设置（Salmon版本中的默认设置） ≥≥0.99），具有GC偏差校正，并写出20 Gibbs样本（这是可选的）。最后三个参数指定输出目录和两个配对的读取文件。
上面的Snakemake文件要求在处创建索引 /path/to/gencode.vVV_salmon_X.Y.Z，其中VV和X，Y，Z应该有助于指定参考成绩单和Salmon的发布。对于人类和小鼠参考转录本，我们建议使用 GENCODE （Frankish et al.2018 ）。
的鲑鱼指数可以很容易地用下面的命令创建：
salmon index -t transcripts.fa.gz -i name_of_index如果成绩单是从GENCODE下载的，则建议使用与以下命令类似的命令（它仅有助于从成绩单名称中剥离额外的信息）：
salmon index --gencode -t gencode.v29.transcripts.fa.gz \-i gencode.v29_salmon_X.Y.Z然后，可以使用上述各种Snakemake文件 以各种方式执行Snakemake，包括将多个作业提交到计算群集或云中。上面的Snakemake文件是在具有SLURM调度的群集上执行的，并将以下行提交到群集中的单独作业中：
snakemake -j 4 --latency-wait 30 --cluster "sbatch -N 1 -n 6"2.3用tximeta读数据
在工作流程的后面，我们将加载一个对象，其中包含所有八个样本在基因级别的定量数据。但是，气道软件包还包含Salmon输出的两个定量目录，以演示如何将这些数据读入R / Bioconductor。为了使数据包更小，quant.sf量化目录中的 文件已压缩，因此在您看到的下面quant.sf.gz，您可能会 quant.sf在自己的计算机上使用。
在演示了使用tximeta进行导入之后，我们将加载与所有样品和所有数据相对应的完整计数矩阵，该矩阵已经在同一包中提供，并将继续对该完整数据对象进行分析。
我们首先使用示例数据加载数据包：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"airway"</span>)</code></span></span>R函数system.file可用于查找软件包中文件在计算机上的安装位置。在这里，我们要求extdata目录的完整路径，其中R包在其中存储外部数据，这是气道包的一部分。
<span style="color:#333333"><span style="color:#333333"><code>dir <- system.file(<span style="color:#dd1144">"extdata"</span>, package=<span style="color:#dd1144">"airway"</span>, mustWork=<span style="color:#990073">TRUE</span>)</code></span></span>在此目录中，我们找到许多文件，其中包括在此工作流程的先前版本中使用的八个BAM文件，这些文件证明了对齐和计数。我们将重点关注目录中的两个目录quants，其中包含Salmon在两个文件上的输出。
<span style="color:#333333"><span style="color:#333333"><code>list.files(dir)</code></span></span>## [1] "GSE52778_series_matrix.txt" "Homo_sapiens.GRCh37.75_subset.gtf"
## [3] "SRR1039508_subset.bam" "SRR1039509_subset.bam"
## [5] "SRR1039512_subset.bam" "SRR1039513_subset.bam"
## [7] "SRR1039516_subset.bam" "SRR1039517_subset.bam"
## [9] "SRR1039520_subset.bam" "SRR1039521_subset.bam"
## [11] "SraRunInfo_SRP033351.csv" "quants"
## [13] "sample_table.csv"<span style="color:#333333"><span style="color:#333333"><code>list.files(file.path(dir, <span style="color:#dd1144">"quants"</span>))</code></span></span>## [1] "SRR1039508" "SRR1039509"通常，我们有一个表格，其中包含每个样品的详细信息，该表将样品链接到关联的FASTQ和Salmon目录。对于您自己的项目，您可以使用文本编辑器或电子表格软件（例如Excel）创建这样的逗号分隔值（CSV）文件。
我们使用read.csv加载这样的CSV文件：
<span style="color:#333333"><span style="color:#333333"><code>csvfile <- file.path(dir, <span style="color:#dd1144">"sample_table.csv"</span>)
coldata <- read.csv(csvfile, row.names=<span style="color:#009999">1</span>, stringsAsFactors=<span style="color:#990073">FALSE</span>)
coldata</code></span></span>## SampleName cell dex albut Run avgLength Experiment
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358
## Sample BioSample
## SRR1039508 SRS508568 SAMN02422669
## SRR1039509 SRS508567 SAMN02422675
## SRR1039512 SRS508571 SAMN02422678
## SRR1039513 SRS508572 SAMN02422670
## SRR1039516 SRS508575 SAMN02422682
## SRR1039517 SRS508576 SAMN02422673
## SRR1039520 SRS508579 SAMN02422683
## SRR1039521 SRS508580 SAMN02422677为了演示将鲑鱼定量数据加载到R中，我们将仅使用气道 包装中提供的两个样本。我们创建一个名为names的列和一个名为的列 files：
<span style="color:#333333"><span style="color:#333333"><code>coldata <- coldata[<span style="color:#009999">1</span>:<span style="color:#009999">2</span>,]
coldata$names <- coldata$Run
coldata$files <- file.path(dir, <span style="color:#dd1144">"quants"</span>, coldata$names, <span style="color:#dd1144">"quant.sf.gz"</span>)
file.exists(coldata$files)</code></span></span>## [1] TRUE TRUE现在，我们加载tximeta软件包并运行其主要功能：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"tximeta"</span>)
se <- tximeta(coldata)</code></span></span>## importing quantifications## reading in files with read_tsv## 1 2
## found matching transcriptome:
## [ GENCODE - Homo sapiens - release 29 ]
## useHub=TRUE: checking for TxDb via 'AnnotationHub'
## using temporary cache /tmp/RtmpE5fMSc/BiocFileCache
## snapshotDate(): 2020-10-26
## found matching TxDb via 'AnnotationHub'
## downloading 1 resources
## retrieving 1 resource
## loading from cache
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
## generating transcript ranges如果参考转录组校验和被tximeta识别（tximeta插图中的详细信息），并且如果我们的互联网连接正常，则tximeta将从各种来源找到并下载相关的注释数据。一些细节：注释数据仅下载和解析一次，随后将根据需要使用元数据的本地缓存版本（如果您第二次加载针对相同参考成绩单量化的数据）。同样，第一次使用tximeta时，它会要求您批准默认的缓存位置（遵循其他R和Bioconductor软件包使用的缓存位置的范例）。您以后可以随时更改此位置。
se在下一节中，我们将讨论对象的结构是什么，但我们首先可以考虑尺寸。请注意， tximeta会在成绩单级别导入数据。
<span style="color:#333333"><span style="color:#333333"><code>dim(se)</code></span></span>## [1] 205870 2<span style="color:#333333"><span style="color:#333333"><code>head(rownames(se))</code></span></span>## [1] "ENST00000456328.2" "ENST00000450305.2" "ENST00000488147.1"
## [4] "ENST00000619216.1" "ENST00000473358.1" "ENST00000469289.1"由于此工作流程涉及基因水平的分析，因此我们现在将转录水平的定量总结到基因水平（内部使用tximport 中的方法（Soneson，Love和Robinson，2015年））。基于存储在se对象中的元数据，会自动创建正确的转录本-基因映射表 。
<span style="color:#333333"><span style="color:#333333"><code>gse <- summarizeToGene(se)</code></span></span>## loading existing TxDb created: 2020-10-30 19:30:40## obtaining transcript-to-gene mapping from database## generating gene ranges## summarizing abundance## summarizing counts## summarizing length现在我们可以检查尺寸是否减小，并且行ID现在是基因ID：
<span style="color:#333333"><span style="color:#333333"><code>dim(gse)</code></span></span>## [1] 58294 2<span style="color:#333333"><span style="color:#333333"><code>head(rownames(gse))</code></span></span>## [1] "ENSG00000000003.14" "ENSG00000000005.5" "ENSG00000000419.12"
## [4] "ENSG00000000457.13" "ENSG00000000460.16" "ENSG00000000938.12"2.4DESeq2导入功能
虽然上述部分中描述的使用鲑鱼和tximeta，有许多可能的输入，以DESeq2，其中的每一个具有自己专用的导入功能。可以使用以下工具来生成或编译用于DESeq2的计数数据： tximport （Soneson，Love和Robinson 2015）， tximeta （Love等人2020）， htseq-count （Anders，Pyl和Huber 2015）， featureCounts （Liao，Smyth，and Shi 2014）， 总结了重叠 （Lawrence et al.2013）。
| 功能 | 包 | 构架 | 输出 | DESeq2输入功能 |
|---|---|---|---|---|
| tximport | tximport | R /生物导体 | 矩阵列表 | DESeqDataSetFromTximport |
| 紫杉醇 | 紫杉醇 | R /生物导体 | 总结实验 | DESeqDataSet |
| htseq计数 | HTSeq | 蟒蛇 | 档案 | DESeqDataSetFromHTSeq |
| featureCounts | 重新读取 | R /生物导体 | 矩阵 | DESeqDataSetFromMatrix |
| summaryOverlaps | 基因组比对 | R /生物导体 | 总结实验 | DESeqDataSet |
接下来，我们将描述由tximeta创建的对象类，该对象类另存为se和gse，以及如何从中创建DESeqDataSet 对象以与DESeq2一起使用（上面的其他函数也创建了DESeqDataSet）。
2.5总结实验
SummarizedExperiment对象的组成部分。的assay（粉块）包含计数的矩阵，所述rowRanges（蓝色块）包含有关的基因组范围和信息colData（绿色方框）包含有关样品的信息。每个块中突出显示的行代表第一行（请注意，第一行colData 与的第一列并排assay）。
所述SummarizedExperiment容器在上述图图解和在最新Bioconductor的纸讨论（Huber等人。2015） 。在我们的案例中，tximeta创建了一个gse具有三个矩阵的对象：“计数”-每个基因和样品的估计片段计数，“丰度” -TPM中的估计转录本丰度，“长度”-包含变化的有效基因长度由于偏见和抄本的使用导致长度的增加。可以使用assayNames检查化验的名称，并将化验本身存储为assays（矩阵列表）。列表中的第一个矩阵可以通过拔出assay。在rowRanges我们的对象是农庄基因（从所有转录本的最左侧位置到所有转录本的最右侧位置）。所述的组成部件SummarizedExperiment都具有相同名称的R函数访问：assay（或者assays），rowRanges和colData。
现在，我们将加载气道包装中提供的与所有样本和所有数据相对应的全计数矩阵，并将继续使用全数据对象进行分析。我们可以研究这种SummarizedExperiment通过在看矩阵对象assays插槽，有关样本的表型数据 colData插槽，以及有关在基因数据rowRanges插槽。
<span style="color:#333333"><span style="color:#333333"><code>data(gse)
gse</code></span></span>## class: RangedSummarizedExperiment
## dim: 58294 8
## metadata(6): tximetaInfo quantInfo ... txomeInfo txdbInfo
## assays(3): counts abundance length
## rownames(58294): ENSG00000000003.14 ENSG00000000005.5 ...
## ENSG00000285993.1 ENSG00000285994.1
## rowData names(1): gene_id
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(3): names donor condition计数是第一个矩阵，因此我们可以使用以下命令检查它们 assay：
<span style="color:#333333"><span style="color:#333333"><code>assayNames(gse)</code></span></span>## [1] "counts" "abundance" "length"<span style="color:#333333"><span style="color:#333333"><code>head(assay(gse), <span style="color:#009999">3</span>)</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14 708.164 467.962 900.992 424.368 1188.295
## ENSG00000000005.5 0.000 0.000 0.000 0.000 0.000
## ENSG00000000419.12 455.000 510.000 604.000 352.000 583.000
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14 1090.668 805.929 599.337
## ENSG00000000005.5 0.000 0.000 0.000
## ENSG00000000419.12 773.999 409.999 499.000<span style="color:#333333"><span style="color:#333333"><code>colSums(assay(gse))</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## 21100805 19298584 26145537 15688246 25268618 31891456 19683767
## SRR1039521
## 21813903的rowRanges，打印时，示出了用于前五个及最后五个基因的范围：
<span style="color:#333333"><span style="color:#333333"><code>rowRanges(gse)</code></span></span>## GRanges object with 58294 ranges and 1 metadata column:
## seqnames ranges strand | gene_id
## <Rle> <IRanges> <Rle> | <character>
## ENSG00000000003.14 chrX 100627109-100639991 - | ENSG00000000003.14
## ENSG00000000005.5 chrX 100584802-100599885 + | ENSG00000000005.5
## ENSG00000000419.12 chr20 50934867-50958555 - | ENSG00000000419.12
## ENSG00000000457.13 chr1 169849631-169894267 - | ENSG00000000457.13
## ENSG00000000460.16 chr1 169662007-169854080 + | ENSG00000000460.16
## ... ... ... ... . ...
## ENSG00000285990.1 chr14 19244904-19269380 - | ENSG00000285990.1
## ENSG00000285991.1 chr6 149817937-149896011 - | ENSG00000285991.1
## ENSG00000285992.1 chr8 47129262-47132628 + | ENSG00000285992.1
## ENSG00000285993.1 chr18 46409197-46410645 - | ENSG00000285993.1
## ENSG00000285994.1 chr10 12563151-12567351 + | ENSG00000285994.1
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genome在rowRanges还包含关于在该序列（在我们的情况下，染色体）的元数据seqinfo时隙：
<span style="color:#333333"><span style="color:#333333"><code>seqinfo(rowRanges(gse))</code></span></span>## Seqinfo object with 25 sequences (1 circular) from hg38 genome:
## seqnames seqlengths isCircular genome
## chr1 248956422 FALSE hg38
## chr2 242193529 FALSE hg38
## chr3 198295559 FALSE hg38
## chr4 190214555 FALSE hg38
## chr5 181538259 FALSE hg38
## ... ... ... ...
## chr21 46709983 FALSE hg38
## chr22 50818468 FALSE hg38
## chrX 156040895 FALSE hg38
## chrY 57227415 FALSE hg38
## chrM 16569 TRUE hg38在colData为SummarizedExperiment反映了data.frame 这是提供给tximeta功能导入定量数据。在这里我们可以看到有几列显示样品名称，供体ID和治疗条件（用地塞米松治疗或未治疗）。
<span style="color:#333333"><span style="color:#333333"><code>colData(gse)</code></span></span>## DataFrame with 8 rows and 3 columns
## names donor condition
## <factor> <factor> <factor>
## SRR1039508 SRR1039508 N61311 Untreated
## SRR1039509 SRR1039509 N61311 Dexamethasone
## SRR1039512 SRR1039512 N052611 Untreated
## SRR1039513 SRR1039513 N052611 Dexamethasone
## SRR1039516 SRR1039516 N080611 Untreated
## SRR1039517 SRR1039517 N080611 Dexamethasone
## SRR1039520 SRR1039520 N061011 Untreated
## SRR1039521 SRR1039521 N061011 Dexamethasone2.6分支点
至此，我们已经计算了在指定基因模型中与基因重叠的片段。这是一个分支点，在这里我们可以使用各种Bioconductor软件包进行计数数据的探索和差异表达，包括 edgeR （Robinson，McCarthy和Smyth 2009）， limma和voom方法（Law等人，2014）， DSS （Wu，Wang和Wu，2013年）， EBSeq （Leng等人，2013年）和 baySeq （Hardcastle和Kelly，2010年）。 Schurch等。（2016） 比较表现 大量生物复制品用于RNA-seq的不同统计方法的比较，可以帮助用户决定使用哪种工具有意义，以及要获得一定的敏感性需要多少个生物复制品。我们将继续使用DESeq2 （Love，Huber和Anders 2014）。该SummarizedExperiment对象是所有我们需要开始我们的分析。在以下部分中，我们将展示如何使用它来创建DESeq2使用的数据对象。
3的DESeqDataSet对象，样品信息和设计配方
生物导体软件包通常定义并使用自定义类来存储数据，以确保一致地提供所有需要的数据插槽并满足要求。此外，Bioconductor具有通用数据类（例如 SummarizedExperiment），可用于在包装之间移动数据。此外，核心Bioconductor类提供有用的功能：例如，对SummarizedExperiment的行或列进行子集设置或重新排序会自动对关联的rowRanges和colData进行子集或重新排序，这有助于防止意外的样本交换，否则可能会导致虚假结果。随着SummarizedExperiment 这一切都在后台进行。
在DESeq2中，自定义类称为DESeqDataSet。它建立在SummarizedExperiment类的顶部，并且很容易将 SummarizedExperiment对象转换为DESeqDataSet对象，如下所示。两个主要区别之一是assay改为使用counts访问器函数访问该插槽，而 DESeqDataSet类强制此矩阵中的值是非负整数。
第二个区别是DESeqDataSet具有关联的 设计公式。实验设计是在分析开始时指定的，因为它将告知许多DESeq2 函数如何在分析中处理样品（一个例外是大小因子估计，即针对不同库大小的调整，但不会取决于设计公式）。设计公式说明样品信息表（colData）中的哪些列指定了实验设计，以及在分析中应如何使用这些因素。
首先，让我们检查的列colData的gse。我们可以$直接在SummarizedExperiment或DESeqDataSet上 使用来查看每个列。
<span style="color:#333333"><span style="color:#333333"><code>gse$donor</code></span></span>## [1] N61311 N61311 N052611 N052611 N080611 N080611 N061011 N061011
## Levels: N052611 N061011 N080611 N61311<span style="color:#333333"><span style="color:#333333"><code>gse$condition</code></span></span>## [1] Untreated Dexamethasone Untreated Dexamethasone Untreated
## [6] Dexamethasone Untreated Dexamethasone
## Levels: Untreated Dexamethasone如果需要，我们可以重命名变量。让我们cell用来表示供体细胞系，并dex表示治疗条件。
<span style="color:#333333"><span style="color:#333333"><code>gse$cell <- gse$donor
gse$dex <- gse$condition</code></span></span>我们还可以更改级别的名称。当一个人重命名关卡而不改变顺序时，这一点至关重要。在这里，我们将重命名 "Untreated"为"untrt"和"Dexamethasone"为"trt"：
<span style="color:#333333"><span style="color:#333333"><code>levels(gse$dex)</code></span></span>## [1] "Untreated" "Dexamethasone"<span style="color:#333333"><span style="color:#333333"><code><span style="color:#999988"><em># when renaming levels, the order must be preserved!</em></span>
levels(gse$dex) <- c(<span style="color:#dd1144">"untrt"</span>, <span style="color:#dd1144">"trt"</span>)</code></span></span>最简单的差异表达设计公式为 ~ condition，其中condition的一列colData(dds)指定了样本属于两个（或更多）组中的哪一个。对于气道实验，我们将指定~ cell + dex 要测试地塞米松（dex）控制不同细胞系（cell）效果的含义。
注意：在R中，因子的第一水平最好是参考水平（例如，对照或未处理的样品）。在这种情况下，当colData组装桌子时，已经将未经处理的样品设置为参考，但是如果不是这种情况，我们可以使用如下所示的整平。虽然levels(...) <-上面只是为了重命名与关卡相关的字符串，但 relevel是一个非常不同的功能，它决定如何对变量进行编码以及如何计算对比度。对于两组比较，使用重新电平来更改参考电平会翻转与两组之间的对比度相关的系数的符号。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"magrittr"</span>)
gse$dex %<>% relevel(<span style="color:#dd1144">"untrt"</span>)
gse$dex</code></span></span>## [1] untrt trt untrt trt untrt trt untrt trt
## Levels: untrt trt%<>%是magrittr软件包中的复合赋值管道运算符，上面的代码行是一种简洁的说法：
<span style="color:#333333"><span style="color:#333333"><code>gse$dex <- relevel(gse$dex, <span style="color:#dd1144">"untrt"</span>)</code></span></span>对于运行DESeq2模型，您可以使用R的公式表示法来表示任何固定效果的实验设计。请注意，DESeq2使用与例如 碱基R的lm函数相同的公式表示法。如果研究目的是确定哪些基因对各组的治疗效果不同，则可以包括相互作用项并使用设计进行测试如 ~ group + treatment + group:treatment。有关?results更多示例，请参见手册页 。在下面的时间过程示例中，我们将展示如何使用交互作用项来测试随时间变化的特定条件。
在以下各节中，我们将从两个起点演示DESeqDataSet的构造 ：
- 来自SummarizedExperiment对象
- 从计数矩阵和样本信息表中
有关使用HTSeq Python软件包进行读取计数的完整示例，请参见pasilla小插图。用于生成的一例DESeqDataSet从由产生的文件 htseq计数，请参阅DESeq2晕影。
3.1从总结实验开始
同样，我们可以快速检查Salmon可以将其映射到基因的数百万个片段（回合的第二个参数表明要保留多少个小数点）。
<span style="color:#333333"><span style="color:#333333"><code>round( colSums(assay(gse)) / <span style="color:#009999">1e6</span>, <span style="color:#009999">1</span> )</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## 21.1 19.3 26.1 15.7 25.3 31.9 19.7
## SRR1039521
## 21.8一旦有了完全注释的SummarizedExperiment对象，就可以从该对象构造一个DESeqDataSet对象，该对象将构成分析的起点。我们为分析添加适当的设计：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"DESeq2"</span>)</code></span></span><span style="color:#333333"><span style="color:#333333"><code>dds <- DESeqDataSet(gse, design = ~ cell + dex)</code></span></span>3.2从计数矩阵开始
在本节中，我们将假设只有一个计数矩阵和一个样本信息表，展示如何构建DESeqDataSet。
注意：如果您准备了SummarizedExperiment，则应跳过此部分。虽然前面的部分将被用于构建 DESeqDataSet从SummarizedExperiment，在这里我们首先提取了单个对象（数矩阵和样本信息） SummarizedExperiment为了建立它备份到一个新的目标-仅用于演示目的。在实践中，计数矩阵要么被在从文件中读取或可能通过类似的R函数生成featureCounts从 Rsubread包（辽，史密斯和施2014） 。
可以使用访问器函数访问SummarizedExperiment对象中的信息。例如，要查看实际数据，即此处的片段计数，我们可以使用化验函数。（head 函数将输出限制为前几行。）
<span style="color:#333333"><span style="color:#333333"><code>countdata <- round(assays(gse)[[<span style="color:#dd1144">"counts"</span>]])
head(countdata, <span style="color:#009999">3</span>)</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14 708 468 901 424 1188
## ENSG00000000005.5 0 0 0 0 0
## ENSG00000000419.12 455 510 604 352 583
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14 1091 806 599
## ENSG00000000005.5 0 0 0
## ENSG00000000419.12 774 410 499在该计数矩阵中，每一行代表一个基因，每一列代表一个测序的RNA文库，这些值给出了Salmon概率分配给每个文库中各个基因的片段的估计计数。我们还具有每个样本的信息（计数矩阵的列）。如果您已经以其他方式导入了计数数据，例如加载了预先计算的计数矩阵，那么手动检查计数矩阵的列是否与样本信息表的行相对应非常重要。
<span style="color:#333333"><span style="color:#333333"><code>coldata <- colData(gse)</code></span></span>现在，我们拥有所有要素来以适合分析的形式准备数据对象，即：
countdata：带有碎片计数的表coldata：带有样本信息的表格
现在要从计数矩阵和样本信息表构造DESeqDataSet对象，我们使用：
<span style="color:#333333"><span style="color:#333333"><code>ddsMat <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ cell + dex)</code></span></span>我们将继续从SummarizedExperiment部分生成的对象 。
4探索性分析和可视化
此工作流程中有两个单独的路径。我们将首先看到的是计数转换， 以便直观地探索样本关系。在第二部分中，我们将返回原始原始计数以进行统计测试。这是至关重要的，因为统计测试方法依赖原始计数数据（未经缩放或转换）来计算测量精度。
4.1预过滤数据集
Our count matrix with our DESeqDataSet contains many rows with only zeros, and additionally many rows with only a few fragments total. In order to reduce the size of the object, and to increase the speed of our functions, we can remove the rows that have no or nearly no information about the amount of gene expression. Here we apply the most minimal filtering rule: removing rows of the DESeqDataSet that have no counts, or only a single count across all samples. Additional weighting/filtering to improve power is applied at a later step in the workflow.
<span style="color:#333333"><span style="color:#333333"><code>nrow(dds)</code></span></span>## [1] 58294<span style="color:#333333"><span style="color:#333333"><code>keep <- rowSums(counts(dds)) > <span style="color:#009999">1</span>
dds <- dds[keep,]
nrow(dds)</code></span></span>## [1] 31604对于某些数据集，执行附加过滤可能很有意义。例如，可以指定至少3个样本的计数为10或更高。关于样本数量的一项建议将设置为最小的组大小。这样的规则可以通过创建一个逻辑向量并按dds上述子集指定。这是我们可以使用的另一条规则的示例（此处未用于过滤）：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#999988"><em># at least 3 samples with a count of 10 or higher</em></span>
keep <- rowSums(counts(dds) >= <span style="color:#009999">10</span>) >= <span style="color:#009999">3</span></code></span></span>4.2方差稳定化变换和rlog
探索性多维数据分析的许多常用统计方法，例如聚类和主成分分析（PCA），最适合通常在不同平均值范围内具有相同方差范围的数据。当方差的预期量为约跨越不同的平均值是相同的，该数据被认为是同方差。但是，对于RNA-seq计数，预期方差随平均值而增加。例如，如果人们直接在计数或归一化计数的矩阵上执行PCA（例如，校正测序深度的差异），则结果图通常主要取决于具有最高序列的基因。计数是因为它们在样本之间显示出最大的绝对差异。避免这种情况的一种简单且经常使用的策略是对归一化计数值加上对数1的对数；但是，根据伪计数的选择，现在计数最低的基因将对结果图产生很大的干扰，因为取小计数的对数实际上会夸大其方差。我们可以使用一些模拟数据来快速显示计数的这种属性（此处，泊松计数的Lambda范围为0.1到100）。我们针对均值绘制每行（基因）的标准偏差：
<span style="color:#333333"><span style="color:#333333"><code>lambda <- <span style="color:#009999">10</span>^seq(from = -<span style="color:#009999">1</span>, to = <span style="color:#009999">2</span>, length = <span style="color:#009999">1000</span>)
cts <- matrix(rpois(<span style="color:#009999">1000</span>*<span style="color:#009999">100</span>, lambda), ncol = <span style="color:#009999">100</span>)
<span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"vsn"</span>)
meanSdPlot(cts, ranks = <span style="color:#990073">FALSE</span>)</code></span></span>
对于对数转换的计数：
<span style="color:#333333"><span style="color:#333333"><code>log.cts.one <- log2(cts + <span style="color:#009999">1</span>)
meanSdPlot(log.cts.one, ranks = <span style="color:#990073">FALSE</span>)</code></span></span>
当数值接近0时，具有较小伪计数的对数会放大差异。具有低信噪比的低计数基因会过度影响样本-样本距离和PCA图。
作为解决方案，DESeq2为计数数据提供了两种转换，可稳定均值上的方差：具有离散均值趋势的负二项式数据的方差稳定转换（VST） （在vst函数中实现）和的正则化对数转换或rlog可 （爱，胡伯，和安德斯2014） 。
对于具有高计数的基因，VST和rlog都将给出与归一化计数的普通log2转换相似的结果。但是，对于计数较低的基因，其值会缩小到中间值。该VST或rlog可变换数据然后成为大致同方差（在更平坦的趋势meanSdPlot），可直接用于计算样本之间的距离，使得PCA图，或作为输入到与同方差数据进行最佳下游方法。
选择哪种转换？ 与rlog相比，VST的计算速度快得多，并且对高计数值异常值的敏感性较低。rlog倾向于在较小的数据集（n <30）上很好地工作，当样本之间的测序深度范围很广时（数量级差），可能会胜过VST。因此，我们建议将VST用于大中型数据集（n> 30）。您可以执行两个转换，并比较meanSdPlot生成的或PCA图，如下所述。
注意的是，提供应用程序通过DESeq2提供的两个转变等比差测试。对于差异测试，我们建议将 DESeq函数应用于原始计数，如本工作流程中稍后所述，该函数还考虑了色散估计步骤中计数方差对平均值的依赖性。
无论VST和rlog可回报DESeqTransform这是基于对象SummarizedExperiment类。转换后的值不再计数，并存储在测定槽中。该colData这是附着dds仍然是可以访问：
<span style="color:#333333"><span style="color:#333333"><code>vsd <- vst(dds, blind = <span style="color:#990073">FALSE</span>)
head(assay(vsd), <span style="color:#009999">3</span>)</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14 10.105781 9.852029 10.169726 9.991545 10.424865
## ENSG00000000419.12 9.692244 9.923647 9.801921 9.798653 9.763455
## ENSG00000000457.13 9.449592 9.312186 9.362754 9.459168 9.281415
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14 10.194490 10.315814 10.002177
## ENSG00000000419.12 9.874703 9.683211 9.845507
## ENSG00000000457.13 9.395937 9.477971 9.477027<span style="color:#333333"><span style="color:#333333"><code>colData(vsd)</code></span></span>## DataFrame with 8 rows and 5 columns
## names donor condition cell dex
## <factor> <factor> <factor> <factor> <factor>
## SRR1039508 SRR1039508 N61311 Untreated N61311 untrt
## SRR1039509 SRR1039509 N61311 Dexamethasone N61311 trt
## SRR1039512 SRR1039512 N052611 Untreated N052611 untrt
## SRR1039513 SRR1039513 N052611 Dexamethasone N052611 trt
## SRR1039516 SRR1039516 N080611 Untreated N080611 untrt
## SRR1039517 SRR1039517 N080611 Dexamethasone N080611 trt
## SRR1039520 SRR1039520 N061011 Untreated N061011 untrt
## SRR1039521 SRR1039521 N061011 Dexamethasone N061011 trt再次，对于rlog：
<span style="color:#333333"><span style="color:#333333"><code>rld <- rlog(dds, blind = <span style="color:#990073">FALSE</span>)
head(assay(rld), <span style="color:#009999">3</span>)</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14 9.482613 9.172197 9.558383 9.346001 9.851349
## ENSG00000000419.12 8.860186 9.150196 9.000042 8.995902 8.951327
## ENSG00000000457.13 8.354790 8.166700 8.236582 8.366693 8.121781
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14 9.587602 9.727248 9.357876
## ENSG00000000419.12 9.091075 8.848782 9.054384
## ENSG00000000457.13 8.282307 8.392384 8.391023在上面的函数调用中，我们指定blind = FALSE，这意味着细胞系和处理之间的差异（设计中的变量）将不会对实验的预期方差平均趋势有所贡献。实验设计不直接用于转换，而仅用于估计计数的总体变化量。对于完全无监督的转换，可以设置 blind = TRUE（默认设置）。
To show the effect of the transformation, in the figure below we plot the first sample against the second, first simply using the log2 function (after adding 1, to avoid taking the log of zero), and then using the VST and rlog-transformed values. For the log2 approach, we need to first estimate size factors to account for sequencing depth, and then specify normalized=TRUE. Sequencing depth correction is done automatically for the vst and rlog.
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"dplyr"</span>)
<span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"ggplot2"</span>)dds <- estimateSizeFactors(dds)df <- bind_rows(as_data_frame(log2(counts(dds, normalized=<span style="color:#990073">TRUE</span>)[, <span style="color:#009999">1</span>:<span style="color:#009999">2</span>]+<span style="color:#009999">1</span>)) %>%mutate(transformation = <span style="color:#dd1144">"log2(x + 1)"</span>),as_data_frame(assay(vsd)[, <span style="color:#009999">1</span>:<span style="color:#009999">2</span>]) %>% mutate(transformation = <span style="color:#dd1144">"vst"</span>),as_data_frame(assay(rld)[, <span style="color:#009999">1</span>:<span style="color:#009999">2</span>]) %>% mutate(transformation = <span style="color:#dd1144">"rlog"</span>))colnames(df)[<span style="color:#009999">1</span>:<span style="color:#009999">2</span>] <- c(<span style="color:#dd1144">"x"</span>, <span style="color:#dd1144">"y"</span>) lvls <- c(<span style="color:#dd1144">"log2(x + 1)"</span>, <span style="color:#dd1144">"vst"</span>, <span style="color:#dd1144">"rlog"</span>)
df$transformation <- factor(df$transformation, levels=lvls)ggplot(df, aes(x = x, y = y)) + geom_hex(bins = <span style="color:#009999">80</span>) +coord_fixed() + facet_grid( . ~ transformation) </code></span></span>
来自两个样本的转换计数散点图。显示的是使用归一化计数的log2转换（左），使用VST（中）和使用rlog（右）的散点图。尽管rlog的大小与log2的大小大致相同，但VST的值越小则向上移动。样本之间的差异（这些散点图中与y = x的偏差）将有助于距离计算和PCA图。
我们可以看到低计数基因（左下角）在普通对数尺度上似乎有太大变化，而VST和rlog压缩了低计数基因的差异，而这些数据所提供的差异表达信息很少。
4.3采样距离
RNA序列分析中有用的第一步通常是评估样品之间的总体相似性：哪些样品彼此相似，哪些不同？这符合实验设计的期望吗？
我们使用R函数dist计算样本之间的欧几里得距离。为了确保我们所有基因的贡献大致相等，我们将其用于VST数据。我们需要使用t转置值矩阵，因为dist函数期望不同的样本为其参数的行，而不同的维度（此处为基因）为列。
<span style="color:#333333"><span style="color:#333333"><code>sampleDists <- dist(t(assay(vsd)))
sampleDists</code></span></span>## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
## SRR1039509 39.42362
## SRR1039512 32.37620 44.93748
## SRR1039513 51.09677 37.18799 41.79886
## SRR1039516 35.59642 47.54671 34.83458 52.05265
## SRR1039517 51.26314 41.58572 46.89609 40.67315 39.74268
## SRR1039520 32.38578 46.96000 30.35980 48.08669 37.07106 50.38349
## SRR1039521 51.49108 37.57383 47.17283 31.45899 52.62276 41.35941
## SRR1039520
## SRR1039509
## SRR1039512
## SRR1039513
## SRR1039516
## SRR1039517
## SRR1039520
## SRR1039521 43.01502我们使用pheatmap包中的pheatmap函数在下图中可视化了热图中的距离 。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"pheatmap"</span>)
<span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"RColorBrewer"</span>)</code></span></span>为了用距离矩阵中的距离排列行/列来绘制样本距离矩阵，我们手动提供 了pheatmap函数sampleDists的clustering_distance参数。否则，pheatmap函数将假定矩阵包含数据值本身，并且将计算距离矩阵的行/列之间的距离，这是不希望的。我们还使用RColorBrewer包中的colorRampPalette函数手动指定了蓝色调色板 。
<span style="color:#333333"><span style="color:#333333"><code>sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- paste( vsd$dex, vsd$cell, sep = <span style="color:#dd1144">" - "</span> )
colnames(sampleDistMatrix) <- <span style="color:#990073">NULL</span>
colors <- colorRampPalette( rev(brewer.pal(<span style="color:#009999">9</span>, <span style="color:#dd1144">"Blues"</span>)) )(<span style="color:#009999">255</span>)
pheatmap(sampleDistMatrix,clustering_distance_rows = sampleDists,clustering_distance_cols = sampleDists,col = colors)</code></span></span>
使用方差稳定化的转换值的样本间距离热图。
请注意，我们已更改了距离矩阵的行名称，以包含治疗类型和患者编号，而不是样本ID，因此在查看热图时，我们可以查看所有这些信息。
计算样本距离的另一种方法是使用PoiClaClu软件包中实现 的 Poisson距离 （Witten 2011）。在计算样本之间的距离时，这种计数差异的度量还考虑了计数的固有方差结构。该PoissonDistance功能将原始计数矩阵（未归一化）的样品为行，而不是列，所以我们需要调换计数。dds
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"PoiClaClu"</span>)
poisd <- PoissonDistance(t(counts(dds)))</code></span></span>我们在下图中绘制热图。
<span style="color:#333333"><span style="color:#333333"><code>samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- paste( dds$dex, dds$cell, sep=<span style="color:#dd1144">" - "</span> )
colnames(samplePoisDistMatrix) <- <span style="color:#990073">NULL</span>
pheatmap(samplePoisDistMatrix,clustering_distance_rows = poisd$dd,clustering_distance_cols = poisd$dd,col = colors)</code></span></span>
使用Poisson距离的样品间距离的热图。
4.4PCA图
可视化样品到样品距离的另一种方法是主成分分析（PCA）。在这种排序方法中，数据点（此处为样本）被投影到2D平面上，以便它们在解释大多数差异的两个方向上展开（下图）。x轴是分隔数据点最多的方向。该方向上的样本值被写入PC1。y轴是第二个分隔数据的方向（必须与第一个方向正交）。该方向上的样本值写入PC2。方向标签中包含的总方差百分比将打印在轴标签中。请注意，这些百分比不会加到100％，因为有更多维度包含剩余方差（尽管这些剩余维度中的每一个将解释的比我们看到的更少）。
<span style="color:#333333"><span style="color:#333333"><code>plotPCA(vsd, intgroup = c(<span style="color:#dd1144">"dex"</span>, <span style="color:#dd1144">"cell"</span>))</code></span></span>
使用VST数据的PCA图。处理和细胞系的每种独特组合都有其自己的颜色。
在这里，我们使用了DESeq2附带的plotPCA函数。指定的两个术语intgroup是用于标记样品的有趣组。他们告诉函数使用它们来选择颜色。我们还可以使用ggplot2软件包从头开始构建PCA图 （Wickham 2009）。这是通过要求plotPCA函数返回用于绘图的数据而不是构建绘图来完成的。见GGPLOT2 文档 有关使用的详细信息ggplot。
<span style="color:#333333"><span style="color:#333333"><code>pcaData <- plotPCA(vsd, intgroup = c( <span style="color:#dd1144">"dex"</span>, <span style="color:#dd1144">"cell"</span>), returnData = <span style="color:#990073">TRUE</span>)
pcaData</code></span></span>## PC1 PC2 group dex cell name
## SRR1039508 -14.311369 -2.6000421 untrt:N61311 untrt N61311 SRR1039508
## SRR1039509 8.058574 -0.7500532 trt:N61311 trt N61311 SRR1039509
## SRR1039512 -9.404122 -4.3920761 untrt:N052611 untrt N052611 SRR1039512
## SRR1039513 14.497842 -4.1323833 trt:N052611 trt N052611 SRR1039513
## SRR1039516 -12.365055 11.2109581 untrt:N080611 untrt N080611 SRR1039516
## SRR1039517 9.343946 14.9115160 trt:N080611 trt N080611 SRR1039517
## SRR1039520 -10.852633 -7.7618618 untrt:N061011 untrt N061011 SRR1039520
## SRR1039521 15.032816 -6.4860576 trt:N061011 trt N061011 SRR1039521<span style="color:#333333"><span style="color:#333333"><code>percentVar <- round(<span style="color:#009999">100</span> * attr(pcaData, <span style="color:#dd1144">"percentVar"</span>))</code></span></span>然后，我们可以使用这些数据在下图中建立第二个图，指定点的颜色应反映地塞米松处理，形状应反映细胞系。
<span style="color:#333333"><span style="color:#333333"><code>ggplot(pcaData, aes(x = PC1, y = PC2, color = dex, shape = cell)) +geom_point(size =<span style="color:#009999">3</span>) +xlab(paste0(<span style="color:#dd1144">"PC1: "</span>, percentVar[<span style="color:#009999">1</span>], <span style="color:#dd1144">"% variance"</span>)) +ylab(paste0(<span style="color:#dd1144">"PC2: "</span>, percentVar[<span style="color:#009999">2</span>], <span style="color:#dd1144">"% variance"</span>)) +coord_fixed() +ggtitle(<span style="color:#dd1144">"PCA with VST data"</span>)</code></span></span>
使用具有自定义ggplot2代码的VST值进行PCA图。 在这里，我们指定细胞系（绘图符号）和地塞米松处理（颜色）。
从PCA图中，我们看到细胞之间的差异（不同的绘制形状）相当大，尽管不强于地塞米松处理（红色与蓝色）所引起的差异。这说明了为什么使用配对设计（“成对”，因为每个经过dex处理的样品与来自同一细胞系的一个未经处理的样品配对）在差异测试中考虑这一点很重要。我们已经通过~ cell + dex较早地分配公式来进行此设计。
4.5使用广义PCA的PCA图
用于对非正态分布的数据（例如，过度分散的计数数据） 执行降维的另一种技术是 通用主成分分析或GLM-PCA （Townes等人，2019年），该方法在CRAN软件包glmpca中实现。该程序包将计数矩阵以及要适应的潜在维数作为输入（此处指定2）。如Townes等人所述。（2019）：
…我们建议使用GLM-PCA，这是PCA对指数族可能性的概括。GLM-PCA对原始计数进行操作，避免了标准化的陷阱。我们还证明了将PCA应用于偏差或Pearson残差可为GLM-PCA提供有用且快速的近似。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"glmpca"</span>)
gpca <- glmpca(counts(dds), L=<span style="color:#009999">2</span>)
gpca.dat <- gpca$factors
gpca.dat$dex <- dds$dex
gpca.dat$cell <- dds$cell</code></span></span><span style="color:#333333"><span style="color:#333333"><code>ggplot(gpca.dat, aes(x = dim1, y = dim2, color = dex, shape = cell)) +geom_point(size =<span style="color:#009999">3</span>) + coord_fixed() + ggtitle(<span style="color:#dd1144">"glmpca - Generalized PCA"</span>)</code></span></span>
4.6MDS图
可以使用底数R中的多维比例缩放（MDS）函数来制作另一个与PCA图非常相似的图 。当我们没有数据矩阵而是只有距离矩阵时，这很有用。在这里，我们针对从VST数据计算出的距离计算MDS，并将其绘制在下图中。
<span style="color:#333333"><span style="color:#333333"><code>mds <- as.data.frame(colData(vsd)) %>%cbind(cmdscale(sampleDistMatrix))
ggplot(mds, aes(x = `1`, y = `2`, color = dex, shape = cell)) +geom_point(size = <span style="color:#009999">3</span>) + coord_fixed() + ggtitle(<span style="color:#dd1144">"MDS with VST data"</span>)</code></span></span>
使用VST数据的MDS图。
在下图中，我们显示了PoissonDistance的相同图：
<span style="color:#333333"><span style="color:#333333"><code>mdsPois <- as.data.frame(colData(dds)) %>%cbind(cmdscale(samplePoisDistMatrix))
ggplot(mdsPois, aes(x = `1`, y = `2`, color = dex, shape = cell)) +geom_point(size = <span style="color:#009999">3</span>) + coord_fixed() + ggtitle(<span style="color:#dd1144">"MDS with PoissonDistances"</span>)</code></span></span>
使用泊松距离的MDS图。
5差异表达分析
5.1运行差异表达式管道
正如我们在创建DESeqDataSet时已经指定实验设计那样，我们可以通过调用DESeq函数对原始计数运行差分表达式管道：
<span style="color:#333333"><span style="color:#333333"><code>dds <- DESeq(dds)</code></span></span>此功能将为执行的各个步骤打印一条消息。这些在DESeq的手册页中有更详细的描述 ，可以通过键入进行访问?DESeq。简要地说，这些是：大小因子的估计（控制样品测序深度的差异），每个基因的分散值估计以及拟合广义线性模型。
返回一个DESeqDataSet，其中包含所有已拟合的参数，以下部分描述了如何从此对象中提取感兴趣的结果表。
5.2建立结果表
不带任何参数的调用结果将提取设计公式中最后一个变量的估计log2倍变化和p值。如果此变量有两个以上级别，则结果 将提取结果表，以比较最后一个级别与第一个级别。比较将显示在输出的顶部： dex trt vs untrt。
<span style="color:#333333"><span style="color:#333333"><code>res <- results(dds)
res</code></span></span>## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 31604 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003.14 739.940717 -0.3611537 0.106869 -3.379419 0.000726392
## ENSG00000000419.12 511.735722 0.2063147 0.128665 1.603509 0.108822318
## ENSG00000000457.13 314.194855 0.0378308 0.158633 0.238479 0.811509461
## ENSG00000000460.16 79.793622 -0.1152590 0.314991 -0.365912 0.714430444
## ENSG00000000938.12 0.307267 -1.3691185 3.503764 -0.390757 0.695977205
## ... ... ... ... ... ...
## ENSG00000285979.1 38.353886 0.3423657 0.359511 0.952310 0.340940
## ENSG00000285987.1 1.562508 0.7064145 1.547295 0.456548 0.647996
## ENSG00000285990.1 0.642315 0.3647333 3.433276 0.106235 0.915396
## ENSG00000285991.1 11.276284 -0.1165515 0.748601 -0.155692 0.876275
## ENSG00000285994.1 3.651041 -0.0960094 1.068660 -0.089841 0.928414
## padj
## <numeric>
## ENSG00000000003.14 0.00531137
## ENSG00000000419.12 0.29318870
## ENSG00000000457.13 0.92255697
## ENSG00000000460.16 0.87298038
## ENSG00000000938.12 NA
## ... ...
## ENSG00000285979.1 0.600750
## ENSG00000285987.1 NA
## ENSG00000285990.1 NA
## ENSG00000285991.1 0.952921
## ENSG00000285994.1 NA我们可以使用以下更具体的命令等效地生成此结果表。因为dex是设计中的最后一个变量，我们可以选择不选择contrast参数来提取的两个级别的比较dex。
<span style="color:#333333"><span style="color:#333333"><code>res <- results(dds, contrast=c(<span style="color:#dd1144">"dex"</span>,<span style="color:#dd1144">"trt"</span>,<span style="color:#dd1144">"untrt"</span>))</code></span></span>由于res是一个数据帧的对象，它带有与所述列的含义信息的元数据：
<span style="color:#333333"><span style="color:#333333"><code>mcols(res, use.names = <span style="color:#990073">TRUE</span>)</code></span></span>## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized c..
## log2FoldChange results log2 fold change (ML..
## lfcSE results standard error: dex ..
## stat results Wald statistic: dex ..
## pvalue results Wald test p-value: d..
## padj results BH adjusted p-values第一列baseMean是DESeqDataSet中所有样本的归一化计数值的平均值除以大小因子 。其余四列涉及特定的对比，即trt，untrt水平与因子变量的水平之间的比较dex。我们将在下面找到如何获得其他对比。
The column log2FoldChange is the effect size estimate. It tells us how much the gene’s expression seems to have changed due to treatment with dexamethasone in comparison to untreated samples. This value is reported on a logarithmic scale to base 2: for example, a log2 fold change of 1.5 means that the gene’s expression is increased by a multiplicative factor of 21.5≈2.8221.5≈2.82.
当然，该估计值具有不确定性，可在lfcSElog2倍变化估计值的标准误差估计值列中找到。我们还可以通过统计检验来表示特定效果大小估计值的不确定性。测试差异表达的目的是测试数据是否提供足够的证据来推断该值确实不同于零。DESeq2对每个基因进行 假设检验，以查看是否有足够的证据来反对无效假设认为治疗对基因的影响为零，并且观察到的治疗和对照之间的差异仅是由实验变异性引起的（即，同一治疗组中不同样本之间可以预期的变异类型）。与统计数据一样，此测试的结果报告为p值，并在列中找到pvalue。请记住，p值表示在原假设所描述的情况下，观察到的倍数变化倍于观察到的倍数甚至更高的概率。
我们还可以使用以下代码行来总结结果，该代码行报告一些其他信息，这些信息将在后面的部分中介绍。
<span style="color:#333333"><span style="color:#333333"><code>summary(res)</code></span></span>##
## out of 31604 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 2373, 7.5%
## LFC < 0 (down) : 1949, 6.2%
## outliers [1] : 0, 0%
## low counts [2] : 14706, 47%
## (mean count < 9)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results注意，由于地塞米松在FDR水平为10％的情况下存在许多差异表达的基因。这是有道理的，因为已知气道的平滑肌细胞会与糖皮质激素类固醇发生反应。但是，有两种方法可以更严格地确定哪些基因被认为是重要的：
- 降低错误发现率阈值（
padj结果表中的阈值） - 使用结果
lfcThreshold参数将log2倍数变化阈值从0提高
如果降低虚假发现率阈值，则还应告知results()函数，以便函数可以将此阈值用于其执行的最佳独立过滤：
<span style="color:#333333"><span style="color:#333333"><code>res.05 <- results(dds, alpha = <span style="color:#009999">0.05</span>)
table(res.05$padj < <span style="color:#009999">0.05</span>)</code></span></span>##
## FALSE TRUE
## 13357 3541如果我们想提高log2倍变化阈值，以便测试因治疗而显示更多实质性变化的基因，我们只需提供log2范围内的值即可。例如，通过指定 lfcThreshold = 1，我们测试显示对基因计数有显着影响的基因多于两倍或少于一半，这是因为21个= 221个=2。
<span style="color:#333333"><span style="color:#333333"><code>resLFC1 <- results(dds, lfcThreshold=<span style="color:#009999">1</span>)
table(resLFC1$padj < <span style="color:#009999">0.1</span>)</code></span></span>##
## FALSE TRUE
## 16687 211有时，其中p值的子集res将是“ NA（不可用）”。这是DESeq的报告该基因的所有计数均为零的方法，因此未进行任何测试。此外，如果基因被排除在分析之外，因为该基因包含一个极端计数异常 值，则可以指定p值NA。有关更多信息，请参见DESeq2小插图的异常值检测部分。
如果在已发表的研究中使用R分析软件包的结果，则可以通过键入来找到该软件的正确引文，用 citation("pkgName")将该软件包的名称替换为pkgName。引用方法的论文有助于支持和奖励那些将时间投入开源软件进行基因组数据分析的人员。
5.3其他比较
通常，可以使用result的contrast参数来 提取变量任意两个级别的比较结果。用户应指定三个值：变量名称，分子级别的名称和分母级别的名称。在这里，我们提取一个细胞系相对于另一个细胞系的倍数变化的log2结果：
<span style="color:#333333"><span style="color:#333333"><code>results(dds, contrast = c(<span style="color:#dd1144">"cell"</span>, <span style="color:#dd1144">"N061011"</span>, <span style="color:#dd1144">"N61311"</span>))</code></span></span>## log2 fold change (MLE): cell N061011 vs N61311
## Wald test p-value: cell N061011 vs N61311
## DataFrame with 31604 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003.14 739.940717 0.270945 0.152171 1.780534 0.0749886
## ENSG00000000419.12 511.735722 -0.071831 0.182817 -0.392912 0.6943842
## ENSG00000000457.13 314.194855 0.179881 0.225122 0.799036 0.4242696
## ENSG00000000460.16 79.793622 -0.119482 0.441594 -0.270570 0.7867217
## ENSG00000000938.12 0.307267 0.000000 4.997580 0.000000 1.0000000
## ... ... ... ... ... ...
## ENSG00000285979.1 38.353886 0.0589757 0.512391 0.1150989 0.908367
## ENSG00000285987.1 1.562508 1.0216804 2.201861 0.4640078 0.642642
## ENSG00000285990.1 0.642315 -3.0956404 4.852715 -0.6379193 0.523526
## ENSG00000285991.1 11.276284 -0.8779628 1.046963 -0.8385804 0.401705
## ENSG00000285994.1 3.651041 -0.0192351 1.513236 -0.0127112 0.989858
## padj
## <numeric>
## ENSG00000000003.14 0.378828
## ENSG00000000419.12 0.936703
## ENSG00000000457.13 0.820733
## ENSG00000000460.16 0.960662
## ENSG00000000938.12 NA
## ... ...
## ENSG00000285979.1 0.98371
## ENSG00000285987.1 NA
## ENSG00000285990.1 NA
## ENSG00000285991.1 NA
## ENSG00000285994.1 NAThere are additional ways to build results tables for certain comparisons after running DESeq once. If results for an interaction term are desired, the name argument of results should be used. Please see the help page for the results function for details on the additional ways to build results tables. In particular, the Examples section of the help page for results gives some pertinent examples.
5.4Multiple testing
In high-throughput biology, we are careful to not use the p values directly as evidence against the null, but to correct for multiple testing. What would happen if we were to simply threshold the p values at a low value, say 0.05? There are 5170 genes with a p value below 0.05 among the 31604 genes for which the test succeeded in reporting a p value:
<span style="color:#333333"><span style="color:#333333"><code>sum(res$pvalue < <span style="color:#009999">0.05</span>, na.rm=<span style="color:#990073">TRUE</span>)</code></span></span>## [1] 5170<span style="color:#333333"><span style="color:#333333"><code>sum(!is.na(res$pvalue))</code></span></span>## [1] 31604Now, assume for a moment that the null hypothesis is true for all genes, i.e., no gene is affected by the treatment with dexamethasone. Then, by the definition of the p value, we expect up to 5% of the genes to have a p value below 0.05. This amounts to 1580 genes. If we just considered the list of genes with a p value below 0.05 as differentially expressed, this list should therefore be expected to contain up to 1580 / 5170 = 31% false positives.
DESeq2 uses the Benjamini-Hochberg (BH) adjustment (Benjamini and Hochberg 1995) as implemented in the base R p.adjust function; in brief, this method calculates for each gene an adjusted p value that answers the following question: if one called significant all genes with an adjusted p value less than or equal to this gene’s adjusted p value threshold, what would be the fraction of false positives (the false discovery rate, FDR) among them, in the sense of the calculation outlined above? These values, called the BH-adjusted p values, are given in the column padj of the res object.
对于许多高通量实验，FDR是有用的统计数据，因为我们经常对报告或关注一组有趣的基因感兴趣，并且我们希望对该组假阳性的百分比设置上限。
因此，如果我们认为一部分10％的假阳性是可以接受的，则可以将所有调整后的p值低于10％= 0.1的基因视为显着。有多少这样的基因？
<span style="color:#333333"><span style="color:#333333"><code>sum(res$padj < <span style="color:#009999">0.1</span>, na.rm=<span style="color:#990073">TRUE</span>)</code></span></span>## [1] 4322我们将结果表子集化为这些基因，然后按log2倍变化估算值对其进行排序，以获得具有最强下调的重要基因：
<span style="color:#333333"><span style="color:#333333"><code>resSig <- subset(res, padj < <span style="color:#009999">0.1</span>)
head(resSig[ order(resSig$log2FoldChange), ])</code></span></span>## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000216490.3 42.3007 -5.72483 1.475652 -3.87952 1.04661e-04
## ENSG00000267339.5 30.5206 -5.39781 0.773017 -6.98278 2.89390e-12
## ENSG00000257542.5 10.0399 -5.25991 1.282001 -4.10289 4.08015e-05
## ENSG00000146006.7 61.6448 -4.49504 0.663821 -6.77147 1.27484e-11
## ENSG00000108700.4 14.6324 -4.09069 0.941842 -4.34328 1.40369e-05
## ENSG00000213240.8 12.0962 -3.87313 1.274133 -3.03981 2.36725e-03
## padj
## <numeric>
## ENSG00000216490.3 9.87853e-04
## ENSG00000267339.5 9.45863e-11
## ENSG00000257542.5 4.30646e-04
## ENSG00000146006.7 3.82632e-10
## ENSG00000108700.4 1.66687e-04
## ENSG00000213240.8 1.44987e-02…以及最强的上调：
<span style="color:#333333"><span style="color:#333333"><code>head(resSig[ order(resSig$log2FoldChange, decreasing = <span style="color:#990073">TRUE</span>), ])</code></span></span>## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000254692.1 62.2302 10.20714 3.37706 3.02250 2.50700e-03
## ENSG00000179593.15 67.0895 9.50515 1.07705 8.82513 1.09334e-18
## ENSG00000268173.3 46.4370 8.40438 3.38506 2.48279 1.30358e-02
## ENSG00000224712.12 35.5607 7.16686 2.16476 3.31070 9.30628e-04
## ENSG00000109906.13 438.1940 6.37750 0.31381 20.32276 8.08934e-92
## ENSG00000257663.1 24.3946 6.34758 2.09531 3.02942 2.45024e-03
## padj
## <numeric>
## ENSG00000254692.1 1.52167e-02
## ENSG00000179593.15 6.81746e-17
## ENSG00000268173.3 5.94385e-02
## ENSG00000224712.12 6.55240e-03
## ENSG00000109906.13 1.24267e-88
## ENSG00000257663.1 1.49151e-026绘图结果
6.1计数图
可视化特定基因计数的一种快速方法是使用plotCounts函数，该函数将DESeqDataSet，基因名称和要绘制计数的组作为参数 （下图）。
<span style="color:#333333"><span style="color:#333333"><code>topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene = topGene, intgroup=c(<span style="color:#dd1144">"dex"</span>))</code></span></span>
治疗组中单个基因的标准化计数。
我们还可以使用ggplot2软件包中的ggplot函数 进行自定义绘图（以下图）。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"ggbeeswarm"</span>)
geneCounts <- plotCounts(dds, gene = topGene, intgroup = c(<span style="color:#dd1144">"dex"</span>,<span style="color:#dd1144">"cell"</span>),returnData = <span style="color:#990073">TRUE</span>)
ggplot(geneCounts, aes(x = dex, y = count, color = cell)) +scale_y_log10() + geom_beeswarm(cex = <span style="color:#009999">3</span>)</code></span></span>
<span style="color:#333333"><span style="color:#333333"><code>ggplot(geneCounts, aes(x = dex, y = count, color = cell, group = cell)) +scale_y_log10() + geom_point(size = <span style="color:#009999">3</span>) + geom_line()</code></span></span>
连接细胞系的线归一化计数。 请注意，DESeq测试实际上考虑了细胞系效应，因此该图更紧密地描绘了所测试的差异。
6.2MA图
An MA-plot (Dudoit et al. 2002) provides a useful overview for the distribution of the estimated coefficients in the model, e.g. the comparisons of interest, across all genes. On the y-axis, the “M” stands for “minus” – subtraction of log values is equivalent to the log of the ratio – and on the x-axis, the “A” stands for “average”. You may hear this plot also referred to as a mean-difference plot, or a Bland-Altman plot.
Before making the MA-plot, we use the lfcShrink function to shrink the log2 fold changes for the comparison of dex treated vs untreated samples. There are three types of shrinkage estimators in DESeq2, which are covered in the DESeq2 vignette. Here we specify the apeglm method for shrinking coefficients, which is good for shrinking the noisy LFC estimates while giving low bias LFC estimates for true large differences (Zhu, Ibrahim, and Love 2018). To use apeglm we specify a coefficient from the model to shrink, either by name or number as the coefficient appears in resultsNames(dds).
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"apeglm"</span>)
resultsNames(dds)</code></span></span>## [1] "Intercept" "cell_N061011_vs_N052611"
## [3] "cell_N080611_vs_N052611" "cell_N61311_vs_N052611"
## [5] "dex_trt_vs_untrt"<span style="color:#333333"><span style="color:#333333"><code>res <- lfcShrink(dds, coef=<span style="color:#dd1144">"dex_trt_vs_untrt"</span>, type=<span style="color:#dd1144">"apeglm"</span>)
plotMA(res, ylim = c(-<span style="color:#009999">5</span>, <span style="color:#009999">5</span>))</code></span></span>
如果有必要，指定在未示出的对比度 resultsNames(dds)任一的其它两个收缩的方法，可以使用，或者在一些情况下，再保相关的变量和运行 nbinomWaldTest，接着lfcShrink是足够的。有关更多详细信息，请参见DESeq2小插图。
由治疗引起的变化的MA图。 在y轴上绘制了特定比较的log2倍数变化，在x轴上显示了通过尺寸系数归一化的计数平均值。每个基因都用一个点表示。调整后的p值低于阈值（此处为默认值0.1）的基因以红色显示。
所述DESeq2包使用贝叶斯方法，中度（或“收缩”）log 2倍数具有非常低的计数和高度可变的计数从基因的变化，这可以通过点的垂直传播在MA-左侧变窄可以看出情节。如上所示， lfcShrink函数执行此操作。有关适度倍数变化原理的详细说明，请参见 DESeq2论文（Love，Huber和Anders 2014）。
如果我们没有使用统计调节来缩小嘈杂的log2倍变化，那么我们将看到以下图：
<span style="color:#333333"><span style="color:#333333"><code>res.noshr <- results(dds, name=<span style="color:#dd1144">"dex_trt_vs_untrt"</span>)
plotMA(res.noshr, ylim = c(-<span style="color:#009999">5</span>, <span style="color:#009999">5</span>))</code></span></span>
我们也可以在MA图上标记单个点。在这里，我们使用 with R函数为结果对象的选定行绘制圆和文本。在with函数中，仅使用baseMean和 log2FoldChange的选定行的值res。
<span style="color:#333333"><span style="color:#333333"><code>plotMA(res, ylim = c(-<span style="color:#009999">5</span>,<span style="color:#009999">5</span>))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {points(baseMean, log2FoldChange, col=<span style="color:#dd1144">"dodgerblue"</span>, cex=<span style="color:#009999">2</span>, lwd=<span style="color:#009999">2</span>)text(baseMean, log2FoldChange, topGene, pos=<span style="color:#009999">2</span>, col=<span style="color:#dd1144">"dodgerblue"</span>)
})</code></span></span>
另一个有用的诊断图是p值的直方图（下图）。该图最好通过排除计数非常少的基因来形成，否则计数会在直方图中产生尖峰。
<span style="color:#333333"><span style="color:#333333"><code>hist(res$pvalue[res$baseMean > <span style="color:#009999">1</span>], breaks = <span style="color:#009999">0</span>:<span style="color:#009999">20</span>/<span style="color:#009999">20</span>,col = <span style="color:#dd1144">"grey50"</span>, border = <span style="color:#dd1144">"white"</span>)</code></span></span>
平均归一化计数大于1的基因的p值直方图。
6.3基因聚类
在先前制作的样本距离热图中，侧面的树状图显示了样本的层次聚类。也可以对基因进行这种聚类。由于聚类仅与实际携带信号的基因有关，因此通常只聚类高度可变基因的子集。在这里，为进行演示，让我们选择样本中方差最大的20个基因。我们将处理VST数据。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"genefilter"</span>)
topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = <span style="color:#990073">TRUE</span>), <span style="color:#009999">20</span>)</code></span></span>如果我们不看绝对表达强度，而是看每个基因在特定样品中偏离所有样品的基因平均值的数量，那么热图就会变得更加有趣。因此，我们将每个基因的值在样本中居中，并绘制一个热图（下图）。我们提供了一个data.frame，它指示 Pheatmap函数如何标记列。
<span style="color:#333333"><span style="color:#333333"><code>mat <- assay(vsd)[ topVarGenes, ]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c(<span style="color:#dd1144">"cell"</span>,<span style="color:#dd1144">"dex"</span>)])
pheatmap(mat, annotation_col = anno)</code></span></span>
样本之间相对VST转换值的热图。 在热图顶部以彩色条显示治疗状态和细胞系信息。跨患者变异的基因块。请注意，热图中的一组基因将N061011细胞株与其他细胞分开，并且还有另一组基因，地塞米松处理的样品具有更高的基因表达。
6.4独立过滤
MA图突出显示了RNA序列数据的重要属性。对于弱表达的基因，我们没有机会看到差异表达，因为低读数计数会遭受如此高的泊松噪声，以至于任何生物效应都会被低速率采样所带来的不确定性淹没。我们还可以通过检查平均归一化计数归类的基因的小p值（例如，小于0.05）的比率来证明这一点。我们将使用经过阈值的结果表来显示在很少有带有小p 值的测试的情况下的结果。
在以下代码块中，我们使用分位数 功能创建bin，使用cut通过基本均值对基因进行bin ，使用中点重命名bin的级别，计算每个bin的p值小于0.05的比率，最后进行绘图这些比率（下图）。
<span style="color:#333333"><span style="color:#333333"><code>qs <- c(<span style="color:#009999">0</span>, quantile(resLFC1$baseMean[resLFC1$baseMean > <span style="color:#009999">0</span>], <span style="color:#009999">0</span>:<span style="color:#009999">6</span>/<span style="color:#009999">6</span>))
bins <- cut(resLFC1$baseMean, qs)
levels(bins) <- paste0(<span style="color:#dd1144">"~"</span>, round(signif((qs[-<span style="color:#009999">1</span>] + qs[-length(qs)])/<span style="color:#009999">2</span>, <span style="color:#009999">2</span>)))
fractionSig <- tapply(resLFC1$pvalue, bins, <span style="color:#990000"><strong>function</strong></span>(p)mean(p < <span style="color:#009999">.05</span>, na.rm = <span style="color:#990073">TRUE</span>))
barplot(fractionSig, xlab = <span style="color:#dd1144">"mean normalized count"</span>,ylab = <span style="color:#dd1144">"fraction of small p values"</span>)</code></span></span>
基因的小p值与平均归一化计数合并的比率。 的p值是来自log 2倍数变化更大的大于1或小于-1的测试。该图表明，具有极低平均计数的基因几乎没有功效，甚至没有功效，因此最好从测试中排除。
乍一看，过滤掉这些基因似乎没有什么好处。毕竟，测试发现它们无论如何都不重要。但是，这些基因会影响多重测试调整，如果删除此类基因，其性能会提高。通过从FDR程序的输入中删除低计数基因，我们可以发现我们保留的基因中有更多重要基因，从而提高了测试的功效。这种方法称为独立过滤。
该DESeq2软件自动进行最大化基因的调节的数目独立过滤p比（默认情况下，临界值值以下alpha设定为0.1）。此自动独立过滤由结果功能执行，并可以由结果功能控制。
独立一词突出了一个重要的警告。仅当我们要过滤的统计量（此处所有样本的归一化计数的平均值）与零假设下的实际检验统计量（p值）无关时，才允许进行此类过滤。否则，过滤将使测试无效，从而使BH程序的假设无效。DESeq2内部使用的独立过滤软件 来自genefilter软件包，其中包含对描述独立过滤统计基础的论文的引用（Bourgon，Gentleman和Huber 2010）。
6.5独立假设权重
A generalization of the idea of p value filtering is to weight hypotheses to optimize power. A Bioconductor package, IHW is available that implements the method of Independent Hypothesis Weighting (Ignatiadis et al. 2016). See the DESeq2 package vignette for an example of using IHW in combination with DESeq2. In particular, the following (here, un-evaluated) code chunk can be used to perform IHW in lieu of independent filtering described above.
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"IHW"</span>)
res.ihw <- results(dds, filterFun=ihw)</code></span></span>7Annotating and exporting results
到目前为止，我们的结果表仅包含Ensembl基因ID，但其他基因名称可能对解释更有用。Bioconductor的注释包有助于将各种ID方案相互映射。我们加载AnnotationDbi包和注释包 org.Hs.eg.db：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"AnnotationDbi"</span>)
<span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"org.Hs.eg.db"</span>)</code></span></span>这是智人（“ Hs”）的有机体注释包（“ org”） ，使用Entrez基因ID（“ eg”）作为主键，组织为AnnotationDbi数据库包（“ db”）。要获取所有可用键类型的列表，请使用：
<span style="color:#333333"><span style="color:#333333"><code>columns(org.Hs.eg.db)</code></span></span>## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GO" "GOALL" "IPI" "MAP" "OMIM"
## [16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" "UNIGENE"
## [26] "UNIPROT"我们可以使用mapIds函数将单独的列添加到结果表中。我们提供结果表的行名称作为键，并指定keytype=ENSEMBL。该column参数告诉 mapIds我们想要的信息功能和multiVals 参数告诉函数，如果有一个单一的输入值的多个可能值做什么。在这里，我们要求将我们在数据库中出现的第一个信息还给我们。要添加基因符号和Entrez ID，我们两次调用mapIds。
<span style="color:#333333"><span style="color:#333333"><code>ens.str <- substr(rownames(res), <span style="color:#009999">1</span>, <span style="color:#009999">15</span>)
res$symbol <- mapIds(org.Hs.eg.db,keys=ens.str,column=<span style="color:#dd1144">"SYMBOL"</span>,keytype=<span style="color:#dd1144">"ENSEMBL"</span>,multiVals=<span style="color:#dd1144">"first"</span>)
res$entrez <- mapIds(org.Hs.eg.db,keys=ens.str,column=<span style="color:#dd1144">"ENTREZID"</span>,keytype=<span style="color:#dd1144">"ENSEMBL"</span>,multiVals=<span style="color:#dd1144">"first"</span>)</code></span></span>现在，结果具有所需的外部基因ID：
<span style="color:#333333"><span style="color:#333333"><code>resOrdered <- res[order(res$pvalue),]
head(resOrdered)</code></span></span>## log2 fold change (MAP): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 7 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000189221.9 2373.805 3.38765 0.136985 1.84494e-137 3.11758e-133
## ENSG00000120129.5 3420.727 2.96335 0.120850 6.35042e-135 5.36547e-131
## ENSG00000101347.9 14125.584 3.74129 0.157927 2.88983e-127 1.62775e-123
## ENSG00000196136.17 2710.217 3.23518 0.143951 3.68488e-114 1.55668e-110
## ENSG00000152583.12 974.737 4.48641 0.201276 2.94551e-113 9.95466e-110
## ENSG00000211445.11 12512.792 3.75875 0.169536 2.36246e-112 6.65348e-109
## symbol entrez
## <character> <character>
## ENSG00000189221.9 MAOA 4128
## ENSG00000120129.5 DUSP1 1843
## ENSG00000101347.9 SAMHD1 25939
## ENSG00000196136.17 SERPINA3 12
## ENSG00000152583.12 SPARCL1 8404
## ENSG00000211445.11 GPX3 28787.1汇出结果
您可以轻松地将结果表保存在CSV文件中，然后可以与电子表格程序（例如Excel）共享或加载。到呼叫 as.data.frame需要将转换数据帧对象（IRanges封装）至data.frame可以由加工对象物write.csv。在这里，我们仅采用前100个基因进行演示。
<span style="color:#333333"><span style="color:#333333"><code>resOrderedDF <- as.data.frame(resOrdered)[<span style="color:#009999">1</span>:<span style="color:#009999">100</span>, ]
write.csv(resOrderedDF, file = <span style="color:#dd1144">"results.csv"</span>)</code></span></span>Bioconductor包ReportingTools （Huntley et al。2013）的输出结果的更复杂方法。 ReportingTools将自动生成动态HTML文档，包括使用基因标识符和框图汇总到各个组的标准化计数的外部数据库链接。有关完整的详细信息，请参见ReportingTools插页。使用以下代码执行创建动态ReportingTools报表的最简单版本：
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"ReportingTools"</span>)
htmlRep <- HTMLReport(shortName=<span style="color:#dd1144">"report"</span>, title=<span style="color:#dd1144">"My report"</span>,reportDirectory=<span style="color:#dd1144">"./report"</span>)
publish(resOrderedDF, htmlRep)
url <- finish(htmlRep)
browseURL(url)</code></span></span>7.2绘制基因组空间的倍数变化
如果我们已使用tximeta函数读取定量数据，则我们的DESeqDataSet对象将建立在指定基因的基因组坐标的即用型Bioconductor对象之上。因此，我们可以轻松地在基因组空间中绘制差异表达结果。虽然结果或lfcShrink函数默认情况下使用参数返回DataFrame，但format我们可以要求 GRanges或GRangesList输出（仅当我们在创建DESeqDataSet的上游使用tximeta包中的addExons函数时，后者才是可能的）。
<span style="color:#333333"><span style="color:#333333"><code>resGR <- lfcShrink(dds, coef=<span style="color:#dd1144">"dex_trt_vs_untrt"</span>, type=<span style="color:#dd1144">"apeglm"</span>, format=<span style="color:#dd1144">"GRanges"</span>)
resGR</code></span></span>## GRanges object with 31604 ranges and 5 metadata columns:
## seqnames ranges strand | baseMean
## <Rle> <IRanges> <Rle> | <numeric>
## ENSG00000000003.14 chrX 100627109-100639991 - | 739.940717
## ENSG00000000419.12 chr20 50934867-50958555 - | 511.735722
## ENSG00000000457.13 chr1 169849631-169894267 - | 314.194855
## ENSG00000000460.16 chr1 169662007-169854080 + | 79.793622
## ENSG00000000938.12 chr1 27612064-27635277 - | 0.307267
## ... ... ... ... . ...
## ENSG00000285979.1 chr16 57177349-57181390 + | 38.353886
## ENSG00000285987.1 chr9 84316514-84657077 + | 1.562508
## ENSG00000285990.1 chr14 19244904-19269380 - | 0.642315
## ENSG00000285991.1 chr6 149817937-149896011 - | 11.276284
## ENSG00000285994.1 chr10 12563151-12567351 + | 3.651041
## log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003.14 -0.3360046 0.105909 0.000726392 0.00531137
## ENSG00000000419.12 0.1783789 0.122440 0.108822318 0.29318870
## ENSG00000000457.13 0.0299361 0.141095 0.811509461 0.92255697
## ENSG00000000460.16 -0.0555061 0.222787 0.714430444 0.87298038
## ENSG00000000938.12 -0.0115799 0.304740 0.695977205 NA
## ... ... ... ... ...
## ENSG00000285979.1 0.15284386 0.257070 0.340940 0.600750
## ENSG00000285987.1 0.02551527 0.300687 0.647996 NA
## ENSG00000285990.1 -0.00018563 0.304465 0.915396 NA
## ENSG00000285991.1 -0.01507882 0.283931 0.876275 0.952921
## ENSG00000285994.1 -0.00684681 0.293399 0.928414 NA
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genome我们需要再次添加符号以标记图上的基因：
<span style="color:#333333"><span style="color:#333333"><code>ens.str <- substr(names(resGR), <span style="color:#009999">1</span>, <span style="color:#009999">15</span>)
resGR$symbol <- mapIds(org.Hs.eg.db, ens.str, <span style="color:#dd1144">"SYMBOL"</span>, <span style="color:#dd1144">"ENSEMBL"</span>)</code></span></span>我们将使用Gviz软件包绘制GRanges和相关的元数据：由于地塞米松处理，对数倍数变化。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"Gviz"</span>)</code></span></span>以下代码块指定了p值最小的基因上游和下游100万个碱基对的窗口。我们为窗口内的基因创建了全部结果的子集。如果符号存在并且在子集中没有重复，则将基因符号添加为名称。
<span style="color:#333333"><span style="color:#333333"><code>window <- resGR[topGene] + <span style="color:#009999">1e6</span>
strand(window) <- <span style="color:#dd1144">"*"</span>
resGRsub <- resGR[resGR %over% window]
naOrDup <- is.na(resGRsub$symbol) | duplicated(resGRsub$symbol)
resGRsub$group <- ifelse(naOrDup, names(resGRsub), resGRsub$symbol)</code></span></span>我们创建一个向量，指定此子集中的基因是否具有的低值padj。
<span style="color:#333333"><span style="color:#333333"><code>status <- factor(ifelse(resGRsub$padj < <span style="color:#009999">0.05</span> & !is.na(resGRsub$padj),<span style="color:#dd1144">"sig"</span>, <span style="color:#dd1144">"notsig"</span>))</code></span></span>然后，我们可以使用Gviz函数绘制结果（如下图）。我们创建一个轴轨迹来指定我们在基因组中的位置，一个轨迹将显示基因及其名称（按重要性进行着色），一个数据轨迹将绘制垂直条以显示由DESeq2产生的温和对数倍数变化，我们知道这是仅当计数中的信息很好地支持效果时才大。
<span style="color:#333333"><span style="color:#333333"><code>options(ucscChromosomeNames = <span style="color:#990073">FALSE</span>)
g <- GenomeAxisTrack()
a <- AnnotationTrack(resGRsub, name = <span style="color:#dd1144">"gene ranges"</span>, feature = status)
d <- DataTrack(resGRsub, data = <span style="color:#dd1144">"log2FoldChange"</span>, baseline = <span style="color:#009999">0</span>,type = <span style="color:#dd1144">"h"</span>, name = <span style="color:#dd1144">"log2 fold change"</span>, strand = <span style="color:#dd1144">"+"</span>)
plotTracks(list(g, d, a), groupAnnotation = <span style="color:#dd1144">"group"</span>,notsig = <span style="color:#dd1144">"grey"</span>, sig = <span style="color:#dd1144">"hotpink"</span>)</code></span></span>
对数最小的p值在基因周围的基因组区域中发生log2倍变化。粉红色突出显示的基因调整后的p 值小于0.1。
8删除隐藏的批处理效果
假设我们不知道实验涉及不同的细胞系，只是知道地塞米松可以治疗。然后，细胞系对计数的影响将代表一些隐藏的和有害的变异，这些变异可能会影响数据集中的许多或所有基因。我们可以使用Bioconductor中sva软件包 （Leek 2014）或RUVSeq软件包 （Risso等人2014）中针对RNA-seq设计的统计方法来检测样本的此类分组，然后将它们添加到 DESeqDataSet设计中，为了解决它们。
该SVA包使用术语替代变量的估计变量，我们要占在我们的分析，而RUV包使用条款不需要变化的因素，与首字母缩写“删除不需要的变化”解释包称号。我们首先使用SVA查找隐藏的批处理效果，然后使用RUV进行查找。
8.1将SVA与DESeq2一起使用
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"sva"</span>)</code></span></span>在下面，我们获得归一化计数的矩阵，其样本间的平均计数大于1。如上所述，假设我们不知道细胞系信息，我们正在尝试恢复任何隐藏的批次效应。因此，我们使用具有dex变量的完整模型矩阵 ，以及仅具有拦截项的简化模型或空模型矩阵。最后，我们指定要估计2个代理变量。有关更多信息，请键入，以了解svaseq函数的手册页 ?svaseq。
<span style="color:#333333"><span style="color:#333333"><code>dat <- counts(dds, normalized = <span style="color:#990073">TRUE</span>)
idx <- rowMeans(dat) > <span style="color:#009999">1</span>
dat <- dat[idx, ]
mod <- model.matrix(~ dex, colData(dds))
mod0 <- model.matrix(~ <span style="color:#009999">1</span>, colData(dds))
svseq <- svaseq(dat, mod, mod0, n.sv = <span style="color:#009999">2</span>)</code></span></span>## Number of significant surrogate variables is: 2
## Iteration (out of 5 ):1 2 3 4 5<span style="color:#333333"><span style="color:#333333"><code>svseq$sv</code></span></span>## [,1] [,2]
## [1,] 0.2465669 -0.51599084
## [2,] 0.2588137 -0.59462876
## [3,] 0.1384516 0.24920662
## [4,] 0.2179075 0.37716083
## [5,] -0.6042910 -0.06305844
## [6,] -0.6138795 -0.03623320
## [7,] 0.1821306 0.30328185
## [8,] 0.1743002 0.28026195因为我们确实知道细胞系，所以我们可以看到SVA方法在恢复这些变量方面做得如何（如下图）。
<span style="color:#333333"><span style="color:#333333"><code>par(mfrow = c(<span style="color:#009999">2</span>, <span style="color:#009999">1</span>), mar = c(<span style="color:#009999">3</span>,<span style="color:#009999">5</span>,<span style="color:#009999">3</span>,<span style="color:#009999">1</span>))
<span style="color:#990000"><strong>for</strong></span> (i <span style="color:#990000"><strong>in</strong></span> <span style="color:#009999">1</span>:<span style="color:#009999">2</span>) {stripchart(svseq$sv[, i] ~ dds$cell, vertical = <span style="color:#990073">TRUE</span>, main = paste0(<span style="color:#dd1144">"SV"</span>, i))abline(h = <span style="color:#009999">0</span>)}</code></span></span>
在细胞系上绘制变量1和2的替代值。 在这里，我们知道了隐藏的变异源（细胞系），因此可以看到SVA程序如何识别与细胞系相关的变异源。
最后，为了使用SVA消除替代变量对计数的任何影响，我们只需将这两个替代变量作为列添加到DESeqDataSet中，然后将它们添加到设计中：
<span style="color:#333333"><span style="color:#333333"><code>ddssva <- dds
ddssva$SV1 <- svseq$sv[,<span style="color:#009999">1</span>]
ddssva$SV2 <- svseq$sv[,<span style="color:#009999">2</span>]
design(ddssva) <- ~ SV1 + SV2 + dex</code></span></span>然后，通过使用新设计运行DESeq，我们可以生成控制替代变量的结果。
8.2将RUV与DESeq2一起使用
我们还可以使用RUVSeq包中的RUV方法来检测隐藏的批处理效果。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"RUVSeq"</span>)</code></span></span>与SVA的替代变量类似，我们可以使用该RUVg函数来估算不必要的变化因素。与上面的SVA过程相比，不同之处在于，我们首先将运行 DESeq和结果以获取分析的p值，而无需了解批次，例如。假设我们有此结果表，然后通过查看不具有小的p值的基因，得出一组经验控制基因。~ dexres
<span style="color:#333333"><span style="color:#333333"><code>set <- newSeqExpressionSet(counts(dds))
idx <- rowSums(counts(set) > <span style="color:#009999">5</span>) >= <span style="color:#009999">2</span>
set <- set[idx, ]
set <- betweenLaneNormalization(set, which=<span style="color:#dd1144">"upper"</span>)
not.sig <- rownames(res)[which(res$pvalue > <span style="color:#009999">.1</span>)]
empirical <- rownames(set)[ rownames(set) %<span style="color:#990000"><strong>in</strong></span>% not.sig ]
set <- RUVg(set, empirical, k=<span style="color:#009999">2</span>)
pData(set)</code></span></span>## W_1 W_2
## SRR1039508 -0.224881168 0.42992983
## SRR1039509 -0.249022928 0.53858506
## SRR1039512 0.001460949 0.01437385
## SRR1039513 -0.175547525 0.08408354
## SRR1039516 0.599387535 -0.02512358
## SRR1039517 0.590516825 -0.02549392
## SRR1039520 -0.241071948 -0.50369551
## SRR1039521 -0.300841739 -0.51265927我们可以绘制RUV估计的因素：
<span style="color:#333333"><span style="color:#333333"><code>par(mfrow = c(<span style="color:#009999">2</span>, <span style="color:#009999">1</span>), mar = c(<span style="color:#009999">3</span>,<span style="color:#009999">5</span>,<span style="color:#009999">3</span>,<span style="color:#009999">1</span>))
<span style="color:#990000"><strong>for</strong></span> (i <span style="color:#990000"><strong>in</strong></span> <span style="color:#009999">1</span>:<span style="color:#009999">2</span>) {stripchart(pData(set)[, i] ~ dds$cell, vertical = <span style="color:#990073">TRUE</span>, main = paste0(<span style="color:#dd1144">"W"</span>, i))abline(h = <span style="color:#009999">0</span>)}</code></span></span>
在细胞系上绘制了不想要的变化的因素。
和以前一样，如果我们想控制这些因素，只需将它们添加到DESeqDataSet和设计中：
<span style="color:#333333"><span style="color:#333333"><code>ddsruv <- dds
ddsruv$W1 <- set$W_1
ddsruv$W2 <- set$W_2
design(ddsruv) <- ~ W1 + W2 + dex</code></span></span>然后，我们将使用新设计运行DESeq，以重新估计参数和结果。
9时程实验
DESeq2可用于分析时程实验，例如，找到与一组基线样本相比随时间以特定于条件的方式发生反应的那些基因。在这里，我们展示了使用裂变数据包进行的基本时程分析，该 数据包包含裂变酵母RNA序列时程的基因计数（Leong等人，2014）。酵母暴露于氧化应激，一半的样品含有atf21基因的缺失。我们使用一个设计公式来模拟时间0处的应变差异，随时间变化的差异以及随时间变化的任何应变特定的差异（交互作用项strain:minute）。
<span style="color:#333333"><span style="color:#333333"><code><span style="color:#990000"><strong>library</strong></span>(<span style="color:#dd1144">"fission"</span>)
data(<span style="color:#dd1144">"fission"</span>)
ddsTC <- DESeqDataSet(fission, ~ strain + minute + strain:minute)</code></span></span>下面的代码块执行似然比测试，其中我们去除了随时间变化的应变特定差异。此测试中p值小的基因是那些在时间0后一个或多个时间点显示出菌株特异性作用的基因。因此请注意，这不会给两个菌株中以相同方式随时间上下移动的基因提供小的p值。
<span style="color:#333333"><span style="color:#333333"><code>ddsTC <- DESeq(ddsTC, test=<span style="color:#dd1144">"LRT"</span>, reduced = ~ strain + minute)
resTC <- results(ddsTC)
resTC$symbol <- mcols(ddsTC)$symbol
head(resTC[order(resTC$padj),], <span style="color:#009999">4</span>)</code></span></span>## log2 fold change (MLE): strainmut.minute180
## LRT p-value: '~ strain + minute + strain:minute' vs '~ strain + minute'
## DataFrame with 4 rows and 7 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## SPBC2F12.09c 174.671 -2.6567195 0.752261 97.2834 1.97415e-19
## SPAC1002.18 444.505 -0.0509321 0.204299 56.9536 5.16955e-11
## SPAC1002.19 336.373 -0.3927490 0.573494 43.5339 2.87980e-08
## SPAC1002.17c 261.773 -1.1387648 0.606129 39.3158 2.05137e-07
## padj symbol
## <numeric> <character>
## SPBC2F12.09c 1.33453e-15 atf21
## SPAC1002.18 1.74731e-07 urg3
## SPAC1002.19 6.48916e-05 urg1
## SPAC1002.17c 3.46682e-04 urg2这只是可以应用于时间序列数据的测试之一。另一个选择是将计数建模为时间的平滑函数，并包括条件与平滑函数的交互项。可以使用R中的样条基函数建立这样的模型，另一种更现代的方法是使用高斯过程（Tonner et al.2017）。
我们可以使用ggplot2绘制随时间变化的组计数 ，用于调整后的p值最小的基因，测试条件相关的时间曲线并考虑时间0的差异（下图）。请记住，在互动方面是 差在占区别在时间0之后给定的时间，两组之间。
<span style="color:#333333"><span style="color:#333333"><code>fiss <- plotCounts(ddsTC, which.min(resTC$padj), intgroup = c(<span style="color:#dd1144">"minute"</span>,<span style="color:#dd1144">"strain"</span>), returnData = <span style="color:#990073">TRUE</span>)
fiss$minute <- as.numeric(as.character(fiss$minute))
ggplot(fiss,aes(x = minute, y = count, color = strain, group = strain)) + geom_point() + stat_summary(fun.y=mean, geom=<span style="color:#dd1144">"line"</span>) +scale_y_log10()</code></span></span>随时间变化具有特定条件条件的基因的归一化计数。
可以使用结果test参数来调查在各个时间点的log2倍数变化的Wald测试：
<span style="color:#333333"><span style="color:#333333"><code>resultsNames(ddsTC)</code></span></span>## [1] "Intercept" "strain_mut_vs_wt" "minute_15_vs_0"
## [4] "minute_30_vs_0" "minute_60_vs_0" "minute_120_vs_0"
## [7] "minute_180_vs_0" "strainmut.minute15" "strainmut.minute30"
## [10] "strainmut.minute60" "strainmut.minute120" "strainmut.minute180"<span style="color:#333333"><span style="color:#333333"><code>res30 <- results(ddsTC, name=<span style="color:#dd1144">"strainmut.minute30"</span>, test=<span style="color:#dd1144">"Wald"</span>)
res30[which.min(resTC$padj),]</code></span></span>## log2 fold change (MLE): strainmut.minute30
## Wald test p-value: strainmut.minute30
## DataFrame with 1 row and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## SPBC2F12.09c 174.671 -2.60047 0.634343 -4.09947 4.14099e-05 0.279931我们还可以通过它们的概况对重要基因进行聚类。我们使用coef 函数提取log2倍数变化的矩阵。请注意，这些是最大似然估计（MLE）。对于缩小的LFC，必须使用来一次获得一个系数lfcShrink。
<span style="color:#333333"><span style="color:#333333"><code>betas <- coef(ddsTC)
colnames(betas)</code></span></span>## [1] "Intercept" "strain_mut_vs_wt" "minute_15_vs_0"
## [4] "minute_30_vs_0" "minute_60_vs_0" "minute_120_vs_0"
## [7] "minute_180_vs_0" "strainmut.minute15" "strainmut.minute30"
## [10] "strainmut.minute60" "strainmut.minute120" "strainmut.minute180"现在，我们可以在热图中绘制log2倍数变化（如下图）。
<span style="color:#333333"><span style="color:#333333"><code>topGenes <- head(order(resTC$padj),<span style="color:#009999">20</span>)
mat <- betas[topGenes, -c(<span style="color:#009999">1</span>,<span style="color:#009999">2</span>)]
thr <- <span style="color:#009999">3</span>
mat[mat < -thr] <- -thr
mat[mat > thr] <- thr
pheatmap(mat, breaks=seq(from=-thr, to=thr, length=<span style="color:#009999">101</span>),cluster_col=<span style="color:#990073">FALSE</span>)</code></span></span>
对于p值最小的调整基因，log2倍变化的热图。 底部的基因组在15-60分钟内对基线样品表现出强烈的表达诱导作用（左下角的红色框），但是突变株的差异却很小（在右下角的框中显示）。
10会议信息
作为本文档的最后一部分，我们调用函数sessionInfo，该函数报告R的版本号以及此会话中使用的所有软件包。良好的做法是始终记录此类记录，因为如果在一个软件包的较新版本中更改了功能，R脚本停止工作或给出不同的结果，它将有助于跟踪发生的情况。通过将其包括在脚本的底部，您的报告将变得更加可重复。
会议信息还应始终 包含在发送给Bioconductor支持网站的任何电子邮件中， 以及分析中使用的所有代码。
<span style="color:#333333"><span style="color:#333333"><code>sessionInfo()</code></span></span>## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] fission_1.10.0 RUVSeq_1.24.0
## [3] edgeR_3.32.0 limma_3.46.0
## [5] EDASeq_2.24.0 ShortRead_1.48.0
## [7] GenomicAlignments_1.26.0 Rsamtools_2.6.0
## [9] Biostrings_2.58.0 XVector_0.30.0
## [11] sva_3.38.0 BiocParallel_1.24.0
## [13] mgcv_1.8-33 nlme_3.1-150
## [15] Gviz_1.34.0 org.Hs.eg.db_3.12.0
## [17] genefilter_1.72.0 apeglm_1.12.0
## [19] ggbeeswarm_0.6.0 glmpca_0.2.0
## [21] PoiClaClu_1.0.2.1 RColorBrewer_1.1-2
## [23] pheatmap_1.0.12 ggplot2_3.3.2
## [25] dplyr_1.0.2 vsn_3.58.0
## [27] DESeq2_1.30.0 magrittr_1.5
## [29] GenomicFeatures_1.42.0 AnnotationDbi_1.52.0
## [31] tximeta_1.8.0 airway_1.10.0
## [33] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [35] GenomicRanges_1.42.0 GenomeInfoDb_1.26.0
## [37] IRanges_2.24.0 S4Vectors_0.28.0
## [39] BiocGenerics_0.36.0 MatrixGenerics_1.2.0
## [41] matrixStats_0.57.0 rmarkdown_2.5
## [43] knitr_1.30 BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.10 Hmisc_4.4-1
## [3] AnnotationHub_2.22.0 aroma.light_3.20.0
## [5] BiocFileCache_1.14.0 plyr_1.8.6
## [7] lazyeval_0.2.2 splines_4.0.3
## [9] digest_0.6.27 ensembldb_2.14.0
## [11] htmltools_0.5.0 magick_2.5.0
## [13] checkmate_2.0.0 memoise_1.1.0
## [15] BSgenome_1.58.0 cluster_2.1.0
## [17] readr_1.4.0 annotate_1.68.0
## [19] R.utils_2.10.1 askpass_1.1
## [21] bdsmatrix_1.3-4 prettyunits_1.1.1
## [23] jpeg_0.1-8.1 colorspace_1.4-1
## [25] blob_1.2.1 rappdirs_0.3.1
## [27] xfun_0.18 crayon_1.3.4
## [29] RCurl_1.98-1.2 jsonlite_1.7.1
## [31] tximport_1.18.0 hexbin_1.28.1
## [33] VariantAnnotation_1.36.0 survival_3.2-7
## [35] glue_1.4.2 gtable_0.3.0
## [37] zlibbioc_1.36.0 DelayedArray_0.16.0
## [39] scales_1.1.1 mvtnorm_1.1-1
## [41] DBI_1.1.0 Rcpp_1.0.5
## [43] htmlTable_2.1.0 xtable_1.8-4
## [45] progress_1.2.2 emdbook_1.3.12
## [47] foreign_0.8-80 bit_4.0.4
## [49] preprocessCore_1.52.0 Formula_1.2-4
## [51] htmlwidgets_1.5.2 httr_1.4.2
## [53] ellipsis_0.3.1 R.methodsS3_1.8.1
## [55] pkgconfig_2.0.3 XML_3.99-0.5
## [57] farver_2.0.3 nnet_7.3-14
## [59] dbplyr_1.4.4 locfit_1.5-9.4
## [61] tidyselect_1.1.0 labeling_0.4.2
## [63] rlang_0.4.8 later_1.1.0.1
## [65] munsell_0.5.0 BiocVersion_3.12.0
## [67] tools_4.0.3 generics_0.0.2
## [69] RSQLite_2.2.1 evaluate_0.14
## [71] stringr_1.4.0 fastmap_1.0.1
## [73] yaml_2.2.1 bit64_4.0.5
## [75] purrr_0.3.4 AnnotationFilter_1.14.0
## [77] mime_0.9 R.oo_1.24.0
## [79] xml2_1.3.2 biomaRt_2.46.0
## [81] rstudioapi_0.11 compiler_4.0.3
## [83] beeswarm_0.2.3 curl_4.3
## [85] png_0.1-7 interactiveDisplayBase_1.28.0
## [87] affyio_1.60.0 tibble_3.0.4
## [89] geneplotter_1.68.0 stringi_1.5.3
## [91] lattice_0.20-41 ProtGenerics_1.22.0
## [93] Matrix_1.2-18 vctrs_0.3.4
## [95] pillar_1.4.6 lifecycle_0.2.0
## [97] BiocManager_1.30.10 data.table_1.13.2
## [99] bitops_1.0-6 httpuv_1.5.4
## [101] rtracklayer_1.50.0 hwriter_1.3.2
## [103] R6_2.5.0 latticeExtra_0.6-29
## [105] affy_1.68.0 bookdown_0.21
## [107] promises_1.1.1 gridExtra_2.3
## [109] vipor_0.4.5 dichromat_2.0-0
## [111] MASS_7.3-53 assertthat_0.2.1
## [113] openssl_1.4.3 withr_2.3.0
## [115] GenomeInfoDbData_1.2.4 hms_0.5.3
## [117] rpart_4.1-15 coda_0.19-4
## [119] biovizBase_1.38.0 bbmle_1.0.23.1
## [121] base64enc_0.1-3 numDeriv_2016.8-1.1
## [123] shiny_1.5.0参考文献
Anders，Simon和Wolfgang Huber。2010年。“序列计数数据的差异表达分析。” 基因组生物学11（10）：R106 +。https://doi.org/10.1186/gb-2010-11-10-r106。
Anders，Simon，Paul T. Pyl和Wolfgang Huber。2015年。“ HTSeq –用于处理高通量测序数据的Python框架。” 生物信息学31（2）：166–69。https://doi.org/10.1093/bioinformatics/btu638。
Benjamini，Yoav和Yosef Hochberg。1995年。“控制错误发现率：一种实用且强大的多重测试方法。” 皇家统计学会杂志。系列B（方法论） 57（1）：289–300。http://www.jstor.org/stable/2346101。
R. Bourgon，R。Gentleman和W. Huber。2010年。“独立过滤提高了高通量实验的检测能力。” 美国国家科学院院刊107（21）：9546–51。https://doi.org/10.1073/pnas.0914005107。
Bray，Nicolas，Harold Pimentel，Pall Melsted和Lior Pachter。2016年。“接近最佳概率的Rna-Seq量化”。自然生物技术34：525–27。http://dx.doi.org/10.1038/nbt.3519。
Dudoit，Rine，Yee H.Yang，Matthew J.Callow和Terence P.Speed。2002年。“在复制的cDNA微阵列实验中鉴定差异表达基因的统计方法。” 统计学报，111–39。
Frankish，Adam，Alexandra Bignell，Andrew Berry，Andrew Yates，Anne Parker，Bianca M Schmitt，Bronwen Aken等。2018年。“人类和小鼠基因组的GENCODE参考注释。” 核酸研究47（D1）：D766–D773。
Hardcastle，Thomas和Krystyna Kelly。2010。“ baySeq：用于识别序列计数数据中差异表达的经验贝叶斯方法。” BMC生物信息学11（1）：422+。https://doi.org/10.1186/1471-2105-11-422。
Himes，Blanca E.，江晓峰，Peter Wagner，胡若希，王奇玉，Barbara Klanderman，Reid M.Whitaker等。2014年。“ RNA-Seq转录组分析将CRISPLD2确定为糖皮质激素响应基因，可调节气道平滑肌细胞的细胞因子功能。” 一九（6）。https://doi.org/10.1371/journal.pone.0099625。
Huber，Wolfgang，Vincent J. Carey，Robert Gentleman，Simon Anders，Marc Carlson，Benilton S. Carvalho，Hector Corrada C. Bravo等。2015年。“使用Bioconductor安排高通量基因组分析。” 自然方法12（2）：115–21。https://doi.org/10.1038/nmeth.3252。
Huntley，Melanie A.，Jessica L. Larson，Christina Chaivorapol，Gabriel Becker，Michael Lawrence，Jason A.Hackney和Joshua S.Kaminker。2013年。“ ReportingTools：用于高通量基因组分析的自动化结果处理和演示工具包。” 生物信息学29（24）：3220–1。https://doi.org/10.1093/bioinformatics/btt551。
Ignatiadis，Nikolaos，Bernd Klaus，Judith Zaugg和Wolfgang Huber。2016年。“数据驱动的假设加权增加了基因组规模多重测试的检测能力。” 自然方法。http://dx.doi.org/10.1038/nmeth.3885。
科斯特，约翰内斯和斯文·拉曼。2012年。“ Snakemake-可扩展的生物信息学工作流引擎。” 生物信息学。https://doi.org/10.1093/bioinformatics/bts480。
罗（Law），慈善机构（W. Charity W.），陈云顺（Yunshun Chen），史维（Wei Shi）和史密斯（Gordon K.Smyth）。2014年。“ Voom：精确的权重可解锁用于RNA-seq读数计数的线性模型分析工具。” 基因组生物学15（2）：R29 +。https://doi.org/10.1186/gb-2014-15-2-r29。
劳伦斯（Lawrence），迈克尔（Michael），沃尔夫冈·胡贝尔（Wolfgang Huber），埃尔维·帕杰斯（HervéPagès），帕特里克·阿博永（Patrick Aboyoun），马克·卡尔森（Marc Carlson），罗伯特·Gentleman，马丁·摩根（Martin T. 2013年。“用于计算和注释基因组范围的软件。” 由Andreas Prlic编辑。PLoS计算生物学9（8）：e1003118 +。https://doi.org/10.1371/journal.pcbi.1003118。
Leek，Jeffrey T.，2014年。“ svaseq：从测序数据中去除批处理效应和其他有害噪声。” 核酸研究42（21）：000 https://doi.org/10.1093/nar/gku864。
Leng，N.，JA Dawson，JA Thomson，V.Ruotti，AI Rissman，BMG Smits，JD Haag，MN Gould，RM Stewart和C.Kendziorski。2013年。“ EBSeq：用于推断RNA-seq实验的经验贝叶斯层次模型。” 生物信息学29（8）：1035–43。https://doi.org/10.1093/bioinformatics/btt087。
Leong，Hui S.，Keren Dawson，Chris Wirth，Yaoyong Li，Yvonne Connolly，Duncan L.Smith，Caroline R.Wilkinson和Crispin J.Miller。2014年。“全球非编码RNA系统响应压力调节裂变酵母蛋白水平。” 自然通讯5. https://doi.org/10.1038/ncomms4947。
李波和科林·杜威。2011年。“ RSEM：在有或没有参考基因组的情况下，从RNA-Seq数据进行准确的转录本定量。” BMC生物信息学12：323+。https://doi.org/10.1186/1471-2105-12-3231。
廖扬，GK史密斯和史威。2014年。“ featureCounts：一个高效的通用程序，用于将序列读取分配给基因组特征。” 生物信息学30（7）：923-30。https://doi.org/10.1093/bioinformatics/btt656。
Love，Michael I.，John B.Hogenesch和Rafael A.Irizarry。2016年。“ Rna-Seq片段序列偏倚的建模减少了转录本丰度估计中的系统错误。” 自然生物技术34（12）：1287–91。http://dx.doi.org/10.1038/nbt.3682。
Love，Michael I.，Wolfgang Huber和Simon Anders。2014年。“使用DESeq2对RNA-seq数据进行倍数变化和分散的适度估计。” 基因组生物学15（12）：550+。https://doi.org/10.1186/s13059-014-0550-8。
Love，Michael I.，Charlotte Soneson，Peter F.Hickey，Lisa K.Johnson，N.Tessa Pierce，Lori Shepherd，Martin Morgan和Rob Patro。2020年。“ Tximeta：用于鉴定RNA序列中来源的参考序列校验和。” PLOS计算生物学。https://doi.org/10.1371/journal.pcbi.1007664。
Patro，Rob，Geet Duggal，Michael I.Love，Rafael A.Irizarry和Carl Kingsford。2017年。“ Salmon提供转录物表达的快速且有偏差意识的定量。” 自然方法。http://dx.doi.org/10.1038/nmeth.4197。
Risso，Davide，John Ngai，Terence P.Speed和Sandrine Dudoit。2014年。“使用对照基因或样品的因子分析对RNA-seq数据进行标准化。” 自然生物技术32（9）：896–902。https://doi.org/10.1038/nbt.2931。
罗伯特，克里斯特尔和米克·沃森。2015年。“ RNA-Seq定量错误影响与人类疾病相关的基因。” 基因组生物学。https://doi.org/10.1186/s13059-015-0734-x。
医学博士罗宾逊，麦卡锡DJ和史密斯（GK Smyth）。2009年。“ edgeR：用于数字基因表达数据差异表达分析的Bioconductor软件包。” 生物信息学26（1）：139–40。https://doi.org/10.1093/bioinformatics/btp616。
Schurch，Nicholas J.，Pieta Schofield，Marek Gierlinski，Christian Cole，Alexander Sherstnev，Vijender Singh，Nicola Wrobel等。2016年。“ Rna-Seq实验需要多少次生物复制，您应该使用哪种差异表达工具？” 22（6）：839-51。https://doi.org/10.1261/rna.053959.115。
Soneson，Charlotte，Michael I.Love和Mark Robinson。2015年。“ RNA序列差异分析：转录本水平的估计改善了基因水平的推论。” F1000研究4（1521）。https://doi.org/10.12688/f1000research.7563.1。
Tonner，Peter D，Cynthia L Darnell，Barbara E Engelhardt和Amy K Schmid。2017年。“通过高斯过程回归检测微生物种群的差异增长。” 基因组研究27：320-33。https://doi.org/10.1101/gr.210286.116。
Townes，F。William，Stephanie C. Hicks，Martin J. Aryee和Rafael A. Irizarry。2019年。“基于多项模型的单细胞Rna-Seq的特征选择和降维。” bioRxiv。https://doi.org/10.1101/574574。
Trapnell，Cole，David G Hendrickson，Martin Sauvageau，Loyal Goff，John L Rinn和Lior Pachter。2013。“利用RNA-seq进行转录分辨率的基因调控差异分析。” 自然生物技术。https://doi.org/10.1038/nbt.2450。
威克姆，哈德利。2009年GGPLOT2。纽约，纽约：纽约斯普林格。https://doi.org/10.1007/978-0-387-98141-3。
Witten，Daniela M.，2011年。“使用泊松模型对测序数据进行分类和聚类。” 应用统计年鉴5（4）：2493–2518。https://doi.org/10.1214/11-AOAS493。
吴浩，王驰和吴志金。2013年。“一种用于分散的新的收缩率估算器改善了RNA-seq数据中差异表达的检测。” 生物统计学14（2）：232–43。https://doi.org/10.1093/biostatistics/kxs033。
朱安琪，约瑟夫·易卜拉欣和迈克尔·爱。2018年。“序列计数数据的重尾先验分布：消除噪声并保留较大差异。” bioRxiv。https://doi.org/10.1101/303255。
RNA-seq工作流程：基因水平的探索性分析和差异表达相关推荐
- android的构成和工作流程,分析Android中View的工作流程
8种机械键盘轴体对比 本人程序员,要买一个写代码的键盘,请问红轴和茶轴怎么选? 在分析View的工作流程时,需要先分析一个很重要的类,MeasureSpec.这个类在View的测量(Measure)过 ...
- 计算机网络---IP数据报组成计及IP模块工作流程
IP数据报组成 IP数据报是由首部与数据部两部分组成. 首部的前一部分是固定的20个字节,后一部分是可选的,是每个数据报必不缺少的一部分.后面的数据部分是来自于传输层. 这是IP数据报首部的结构图,上 ...
- Android 7.0 Vold工作流程
一.Vold工作机制 Vold是Volume Daemon的缩写,它是Android平台中外部存储系统的管控中心,是管理和控制Android平台外部存储设备的后台进程.其功能主要包括:SD卡的插拔事件 ...
- JMP 17 全新发布，让分析工作流程更简便
虽然工艺工程师.技术创新者.科学家拥有不同的专业领域知识,但他们都依赖JMP系列统计发现软件来应对现代数据所带来的独特挑战. JMP 产品家族最新版本 - JMP 17.JMP Pro 17 及JMP ...
- 一文掌握RNA seq，RNA seq课程大汇总
RNA测序(RNA-seq)在过往十年里逐渐成为全转录组水平分析差异基因表达和研究mRNA差异剪接必不可少的工具.RNA-seq帮助大家对RNA生物学的理解会越来越全面:从转录本在何时何地转录到RNA ...
- 神经肿瘤组学基础知识、工作流程及应用
在过去几年中,从脑肿瘤患者中采集的具有常规临床用途的神经影像数据的数量.种类和复杂性以及由此产生的影像参数数量都大幅增加.因此,如果没有人工智能(AI)领域的方法支持,对成像数据进行及时且成本效益高的 ...
- iMeta | 大连海洋大学傅松哲和根特大学杨倩开发宏基因组测序和流式细胞术相结合的工作流程...
点击蓝字 关注我们 宏基因组测序结合流式细胞术为城市污水中细菌病原体的微生物风险评估提供了新的框架 原文链接DOI: https://doi.org/10.1002/imt2.77 COMMENTAR ...
- 笔记：机器学习的工作流程
1.机器学习工作流程介绍 本文学习如何思考机器学习问题.机器学习(ML)不仅关于数学和建模,而且关于选择问题的设置,识别客户需求和长期目标.本文的组织如下: 1.通过识别关键利益相关者,选择正确的度量 ...
- QIIME 2教程. 02插件工作流程Plugin Workflows(2021.2)
QIIME 2插件工作流程概述 Overview of QIIME 2 Plugin Workflows 英文原文见:https://docs.qiime2.org/2021.2/tutorials/ ...
最新文章
- WebSocket实战之————Workerman服务器的安装启动
- 谭浩强课后题之----求最大公约数和最小公倍数
- 周报速递丨百度与美团等互联互通；京东Shopify达成合作
- 神奇的sstream头文件(整型与字符串自由转换)
- HDU 3641 Treasure Hunting(阶乘素因子分解+二分)
- 乌云挂了，知识库的文章却在流传
- from Crypto.Cipher import AES报错
- [转] 虚拟机VMware3种网络模式（桥接、nat、Host-only）的工作原理
- Eclipse启动时总是提示“subversive connector discovery”解决方案
- CH372或CH375的USB通讯问题解答
- Hybird App 应用开发中5个必备知识点复习
- 【脑图制作】万彩脑图大师教程 | 怎么制作思维导图
- UVALive 7269 Snake Carpet
- pip install清华镜像源
- 面试官：生产环境中 CPU 利用率飙高怎么办？
- 阿里云OSS上传图片、PDF设置链接预览
- 计算机控制面板属性在哪里打开,win7如何打开系统属性面板|win7打开系统属性面板的方法...
- 系统集成项目管理工程师10大管理
- android谷歌手机刷机教程
- 开源项目CIIP(企业信息管理系统框架).2018.1.0910版更新介绍-上周工作总结