RNA-seq工作流程：基因水平的探索性分析和差异表达_标书中对于rna-seq差异转录本下游蛋白的验证

1.1实验数据

该工作流程中使用的数据存储在气道包装中，该数据 总结了RNA-seq实验，其中用地塞米松（一种具有抗炎作用的合成糖皮质激素类固醇）处理了气道平滑肌细胞 （Himes等，2014）。糖皮质激素被例如哮喘患者用来减轻气道炎症。在该实验中，用1微摩尔地塞米松处理了四个主要的人气道平滑肌细胞系18小时。对于这四个细胞系中的每一个，我们都有一个处理过的样品和一个未处理过的样品。有关实验的更多说明，请参见 PubMed条目24926665 ，有关原始数据，请参见 GEO条目GSE52778。

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）。

接下来，我们将描述由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 Dexamethasone

2.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使用的数据对象。

<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

<span style="color:#333333"><span style="color:#333333"><code>gse$cell <- gse$donor
gse$dex <- gse$condition</code></span></span>

<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>

<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

<span style="color:#333333"><span style="color:#333333"><code>gse$dex <- relevel(gse$dex, <span style="color:#dd1144">"untrt"</span>)</code></span></span>

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.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.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         NA

There 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] 31604

Now, 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-02