内容

作者:Sonali Arora (sarora@fredhutch.org）
时间:2015年7月20-22日

本课程的材料要求R版本3.2.1和Bioconductor版本3.2

0．1生物导体rna序列分析高级实验室

本实验室将引导您通过端到端RNA-Seq差异表达工作流程，使用DESeq2与其他Bioconductor包。

注:其他若干Bioconductor包也可以用于基因水平差异表达的统计推断，包括刨边机，BaySeq，DSS而且limma．

0.2锻炼

使用来自气道、设计及实施
端到端RNA-Seq差异表达分析，使用DESeq2

步骤包括-

0.3分析数据

本实验室使用的数据是经地塞米松(一种具有抗炎作用的合成糖皮质激素类固醇)处理的气道平滑肌细胞的RNA-Seq实验。例如，糖皮质激素被用于哮喘患者，以预防或减轻呼吸道炎症。在实验中，四株原代人气道平滑肌细胞系用1微摩尔地塞米松处理18小时。对于每一种细胞系，我们都有一个处理过的和一个未处理过的样本。

实验参考如下:

Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q.“RNA-Seq转录组分析发现CRISPLD2是调节气道平滑肌细胞细胞因子功能的糖皮质激素响应基因。”公共科学图书馆，2014年6月13日;9(6):e99625。PMID:24926665．地理:GSE52778．

对于我们的分析，我们将使用数据包中的数据气道．

库(“气道”)数据(气管)

0.4解决方案

0.4.1答案1:加载数据

存储在里面的数据气道是一个SummarizedExperiment对象。

库(“气道”)数据(气道)se <-气道se
##类:rangedsummarizeexperimental ## dim: 64102 8 ##元数据(1):“## assays(1):计数## rownames(64102): ENSG00000000003 ENSG00000000005…LRG_98 LRG_99 ## rowRanges元数据列名(0):## colnames(8): SRR1039508 SRR1039509…SRR1039520 SRR1039521 ## colData names(9): SampleName cell…样本BioSample

0.4.2答案2:创建DESeqDataSet

一旦我们完成了完整的注释SummarizedExperiment对象，我们可以构造一个DESeqDataSet对象，然后它将形成实际的起点DESeq2包，将在下面的部分中描述。我们为分析添加了适当的设计。

库("DESeq2") dds <- DESeqDataSet(se, design = ~ cell + dex)

0.4.3答案3:差分表达式管道

这样做很方便untrt第一层是在敏捷因子，这样默认的log2倍的变化会被计算为处理后的变化(默认R将选择第一个字母级别，记住:计算机不知道该做什么，除非你告诉它们)。这个函数relevel实现:

Dds $dex <- relevel(Dds $dex， "untrt")

最后，我们准备运行差分表达式管道。准备好数据对象后，DESeq2现在可以通过对函数的一次调用来运行分析DESeq：

dds <- DESeq(dds)
##估计大小因素##估计分散##基因分散估计##均值-分散关系##最终分散估计##拟合模型和测试

这个函数将为它执行的各个步骤打印一条消息。这些在手册页中有更详细的描述DESeq，可通过键入来访问DESeq ?．简单地说，这些是:大小因素的估计(控制测序实验文库大小的差异)，每个基因离散度的估计，以及拟合广义线性模型。

一个DESeqDataSet返回，其中包含所有拟合信息，下一节描述如何从该对象中提取感兴趣的结果表。

0.4.4答案4:构建结果表

调用结果没有任何参数将提取估计的log2倍变化和p设计公式中最后一个变量的值。如果这个变量有2个以上的级别，结果将提取结果表，以进行最后一级与第一级的比较。

(res <- results(dds))
## log2 fold change (MAP): dex trt vs untrt Wald测试p-value:eng00000000003 708.60217 -0.37424998 0.09873107 -3.7906000 0.0001502838 0.001251416 ## ENSG00000000005 0.00000 NA NA NA NA NA NA NA ## eng00000000419 520.29790 0.20215551 0.10929899 1.8495642 0.064373763851 0.13684258 0.2648902 0.7910940556 0.910776144 ## eng0000000046057.93263 -0.08523371 0.24654402 -0.3457140 0.7295576905 0.878646940  ## ... ... ... ... ... ... ...## lrg_94 0 na na na na na na ## lrg_96 0 na na na na na na ## lrg_97 0 na na na na na na ## lrg_98 0 na na na na na ## lrg_99 0 na na na na na na na

作为res是一个DataFrame对象时，它携带包含列含义信息的元数据:

mcols (res use.names = TRUE)
## 6行2列的数据帧##类型描述## <字符> <字符> ## baseMean所有样本归一化计数的中间平均值## log2FoldChange结果log2 fold change (MAP): dex trt vs untrt# # lfcSE结果标准错误:dex trt vs untrt# #统计结果Wald统计:dex trt vs untrt# # pvalue结果Wald测试p-value: dex trt vs untrt# # padj结果BH调整p值

第一列，baseMean，是对所有样本进行归一化计数值的平均值，除以大小因子。其余四列是指具体的对比，即比较泰爱泰党水平高于untrt水平为因素变量敏捷．参见帮助页结果(通过输入结果呢?)，以获取其他对比资料。

列log2FoldChange是效应量估计。它告诉我们，与未处理的样本相比，由于使用地塞米松治疗，基因的表达似乎发生了多大的变化。该值以2为基数的对数刻度报告:例如，1.5的log2倍变化意味着基因的表达增加了一个乘法因子\(2^{1.5} \约2.82\)．

我们还可以用下面这行代码总结结果，它报告了一些附加信息

总结(res)
## ##出33469非零总读计数##调整p值< 0.1 ## LFC > 0(上):2641,7.9% ## LFC < 0(下):2242,6.7% ##异常值[1]:0,0% ##低计数[2]:15441,46% ##(平均计数< 5)##[1]见'cooksCutoff'参数?results ##[2]见'independentFiltering'参数?results

0.4.5答案5:可视化结果

一种快速可视化特定基因计数的方法是使用plotCounts函数作为参数DESeqDataSet，一个基因名称，以及用来绘制计数图的组。

topGene <- rownames(res)[it .min(res$padj)] plotCounts(dds, gene=topGene, intgroup=c("dex"))

0．5参考文献

有关更详细的分析，请参阅
-案例研究-如何建立一个总结实验气道数据集
-差异表达实验室

0.6在BioC2015上不要错过什么!

如果你喜欢这个实验室，想要了解更多这方面的知识，不要错过以下的BioC2015实验室

0.7sessionInfo ()

sessionInfo ()
## R版本3.2.1(2015-06-18)##平台:x86_64-unknown-linux-gnu(64位)##运行在Ubuntu 14.04.2 LTS下## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC=C LC_TIME=en_US。UTF-8 ## [4] LC_COLLATE=C LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME=C LC_ADDRESS= c# # [10] lc_phone =C LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4并行统计图形grDevices utils数据集方法基础## ##其他附加包:# # # # [1] RColorBrewer_1.1-2 gplots_2.17.0 [3] DESeq2_1.9.23 RcppArmadillo_0.5.200.1.0 # # [5] Rcpp_0.11.6 airway_0.103.1 # # [7] Rattus.norvegicus_1.3.1 org.Rn.eg.db_3.1.2 # # [9] GO.db_3.1.2 OrganismDbi_1.11.42 # # [11] BSgenome.Rnorvegicus.UCSC.rn5_1.4.0 BSgenome_1.37.3 # # [13] rtracklayer_1.29.12 TxDb.Rnorvegicus.UCSC.rn5.refGene_3.1.3 # # [15] org.Hs.eg.db_3.1.2 RSQLite_1.0.0 # # [17] DBI_0.3.1 TxDb.Hsapiens.UCSC.hg19.knownGene_3.1.3 # # [19] GenomicFeatures_1.21.13 AnnotationDbi_1.31.17 # # [21][23] GenomicAlignments_1.5.11 Rsamtools_1.21.14 ## [25] Biostrings_2.37.2 XVector_0.9.1 ## [27] SummarizedExperiment_0.3.2 Biobase_2.29.1 ## [31] iranges_21.16 GenomeInfoDb_1.5.8 ## [31] IRanges_2.3.14 S4Vectors_0.7.10 ## [33] BiocGenerics_0.15.3 ggplot2_1.0.1 ## [35] BiocStyle_1.7.4 ## ##通过命名空间加载(并且没有附加):## [1] bitops_1.0-6 hty_0.12.1 rpart_4.1-10 ## [7] Hmisc_3.16-0 colorspace_1.2-6 nnet_7.3-10 ## [10] gridExtra_2.0.0 curl_0.9.1 graph_1.47.2 ## [13] formatR_1.2 labeling_0.3 caTools_1.17.1 ## [19] string_1 .0.0 digest_0.6.8 foreign_0.8-65 ## [22] rmarkdown_0.7 htmltools_0.2.6 BiocInstaller_1.19.8 ## [28] acepack_1. 3.0 RCurl_1.95-4.7[40] grid_3.2.1 gdata_2.17.0 lattice_0.20-33 ## [43] splines_3.2.1 annotate_1.47.1 locfit_1.5-9.1 ## [46] knitr_1.10.5 geneplotter_1.47.0 reshape2_1.4.1 ## [52] xml_98 -1.3 evaluate_0.7 latticeExtra_0.6-26 ## [55] lambda.r_1.1.7 httpuv_1.3.2 gtable_0.1.2 ## [58] mime_0.3xtable_1.7-4 survivval_2 .38-3 ## [61] cluster_2.0.2 interactiveDisplayBase_1.7.0