内容

1介绍

描述:本次研讨会将向您介绍Bioconductor收集用于统计分析和理解高通量基因组数据的R包。重点是数据探索，使用rna序列基因表达实验作为一个激励的例子。如何从R中访问常见序列数据格式?如何在分析中使用有关基因模型或基因注释的信息?数据的属性如何影响我应该执行的统计分析?我可以用R和Bioconductor？如何在R中处理非常大的数据集?这些都是本次研讨会要解决的问题。

要求:你需要自带笔记本电脑。研讨会将使用基于云的资源，所以你的笔记本电脑需要一个网络浏览器和WiFi功能。参加者应使用R而且RStudio对于本周早些时候介绍性研讨会中涉及的任务。基因表达生物学的一些知识和统计学第一课中所学的概念将是有帮助的。

相关性:本次研讨会与任何渴望探索基因组数据的人有关R．研讨会将帮助连接“核心”R处理数据的概念(例如，数据管理通过data.frame ()，统计模型lm ()或t.test ()，可视化使用图()或ggplot ())，以应对处理大型基因组数据集的特殊挑战。对于那些已经或将要拥有自己的基因组数据，并有兴趣更全面地了解如何使用这些数据的人来说，这将特别有帮助R．

1．1我们的目标

RNA-seq

  • 设计实验，例如，暴露于两种处理的4个细胞系的8个样本(基于Himes等人，PMID:24926665；详情请参阅气道包小插图)。
cell dex SRR1039508 N61311 untrt SRR1039509 N61311 trt SRR1039512 N052611 untrt SRR1039513 N052611 trt SRR1039516 N080611 untrt SRR1039517 N080611 trt SRR1039520 N061011 untrt SRR1039521 N061011 trt
  • 文库制备:mRNA到稳定双链DNA
  • “短”mrna衍生片段的DNA测序
  • 参照基因组或转录组的比对

来源:http://bio.lundberg.gu.se/courses/vt13/rnaseq.html

  • 最终结果:一个“计数”矩阵-读数与样本中的基因对齐。
Srr1039508 srr1039509 srr1039512 srr1039513 srr1039516 ensg00000000003 679 448 873 408 1138 ensg00000000005 00000 ensg00000000419 467 515 621 365 587 ... ... ... ... ... ...Srr1039517 srr1039520 srr1039521 ensg00000000003 1047 770 572 ensg00000000005 000 ensg00000000419 799 417 508 ... ... ... ...

研究问题

  • 哪种基因计数在地塞米松之间差异最大untrt而且泰爱泰党实验治疗?
  • 我们会试着去理解如何为了做到这一点，不需要统计细节。

我们的目标

  • 将20个表达差异最大的基因可视化为热图。

2数据收集、输入、表示和清理

2.1基地R数据结构

样品信息

  • 简单的“制表符分隔值”文本文件，例如，从Excel导出。
  • 使用基数输入R命令read.table ()
  • “原子”向量，例如:整型()
  • 因素()而且NA
  • data.frame ():协调管理
    • 列访问＄
    • 与子集[，]
Samples_file <- system。file(package="BiocIntro"， "extdata"， "samples.tsv") samples <- read.table(samples_file) samples ## cell dex avgLength ## SRR1039508 N61311 untrt 126 ## SRR1039513 N61311 untrt 126 ## SRR1039513 N052611 trt 87 ## SRR1039516 N080611 untrt 120 ## SRR1039517 N080611 trt 126 ## SRR1039520 N061011 untrt 101 ## SRR1039521 N061011 trt 98 samples$dex <- relevel(samples$dex， "untrt")

计数

  • 另一个tsv文件。很多行，所以使用头()要查看前几个。
  • 行名:基因标识符。列名:示例标识符。
  • 所有列都是相同的(数字)类型;表示为矩阵()而不是data.frame ()
Counts_file <- system。文件(包= " BiocIntro”、“extdata”,counts.tsv)计数< - read.table (counts_file)暗(计数)# #[1]63677 8头(计数)# # SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 # # ENSG00000000003 679 448 873 408 1138 # # ENSG00000000005 0 0 0 0 0 # # ENSG00000000419 467 515 621 365 587 260 211 263 164 245 # # ENSG00000000457 60 55 40 35 78 # # # # ENSG00000000460 ENSG00000000938 0 0 2 0 1 # # SRR1039517 SRR1039520 SRR1039521 # # # # ENSG00000000005 ENSG00000000003 1047 770 572 0 0 0 # # # # ENSG00000000457 ENSG00000000419 799 417 508331 233 229 ## ENSG00000000460 63 76 60 ## ENSG00000000938 000计数<- as.matrix(计数)

2．2基因组范围(农庄）

行注释。

  • “GTF”文件包含有关基因模型的信息。
  • 与本实验相关的GTF文件-相同的生物体(智人)，基因组(GRCh37)和基因模型注释(Ensembl release-75)用于对齐和计数步骤-是
url <- "ftp://ftp.ensembl.org/pub/release-75/gtf/homo_sapiens/Homo_sapiens.GRCh37.75.gtf.gz"
  • 使用BiocFileCache将资源一次下载到贯穿始终的位置R会话。
库(BiocFileCache)
gtf_file <- bfcrpath(rnames = url) ##使用临时缓存/var/folders/ n/gmsh_22s2c55v816r6d51fx1tnyl61/T//RtmpLUzYG5/BiocFileCache ##添加rname 'ftp://ftp.ensembl.org/pub/release-75/gtf/homo_sapiens/Homo_sapiens.GRCh37.75.gtf.gz'
库(rtracklayer) GTF <-导入(gtf_file)

  • 一个农庄对象

    • Range-specific信息
    • 每个范围上的注释
    • 使用函数访问核心元素:seqnames ()(例如,染色体),start ()/结束()/宽度()，链()等。
    • 使用＄或mcols()美元访问范围上的注释。
    • Bioconductor惯例:基于1的封闭间隔(如Ensembl)，而不是基于0的1/2开放间隔(如UCSC)。
  • 将信息过滤为基因级注释，只保留每个基因组范围的部分信息。使用gene_id列名称()．
rowidx <- gtf$type == "gene" colidx <- c("gene_id"， "gene_name"， "gene_biotype") genes <- gtf[rowidx, colidx] names(基因)<- genes$gene_id genes$gene_id <- NULL genes ## GRanges对象具有63677个范围和2个元数据列:## seqnames ranges strand | gene_name gene_biotype ##    |   ## ENSG00000223972 1 11869-14412 + | DDX11L1假基因## ENSG00000227232 1 14363-29806 - | WASH7P假基因## ENSG00000243485 1 29554-31109 + | MIR1302-10 lincRNA ## ENSG00000237613 1 34554-36081 - | FAM138A lincRNA ## ENSG00000268020 1 52473-54936 + | OR4G4P假基因## ... ... ... ... . ... ...## ENSG00000210194 MT 14674-14742 - | MT- te Mt_tRNA ## ENSG00000198727 MT 14747-15887 + | MT- cyb蛋白_coding ## ENSG00000210195 MT 15888-15953 + | MT- tt Mt_tRNA ## ENSG00000210196 MT 15956-16023 - | MT- tp Mt_tRNA ## ------- # seqinfo: 265个来自未指定基因组的序列;没有seqlengths

2.3协调管理(SummarizedExperiment）

三个不同的数据集

  • 计数: RNAseq工作流的结果
  • 样品:样品和实验设计信息
  • 基因:我们检测过的基因信息。

协调我们的操作

  • 避免“簿记”错误，例如，我们将数据的一部分以不同的方式划分为子集。

  • 使用SummarizedExperiment包和数据表示。

    • 二维结构，所以子集用[，]
    • 使用函数访问组件:分析()，rowData ()，rowRanges ()，colData ()等。

库(SummarizedExperiment)
  • 确保顺序样品的列中样本的顺序匹配计数矩阵的顺序基因的行顺序匹配计数矩阵。
  • 创建一个SummarizedExperiment来协调我们的数据操作。
samples <- samples[colnames(counts)，] genes <- genes[rownames(counts)，] se <- summarizeexperiment (assays = list(counts)， rowRanges = genes, colData = samples) se ## class: rangedsummarizeexperiment ## dim: 63677 8 ## metadata(0): ## assays(1): counts ## rownames(63677): ENSG00000000003 ENSG00000000005…ENSG00000273492 ## ENSG00000273493 ## rowData names(2): gene_name gene_biotype ## colnames(8): SRR1039508 SRR1039509…SRR1039520 SRR1039521 ## colData names(3): cell dex avgLength

3.分析与可视化

3．1微分表达式分析

完形

  • 执行一个t.test ()对于计数矩阵的每一行，询问是否泰爱泰党样品的平均计数不同于untrt样本。
  • 许多微妙的统计问题

的DESeq2pacakge

  • 实现高效，“正确”，稳健的算法执行RNA-seq差异表达分析的中等规模的实验。
库(DESeq2)
  • 指定我们的实验设计，考虑到细微的统计问题进行分析，并得到结果的总结。这一步的细节超出了本次研讨会的范围。
dds <- DESeqDataSet(se， ~ cell + dex)拟合<- DESeq(dds) ##估计大小因子##估计分散##基因方面的离散度估计##均值-离散度关系##最终离散度估计##拟合模型和测试去统计值<-结果(拟合)去统计值## log2倍变化(MLE): dex trt vs untrt ## Wald测试p-值:敏捷泰爱泰党vs untrt # # DataFrame 63677行6列# # baseMean log2FoldChange lfcSE # # <数字> <数字> <数字> # # ENSG00000000003 708.602169691234 -0.38125388742934 0.100654430181804 # # ENSG00000000005 0 NA NA # # ENSG00000000419 # ENSG00000000457 237.163036796015 0.0379205923946151 0.14344471633862 520.297900552084 0.206812715390398 0.112218674568195 # # # ENSG00000000460 57.9326331250967 -0.0881676962628265 0.287141995236272  ## ... ... ... ...# # ENSG00000273489 0.275899382507797 1.48372584344306 3.51394515550546 # # ENSG00000273490 0 NA NA NA NA # # # # ENSG00000273491 0 ENSG00000273492 # ENSG00000273493 0.106141666408122 -0.521381077922898 3.53139001322807 0.105978355992386 -0.463691271907546 3.52308373749196 # # #统计pvalue padj # # <数字> <数字> <数字> # # ENSG00000000003 -3.78775069056286 0.000152017272514002 0.00128292609656079 # # ENSG00000000005 NA NA NA # # ENSG00000000419 1.84294384322566 - 0.06533721006625810.196469601297369 ## ensg00000000457 0.26435684326705 0.791504962999781 0.91141814384918 ## ensg00000000460 -0.307052600196215 0.75880333554496 0.895006448013164 ## ... ... ... ...ensg00000273492 -0.131615171950935 0.895288684444562 na ## ensg00000273493 -0.147641884914972 0.88262539793309 na
  • 将结果添加到我们的SummarizedExperiment，这样我们也能以一种协调的方式来操作它们。
rowData(se) <- cbind(rowData(se)， destats)

3.2的热图

  • 使用订单()而且头()以确定表达差异最大的前20个基因的行指数(基于调整后的p值)。
  • 子集SummarizedExperiment只包含这些行。
  • 显示的分析()我们子集的热图
top20idx <- head(order(rowData(se)$padj)， 20) top20 <- se[top20idx，]热图(分析(top20))

更新行标签并添加有关治疗组的信息。

  • 提取前20个矩阵。
  • 用相应的基因名更新矩阵的行名。
  • 创建颜色向量，每个样本一个，颜色由地塞米松处理确定。
$gene_name trtcolor <- hcl. m <- assay(top20) rownames(m) <- rowData(top20)colors(2， "Pastel 1")[colData(top20)$dex]热图(m, ColSideColors = trtcolor)

3．3火山的阴谋

  • y轴表示“统计学意义”，x轴表示“生物学意义”。
  • 使用图()来创建点
  • 使用text ()给两个最重要的基因加上标签。
plot(-log10(pvalue) ~ log2FoldChange, rowData(se)) label <- with(rowData(se)， ifelse(-log10(pvalue) > 130, gene_name， "")) text(-log10(pvalue) ~ log2FoldChange, rowData(se)， label = label, pos = 4, srt=-15)

3.4顶级表格和整齐的数据

目标

  • 提供20个最具差异表达基因的简明总结。

dplyr以及“整洁”的数据

  • 一种方便的方式显示和操作严格的表格数据。
    • “长形式”表，其中每行表示一个观察结果，每列表示在观察结果上测量的属性。
  • 宠物猫:一个data.frame具有更好的显示属性
  • %＞%,例如,Mtcars %>%计数(cyl):接收tibble(或data.frame)的管道资源描述并将其用作少数函数的参数，如count ()在右边。
  • ' Tidy '函数通常返回一个宠物猫()，因此可以被连接在一起。
库(dplyr)

  • 下面的步骤:

    • as_tibble ():创建一个tibble fromrowData (se)
    • select ():选择具体列。
    • 安排():将所有行从小到大排列padj
    • 头():过滤到前20行
rowData(se) %>% as_tibble(rownames = "ensembl_id") %>% select(ensembl_id, gene_name, baseMean, log2FoldChange, padj) %>% arrange(padj) %>% head(n = 20) ## # A tibble: 20 x 5 ## ensembl_id gene_name baseMean log2FoldChange padj ##      ## 1 ENSG00000152583 SPARCL1 997。4.57 4.00e-132 ## 2 ENSG00000165995 CACNB2 495。3.29 7.06e-131 ## 3 ENSG00000120129 DUSP1 3409。2.20e-126 ## 4 ENSG00000101347 SAMHD1 12703。3.77 4.32e-126 ## 5 ENSG00000189221 MAOA 2342。3.35 3.96e-120 ## 6 ENSG00000211445 GPX3 12286。3.73 1.39e-108 ## 7 ENSG00000157214 STEAP2 3009。1.98 1.48e-103 ## 8 ENSG00000162614 NEXN 5393。2.04 2.98e-100 ## 9 ENSG00000125148 MT2A 3656。2.21 5.81e- 94 ## 10 ENSG00000154734 ADAMTS1 30315。 2.35 5.87e- 88 ## 11 ENSG00000139132 FGD4 1223. 2.23 1.24e- 83 ## 12 ENSG00000162493 PDPN 1100. 1.89 1.32e- 83 ## 13 ENSG00000134243 SORT1 5511. 2.20 2.01e- 82 ## 14 ENSG00000179094 PER1 777. 3.19 2.73e- 81 ## 15 ENSG00000162692 VCAM1 508. -3.69 6.78e- 81 ## 16 ENSG00000163884 KLF15 561. 4.46 2.51e- 77 ## 17 ENSG00000178695 KCTD12 2650. -2.53 7.07e- 76 ## 18 ENSG00000198624 CCDC69 2057. 2.92 3.58e- 70 ## 19 ENSG00000107562 CXCL12 25136. -1.91 4.54e- 70 ## 20 ENSG00000148848 ADAM12 1365. -1.81 6.14e- 70

4总结

4．1我们学到了什么

包

数据结构

  • 表示和协调“复杂”数据。
  • 已经普遍存在于基地R,例如,data.frame ()，矩阵()．

  • 农庄表示基因组范围

    • 访问器的功能seqnames ()，start ()等核心组件
    • ＄或mcols()美元的注释

  • SummarizedExperiment用于行和列注释的分析数据的协调操作。

    • [，]以协调的方式将分析和注释子集化。
    • 分析()，rowRanges ()，rowData ()，colData ()访问组件。

分析

  • 成熟的软件包，如DESeq2提供优秀的小插图、定义良好的分析步骤、与其他工作流步骤的集成以及非常健壮的支持。
  • 新兴领域通常由实现分析中不太完整或特定步骤的几个包表示。

4.2下一个步骤

单细胞rna序列:一个惊人的资源:进行单细胞分析与Bioconductor，包括食物而且嘘包。

  • 质量控制
  • 归一化
  • 特征选择
  • 降维
  • 聚类
  • 标记基因检测
  • 单元格类型注释(单)(这个包有一个伟大的小插图!)
  • 轨迹分析(命运）
  • 基因集富集
  • 等等!

其他突出的分析领域(查看biocViews）

  • 微阵列-表观基因组，经典表达，拷贝数
  • 带注释的变体
  • 基因集富集
  • 流式细胞术
  • 蛋白质组学

参与!

5确认

本报告中报告的研究得到了NCI和国家卫生研究院NHGRI的支持，资助编号为U24CA232979、U41HG004059和U24CA180996。内容仅为作者的责任，并不一定代表美国国立卫生研究院的官方观点。

这项工作的一部分得到了陈-扎克伯格倡议DAF的支持，这是硅谷社区基金会的一个顾问基金。

## R version 3.6.1 Patched (2019-12-01 r77489) ##平台:x86_64-apple-darwin17.7.0(64位)##运行在macOS High Sierra 10.13.6 ## ##矩阵产品:default ## BLAS: /Users/ma38727/bin/R-3-6-branch/lib/libRblas。/Users/ma38727/bin/R-3-6-branch/lib/libRlapack。dylib # # # #语言环境:# # [1]en_US.UTF-8 / en_US.UTF-8 en_US.UTF-8 / C / en_US.UTF-8 / en_US。UTF-8 ## ##附加的基本包:##[1]并行stats4统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:[9] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0 ## [11] IRanges_2.20.1 S4Vectors_0.24.0 ## [13] BiocGenerics_0.32.0 BiocFileCache_1.10.2 ## [15] dbplyr_1.4.2 BiocStyle_2.14.0 ## ##通过命名空间加载(并且没有附加):## [1] bitops_1.0-6 bit64_0.9-7 RColorBrewer_1.1-2 ## [4] httr_1.4.1 tools_backports_1.1.5 ## [7] utf8_1.1.4 R6_2.4.1 rpart_1 .1-15 ## [10] Hmisc_4.3-0 DBI_1.0.0 lazyeval_0.2.2 ## [13] colorspace_1.4-1 nnet_7.3-12 tidyselect_0.2.5 ## [13] gridExtra_2.3 bit_1.1-14 curl_4.2 ## [19] compiler_3.6.1 cli_1.1.0 htmlTable_1.13.2 ## [22] bookdown_0.16 scales_1.1.0 checkmate_1.9.4 ## [25] genefilter_1.68.0 rappdirs_0.3.1 string_1 .4.0 ## [28] digest_0.6.23 Rsamtools_2.2.1 foreign_0.8-72 ## [31][46] Rcpp_1.0.3 munsell_0.5.0 fansi_0.4.0 ## [52] zlibbioc_1.32.0 grid_3.6.1 blob_1.2.0 ## [55] crayon_1.3.4 lattice_0.20-38 Biostrings_2.54.0 ## [37] rlang_0.4.2 rstudioapi_0.10 RSQLite_2.1.2 ## [40] acepack_1.4.1 RCurl_1.95-4.12 magrittr_1. 1.5 ## [43] GenomeInfoDbData_1.2.2 Formula_1.2-3 Matrix_1.2-18 ## [46] Rcpp_1.0.3 munsell_0.5.0 fansi_0.4.0 ## [49] lifecycle_0.1.0 stringi_1.4.3 yaml_2.2.0 ## [55] crayon_1.3.4 lattice_0.20-38 Biostrings_2.54.0 ## [58] splines_3.6.1 annotate_1.64.0locfit_1. 1.5-9.1 ## [61] zeallot_0.1.0 knitr_1.26 pillar_1.4.2 ## [64] geneplotter_1.64.0 codetools_0.2-16 XML_3.98-1.20 ## [67] glue_1.3.1 evaluate_0.14 latticeExtra_0.6-28 ## [70] data.table_1.12.6 BiocManager_1.30.10 vctrs_0.2.0 ## [73] gtable_0.3.0 purrr_0.3.3 assertthat_0.2.1 ## [76] ggplot2_3.2.1 xfun_0.11 xtable_1.8-4 ## [79] survival_1 .1-7 tibble_2.1.3 genome alignments_1.22.1 ## [82] AnnotationDbi_1.48.0 memoise_1.1.0 cluster_2.1.0