介绍

这个研讨会使用了来自一个涉及气道平滑肌细胞的基因水平RNA-seq实验的数据;附带的小插图提供了详细信息气道包中。原始数据来自Himes等人，“RNA-Seq转录组分析鉴定CRISPLD2是一种糖皮质激素应答基因，可调节气道平滑肌细胞中的细胞因子功能。”公共科学图书馆。2014年6月13日;9(6):e99625。PMID:24926665地理:GSE52778。从摘要：“使用RNA-SEQ [......]我们在用地塞米松治疗的四个一次人体ASM细胞系中表征转录组变化 - 一种有效的合成糖皮质激素（1微摩尔18小时）。”

在排序，读取征收和对齐的概要之后，我们加入分析，并将对齐的读取到每个样本中的读取重叠区域（基因）的读取区域的计数表。我们专注于实验的子集，其中4种细胞系各种细胞系用地塞米松或对照治疗。

设置（工作室期间不是必需的）

我们将使用来自生物体。这些依赖于其他几个包。软件包及其依赖项已安装在课程中使用的Amazon Machine实例（AMI）上。在课程后为您自己的计算机，安装包

源(“//www.andersvercelli.com/biocLite.R”)biocLite (c(“DESeq2”、“org.Hs.eg.db”))

安装需要每台计算机执行一次，而不是每次使用包都会执行一次。

我们在分析中使用两个数据文件。数据文件在于Sydney2016github库。安装github存储库

Bioclite（“Bioconductor / Syndey2016”）

安装了包后，文件的位置（用于使用file.choose ()，下)由

系统。文件(包= " Sydney2016”、“extdata”)

数据输入

第一个挑战是输入数据。我们从“表型”数据开始，描述实验中使用的样本。数据是一个简单的8行几列的表;它可以在Excel中创建，并导出为以tab分隔的文件。找到这个文件在Amazon机器实例中的位置

colDataFile <- file. choose() #查找“airway-colData.tab”

并将数据读入R使用read.table ()功能。数据足够小以便在其中看R会话，输入变量名)或inRStudio（通过使用视图()或者双击“Environment”选项卡中的变量)。

colData <- read.table(colDataFile) colData
## 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 ## SRR1039516GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353 ## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354 ## SRR1039520 GSM1275874 N061011 untrt SRR1039520 101 SRX384357 ## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358 ## Sample bisample # SRR1039508 SRS508568 SAMN02422669 ## SRR1039509 SRS508567SAMN02422675 ## SRR1039512 SRS508571 srs02422678 SRS508572 srs02422670 SRR1039516 SRS508575 srs02422682 SRR1039517 SRS508576 srs02422673 srs1039520 SRS508579 srs02422683 SRS508580 srs02422677

这应该会顺利进行;对于真实的数据，人们通常需要花费相当多的时间来调整论点read.table ()说明页眉、行名、注释等的存在。

下一个挑战是输入表达式估计数。这是一个矩阵，行表示感兴趣的区域(基因)，列表示样本。矩阵中的条目是每个样本中每个区域重叠的读取数。重要的是，这些值是原始计数，而不是FPKM之类的可伸缩度量。找到文件

- file. choose() #查找“airway-assay.tab”

输入数据并使用头()查看数据的前几行。

<- read.table(assayFile)头(化验)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 ## ENSG00000000003 679 448 873 408 1138 ## 467 ENSG00000000419 515 621 365 587 ## 260 ENSG00000000457 211 263 164 245 ## 60 ENSG00000000460 55 40 35 78 ## ENSG00000000938 0 0 2 0 1 ## ENSG000000009713251 7606/127/127 76/127/127/127/127/797/797/127/797/127/127/797

探索

计算每个样本的“库大小”(映射读取的总数量)colSums ()

colSums(化验)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 srr1039571 18809481 25348649 15163415 24448408 30818215 ## SRR1039520 SRR1039521 ## 19126151 21164133

使用以下一系列命令创建每个基因的平均asinh转化的密度图(asinh是类似于log的，期望接近于零)读取计数。

绘图（密度（Rowmeans（AsinH（AsinH）））））

多维标度(MDS)是一种降维方法，将n空间中的向量投影到两个(或更多)维中。使用dist（）函数来计算每个样本之间的(欧几里得)距离和基数R函数cmdscale ()对距离矩阵执行MDS。我们可以使用图()将结果可视化，并看到8个样本中每个样本的大致位置。使用参数col根据单元格线(美元colData细胞)或实验处理美元colData敏捷。

d <- dist(t(asinh(化验)))plot(cmdscale(d)， pch=19, cex=2)

plot(cmdscale(d)， pch=19, cex=2, col=colData$cell)

绘图（cmdscale（d），pch = 19，cex = 2，col = coldata $ dex）

注意，细胞系之间是相对相似的。这表明细胞系在后续分析中应作为协变量使用。

微分表达式

我们将使用DESeq2差异表达分析包装;值得注意的是，其他选择是可能的edger.和limma。

通过提供表达计数数据，实验说明和“模型”来开始分析，描述我们想要估算的统计关系。对于我们的模型和基于前一部分的探索性分析，我们假设计数是确定的细胞系和地塞米松治疗。我们主要包括细胞系作为协变量;我们的主要兴趣是在地塞米松。

图书馆（DESEQ2）DDS < -  DESQDATASETFRommatrix（测定，COLDATA，〜CELL + DEX）

分析非常容易调用，但计算涉及许多复杂的统计问题，包括:

代码被调用为:

dds < - DESeq (dds)
##估计大小因子
# #估计分散
##基因明智色散估计
##均值分散关系
##最终离散度估计
拟合模型及测试
dds
##类：DESQDATASET ## DIM：33469 8 ## METADATA（1）：版本##测定（3）：COUNTS MU COOKS ## Rowname（33469）：ENSG00000000419 ... ## ensg00000273493 ## rowdata名称（46）：底座束缚...偏差MaxCooks ## Colnames（8）：SRR1039508 SRR1039509 ... SRR1039509 ... SRR1039520 SRR1039521 ## COLDATA名称（10）：SAMPLENAME CELL ... BIOSAMPLE SOUDITFACTOR

的DESeq ()函数返回一个对象，可以作为进一步分析的起点，例如生成差异表达基因的“顶表”，按调整顺序排序(用于多次比较)。P价值观。

结果< - 结果（DDS）结果
## log2折叠更改（地图）：DEX UNTRT VS TRT ## WALD TEST P值：DEX UNTRT VS TRT ##与33469行和6列## BASEMEAL LOG2FOLDCHANGE LFCSE STAT ##    <数值> ## ENSG00000000003 708.6021697 0.37415246 0.09884435 3.7852692 ## ENSG00000000419 520.2979006 -0.20206175 0.10974241 -1.8412367 ## ENSG00000000457 237.1630368 -0.03616686 0.13834540 -0.2614244 ## ENSG00000000460 57.9326331 0.08445399 0.24990709 0.3379415 ## ENSG00000000938 0.3180984 0.08413901 0.15133424 0.5559813 ## ... ..。... ... ... ## ENSG00000273487 8.1632350 -0.55007238 0.3725061 -1.4766803 ## ENSG00000273488 8.5844790 -0.10513006 0.3683837 -0.2853820 ## ENSG00000273489 0.2758994 -0.06947899 0.1512520 -0.4593591 ## ENSG00000273492 0.1059784 0.02314357 0.1512520 0.1530133 ## ENSG00000273493 0.1061417 0.02314357 0.15125200.1530133 ## PVALUE PADJ ##   ## ENSG00000000003 0.011279829 0.011279829 0.01279829 0.01279829 0.01279829 0.01279829 0.01279829 0.01279829 0.01279829 0.01279829 ## ESG0000000030 0.1901279829＃937652416 0.912936622 ## ENSG00000000460 0.7354072415 0.883203048 ## ENSG00000000938 0.5782236287 NA ## ...... ## ENSG00000273487 0.1397614 0.3376252 ## ENSG00000273488 0.7753515 0.9032267 ## ENSG00000273489 0.6459763 NA ## ENSG00000273492 0.8783878 NA ## ENSG00000273493 0.8783878 NA
<- head(order(result$padj)， 10) top = result[ridx，] top
## log2 fold change (MAP): dex untrt vs trt ## Wald test p-value:敏捷untrt vs泰爱泰党# # DataFrame 10行6列# # baseMean log2FoldChange lfcSE stat # # <数字> <数字> <数字> <数字> # # ENSG00000152583 # # ENSG00000165995 495.0929 -3.186823 0.12815654 -25.06116 997.4398 -4.313962 0.17213733 -24.86664 # # ENSG00000101347 12703.3871 -3.618734 0.14894336 -24.29604 3409.0294 # # ENSG00000120129## ENSG00000157214 12285.6151 -3.553360 0.15798211 -22.49217 ## ENSG00000157214 3009.2632 -1.948723 0.08867432 -21.97618 ## ENSG00000162614 5393.1017 -2.003487 0.09269629 -21.61345 ## ENSG00000125148 3656.2528 -2.167122 0.10354724 -20.92882 ## ENSG0000015473430315.1355 -2.286778 0.11308412 -20.22192 ## pvalue padj ##   ## ENSG00000152583 1.319237e-138 2.360906e-134 ## ENSG00000165995 1.708565e-136 1.528824e-132 ## ENSG00000101347 2.158637e-130 1.287699e-126 ## ENSG00000120129 2.937247e-130 1.314124e-126 ## ENSG00000189221 1.656535e-123 5.929070e-120 ## ENSG00000211445 4.952260e-112 1.477094e-108 ## ENSG00000157214 4.867315e-107 1.244364e-103 ## ENSG00000162614 1.342345e-103 3.002826e-100 ## ENSG00000125148 2.926277e-97 5.818739e-94 ## ENSG00000154734 6.278709e-91 1.123638e-87

理解

将统计结果放入生物学环境中有很多机会。初始步骤是映射用于标记感兴趣区域的密码联合集基因标识符，以更具食物HGNC基因符号。为此，我们使用org.Hs.eg.db包，一个例子生物体“注释”包包含来自公共领域资源的策划数据，每半年更新一次。的mapId ()函数在标识符类型之间映射，在我们的例子中映射到来自集合id的符号基因id。我们将这些数据添加到顶级表的结果中，以便它们可以与统计结果一起处理。

库（org.hs.eg.db）top $ symbol < -  mapids（org.hs.eg.db，rownames（top），“符号”，“ensembl”）
##'select（）'返回1：1键之间映射键和列之间
最佳
## log2 fold change (MAP): dex untrt vs trt ## Wald test p-value:敏捷untrt vs泰爱泰党# # DataFrame 10行和7列# # baseMean log2FoldChange lfcSE stat # # <数字> <数字> <数字> <数字> # # ENSG00000152583 # # ENSG00000165995 495.0929 -3.186823 0.12815654 -25.06116 997.4398 -4.313962 0.17213733 -24.86664 # # ENSG00000101347 12703.3871 -3.618734 0.14894336 -24.29604 3409.0294 # # ENSG00000120129## ENSG00000157214 12285.6151 -3.553360 0.15798211 -22.49217 ## ENSG00000157214 3009.2632 -1.948723 0.08867432 -21.97618 ## ENSG00000162614 5393.1017 -2.003487 0.09269629 -21.61345 ## ENSG00000125148 3656.2528 -2.167122 0.10354724 -20.92882 ## ENSG0000015473430315.1355 -2.286778 0.11308412 -20.22192 ## pvalue padj Symbol ##    ## ENSG00000152583 1.319237e-138 2.360906e-134 SPARCL1 ## ENSG00000165995 1.708565e-136 1.528824e-132 CACNB2 ## ENSG00000101347 2.158637e-130 1.287699e-126 SAMHD1 ## ENSG00000120129 2.937247e-130 1.314124e-126 DUSP1 ## ENSG00000189221 1.656535e-123 5.929070e-120 MAOA ## ENSG00000211445 4.952260e-112 1.477094e-108 GPX3 ## ENSG00000157214 4.867315e-107 1.244364e-103 STEAP2 ## ENSG00000162614 1.342345e-103 3.002826e-100 NEXN ## ENSG00000125148 2.926277e-97 5.818739e-94 MT2A ## ENSG00000154734 6.278709e-91 1.123638e-87 ADAMTS1

再现性

这里的计算通过报告分析中使用的软件版本，如下所示：

sessioninfo（）
## R版本3.3.1补丁(2016-10-12 r71512) ##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 16.04.1 LTS ## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# ## [3] LC_TIME=en_US。utf - 8 LC_COLLATE = en_US。UTF-8 ## [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] parallel stats4 stats graphics grDevices utils datasets ## [8] methods base ## ##其他attached packages:# # # # [1] org.Hs.eg.db_3.4.0 AnnotationDbi_1.35.5 [3] DESeq2_1.13.16 SummarizedExperiment_1.3.82 # # [5] Biobase_2.33.4 GenomicRanges_1.26.1 # # [7] GenomeInfoDb_1.10.0 IRanges_2.8.0 # # [9] S4Vectors_0.12.0 BiocGenerics_0.20.0 # # [11] BiocInstaller_1.24.0 # # # #通过加载一个名称空间(而不是附加):# # # # [1] Rcpp_0.12.7 formatR_1.4 RColorBrewer_1.1-2 [4] plyr_1.8.4 XVector_0.14.0 bitops_1.0-6 # # [7] tools_3.3.1 zlibbioc_1.20.0 rpart_4.1-10 # # [10] digest_0.6.10 RSQLite_1.0.0 annotate_1.51.1 # # [13] evaluate_0.10 tibble_1.2 gtable_0.2.0 # # [16] lattice_0.20-34 matrix_1.2 - 7.1 DBI_0.5-1 # # [19] yaml_2.1.13 gridExtra_2.2.1 genefilter_1.55.2 # #[22] stringr_1.1.0 knitr_1.14 cluster_2.0.5 # # [25] locfit_1.5 - 9.1 nnet_7.3-12 grid_3.3.1 # # [28] data.table_1.9.6 xml_3.98 - 1.4 survival_2.39-5 # # [31] BiocParallel_1.8.1 foreign_0.8 - 67 rmarkdown_1.1 # # [34] latticeExtra_0.6-28 Formula_1.2-1 geneplotter_1.51.0 # # [37] ggplot2_2.1.0 magrittr_1.5 Hmisc_3.17-4 # # [40] scales_0.4.0 htmltools_0.3.5splines_1 .3.1 ## [43] assertthat_0.1 xtable_1.8-2 colorspace_1.2-7 ## [46] stringi_1.1.2 acepack_1.3-3.3 RCurl_1.95-4.8 ## [49] munsell_0.4.3 chron_2.3-47