作者:马丁·摩根(mtmorgan@fredhutch.org）
时间:2015年9月7日
回车间大纲

本文件中的材料要求R3.2版及Bioconductor版本3.1

stopifnot(getRversion() >= '3.2' && getRversion() < '3.3'， BiocInstaller::biocVersion() >= "3.1")

Bioconductor序列分析的“基础设施”

类、方法和包

本节重点介绍类、方法和包，目标是学习如何导航帮助系统和交互式发现工具。

动机

序列分析是专门的

• 大数据需要以内存和时间效率高的方式进行处理
• 针对序列数据的独特特征，开发了特定的算法

额外的注意事项

• 重用现有的、经过测试的代码比重新发明轮子更容易，更不容易出错。
• 当包共享相似的数据结构时，包之间的互操作性更容易。

解决方案:使用定义良好的类表示复杂数据;方法操作类以执行有用的功能。类和方法放在一起并作为包这样我们就可以从其他人的辛勤工作和测试代码中受益。

核心包

VariantAnnotation | v GenomicFeatures | v BSgenome | v rtracklayer | v GenomicAlignments | | v v summarizeexperimental Rsamtools ShortRead | | | | v v v v v v v v GenomicRanges Biostrings | | v v GenomeInfoDb (XVector) | | v v IRanges | v (S4Vectors)

核心类

案例研究:IRanges而且农庄

的IRangesPackage定义了一个用于指定整数范围的重要类，例如:

- IRanges(start=c(10,20,30)， width=5

长度3 ##开始结束宽度## [1]10 14 5 ## [2]20 24 5 ## [3]30 34 5

在范围上可以执行许多有趣的操作，例如，侧面()识别相邻范围

侧面(ir, 3)

长度3 ##开始结束宽度## [1]7 9 3 ## [2]17 19 3 ## [3]27 29 3

的IRanges类是类层次结构的一部分。要看到这个，向R询问类红外类的类定义IRanges类

类(红外)

# #[1]“IRanges”# # attr(“包”)# #[1]“IRanges”

getClass(类(ir))

##类“IRanges”[包“IRanges”]## ##槽位:## ##名称:起始宽度NAMES elementType elementMetadata ##类:整数整数characterORNULL字符DataTableORNULL ## ##名称:元数据##类:列表## ##扩展:##类“范围”，直接##类“IntegerList”，由类“范围”，距离2 ##类“RangesORmissing”，由类“范围”，距离2 ##类“AtomicList”，由类“范围”，距离3 ##类“列表”，由类“范围”，距离4 ##类“矢量”，由类“范围”，距离5 ##类“注解”，由类“范围”，距离6 ## ##已知子类:“NormalIRanges”

请注意,IRanges扩展了范围类。现在试着进入?旁边（?”旁边,<选项卡> "如果不使用_RStudio，其中<选项卡>意思是按TAB键要求TAB补全)。你可以看到有帮助页面侧面操作几个不同的类。选择补全

?”旁边,Ranges-method”

并验证您所处的页面描述的方法与IRanges实例。探索其他基于范围的操作。

的GenomicRangesPackage扩展了范围的概念，以包括与序列分析中范围应用相关的特性，特别是将范围与序列名(例如，染色体)和链相关联的能力。创建一个农庄实例，该实例基于IRanges实例，如下所示

库(GenomicRanges) gr < -农庄(c(“chr1”、“chr1”,“chr2”),红外光谱、链= c ("+", "-", "+")) gr

## seqnames ranges strand ##    ## [1] chr1 [10,14] + ## [2] chr1 [20,24] - ## [3] chr2 [30,34] + ## ------- ## seqinfo:来自未指定基因组的2个序列;没有seqlengths

侧翼序列的概念在生物学上有更微妙的含义。特别是我们可能期望侧翼序列在+链在值域之前，负链在值域之后。验证侧面适用于农庄对象具有此行为。

侧面(gr, 3)

## seqnames ranges strand ##    ## [1] chr1 [7,9] + ## [2] chr1 [25,27] - ## [3] chr2 [27,29] + ## ------- ## seqinfo:来自未指定基因组的2个序列;没有seqlengths

发现哪些类农庄扩展，找到记录的行为的帮助页侧面当应用于农庄对象，并验证帮助页是否记录了我们刚刚观察到的行为。

类(gr)

# #[1]“农庄”# # attr(“包”)# #[1]“GenomicRanges”

getClass(类(gr))

##类“GRanges”[包“GenomicRanges”]## ##插槽:## ##名称:seqnames ranges strand elementMetadata seqinfo ##类:Rle IRanges Rle DataFrame seqinfo ## ##名称:元数据##类:列表## ##扩展:直接# # # #“GenomicRanges”类,类“向量”,由“GenomicRanges”类,距离2 # #类“GenomicRangesORmissing”,由“GenomicRanges”类,距离2 # #类“GenomicRangesORGRangesList”,由“GenomicRanges”类,距离2 # #类“GenomicRangesORGenomicRangesList”,按类“GenomicRanges”,距离2 # #类“RangedDataORGenomicRanges”,按类“GenomicRanges”,距离2 # #类“注释”,由“GenomicRanges”类,距离3

?”旁边,GenomicRanges-method”

注意可用的侧面()方法中定义的方法进行了扩展GenomicRanges包中。

似乎有一些有用的方法可以用于基因组范围的研究;我们可以从命令行中发现其中一些，这表明这些方法应该是当前的search ()路径

方法(class = "农庄”)

= $ $<- %in% ## [5] < <= == > [9] >= BamViews GenomicFiles NROW ## [13] Ops ROWNAMES ScanBamParam ScanBcfParam ## [17] [[<- aggregate anyNA ## [21] appas .]字符为。complex as。data.frame ## [25] as。Env as.integer as。列表。逻辑## [29]as。数字。raw bamWhich<- blocks ## [33] browseGenome c chrom chrom<- ## [37] coerce coerce<- compare countOverlaps ## [41] coverage disjoin disjointBins distance ## [45] distanceToNearest duplestelementmetadata elementMetadata<- ## [49] end end<- eval export ## [53] extractROWS extractUpstreamSeqs findOverlaps侧边## [57]follow gaps getPromoterSeq granges ## [61] head high2low intersect isDisjoint ## [65] length liftOver mapCoords mapFromTranscripts ## [69] mapToTranscripts match mcols mcols<- ## #[73]元数据元数据<- mstack names ## [77] names<- narrow nearest order ## [81] overlapsAny pack parallelSlotNames parallelVectorNames ## [85] pgap pintersect pmapcods pmapFromTranscripts ## [89] pmapToTranscripts precede promoters psetdiff ## [93] punion range ranges ranges<- ## [97] rank reduce reduceByFile reduceByRange ## [101] relist relistToClass rename rep ## [105] rep.int replaceROWS resize restrict ## [109] rev rowRanges<- scanFa scanTabix ## [113] score score<- seqinfo seqinfo<-## [117] seqlevelsInUse seqnames seqnames<- setdiff ## [121] shift shiftApply show showAsCell ## [125] sort split split<- start ## [129] start<- strand strand<-子集## [133]subsetByOverlaps summarizeOverlaps表尾## [137]tapply tile trim union ## [141] unique update updateObject values ## [145] values<- width width<- window ## [149] window<- with xtfrm ## see '?方法来访问帮助和源代码

使用帮助()中列出帮助页GenomicRanges包,小插曲()查看和访问可用的小插图;这些也可以在Rstudio的“帮助”选项卡中找到。

帮助(包=“基因组范围”)小插图(包=“基因组范围”)小插图(包=“基因组范围”，“基因组范围”)

GenomicRanges

的农庄而且GRangesList类

另外:“TxDb”包提供了基因模型的R表示

txdb <- txdb . hsapiens . ucsc .hg19. knowngene . library(txdb . hsapiens . ucsc .hg19. knowngene

外显子()：农庄

外显子(txdb)

## seqnames ranges strand | exon_id ##    |  ## [1] chr1 [11874, 12227] + | 1 ## [2] chr1 [12595,12721] + | 2 ## [3] chr1 [126136,12721] + | 3 ## [4] chr1 [12646,12697] + | 4 ## [5] chr1 [13221, 14409] + | 5 ## ... ... ... ... ... ...## [289965] chrY [27607404, 27607432] - | 277746 ## [289966] chrY [27635919, 27635954] - | 277747 ## [289967] chrY [59358329, 59359508] - | 277748 ## [289968] chrY [59360007, 59360115] - | 277749 ## [289969] chrY [59360501, 59360854] - | 277750 ## ------- # seqinfo: hg19基因组93个序列(1个循环)

exonsBy ()：GRangesList

exonsBy (txdb tx)

## seqnames ranges strand | exon_id exon_name exon_rank ##    |    ## [1] chr1 [11874,12227] + | 1  1 ## [2] chr1 [12613,12721] + | 3  2 ## [3] chr1 [13221, 14409] + | 5  3 ## ## $2 ## # GRanges对象有3个范围和3个元数据列:# # seqnames范围链| exon_id exon_name exon_rank # # [1] chr1(11874、12227)+ | 1 < NA > 1 # # [2] chr1(12595、12721)+ | 2 < NA > 2 # # [3] chr1(13403、14409)+ | 6 < NA > 3 # # # # # # 3美元农庄对象与范围和3元数据列:# # seqnames范围链| exon_id exon_name exon_rank # # [1] chr1(11874、12227)+ | 1 < NA > 1 # # [2] chr1(12646、12697)+ | 4 < NA > 2 # # [3] chr1(13221、14409)+ | 5 < NA > 3  ## ## ...## <82957更多元素> ## ------- ## seqinfo: 93个序列(1个循环)来自hg19基因组

农庄/GRangesList非常有用

• 代表注释-基因、变异、调控元件、拷贝数区域、…
• 代表数据-对齐读取，ChIP峰值，称为变体，…

基因组范围代数

许多生物学上有趣的问题代表了范围上的运算

• 计数对齐reads和已知基因之间的重叠GenomicRanges: summarizeOverlaps ()
• 最接近调控区域的基因GenomicRanges::最近的()[ChIPseeker] []
• 称为与临床表型相关的变异VariantFiltering

农庄代数

• Intra-range方法
  • 独立于同一对象中的其他范围
  • GRanges变体链感知
  • 转变()，狭窄的()，侧面()，发起人()，调整()，限制()，削减()
  • 看到" ? intra-range-methods
• Inter-range方法
  • 取决于同一对象中的其他范围
  • range ()，reduce ()，空白()，分离()
  • 覆盖()(！)
  • 看到" ? inter-range-methods
• Between-range方法
  • 两个(或多个)范围对象的函数
  • findOverlaps ()，countOverlaps ()、……% / %，%在%，% %外；联盟()，相交()，setdiff ()，punion ()，pintersect ()，psetdiff ()

Biostrings(DNA或氨基酸序列)

类

• XString, XStringSet，例如，DNAString(基因组)，DNAStringSet(读取)

方法- - - - - -

• 作弊剪切
• 操作,例如,reverseComplement ()
• 总结,例如,letterFrequency ()
• 匹配,例如,matchPDict ()，matchPWM ()

相关的包

• BSgenome
  • 全基因组表达
  • 模型和定制
• ShortRead
  • FASTQ文件

例子

• 全基因组序列由ENSEMBL、NCBI和其他机构以FASTA文件的形式分发;模式生物全基因组序列被包装成更用户友好的BSgenome包。下面计算跨chr14的GC内容。

  库(BSgenome.Hsapiens.UCSC.hg19) chr14_range = GRanges("chr14"， IRanges(1, seqlength (Hsapiens)["chr14"])) chr14_dna <- getSeq(Hsapiens, chr14_range) letterFrequency(chr14_dna， "GC"， as.prob=TRUE)

  ## g | c ## [1，] 0.336276

GenomicAlignments(对齐读取)

类——像基因组一样的行为

• GAlignments, GAlignmentPairs, GAlignmentsList

方法

• readGAlignments ()，readGAlignmentsList ()
  • 容易限制输入，迭代块
• summarizeOverlaps ()

例子

• 在14号染色体的19653707 + 66M = 19653773位置找到支持上述连接的reads

  library(GenomicRanges) library(GenomicAlignments) library(Rsamtools) ##我们的“感兴趣的区域”roi <- GRanges(“chr14”，IRanges(19653773, width=1)) ##样本数据库('RNAseqData.HNRNPC.bam.chr14') bf <- BamFile(RNAseqData.HNRNPC.bam。chr14_BAMFILES[[1]]， asMates=TRUE) ##对齐，连接，重叠我们的roi paln <- readGAlignmentsList(bf) j <- summarizejoins (paln, with.revmap=TRUE) j_overlap <- j[j %over% roi] ##支持读取paln[j_overlap$revmap[[1]]]

  ##长度为8的GAlignmentsList对象:## [[1]]## GAlignments对象，2对齐和0元数据列:## seqnames strand cigar qwidth开始结束宽度njunc# # [1] chr14 - 66M120N6M 72 19653707 19653898 192 1 ## [2] chr14 + 7m1270n65m72 19652348 19653689 1342 1 ## ## [[2]] ## GAlignments对象，2对齐和0元数据列:## seqnames绞线雪茄qwidth开始结束宽度njunc# # [1] chr14 - 66M120N6M 72 19653707 19653898 192 1 ## [2] chr14 + 72M 72 19653686 19653757 72 0 ## ## [[3]] ## GAlignments对象2对齐和0元数据列:## seqnames绞线雪茄qwidth开始结束宽度njunc# # [1] chr14 + 72M 72 19653675 19653746 72 0 ## [2] chr14 - 65M120N7M 72 19653708 19653899 192 1 ## ##…## <5个更多的元素> ## ------- ## seqinfo:来自未知基因组的93个序列

VariantAnnotation(称为变种)

类——类似基因组范围的行为

• VCF -“宽”
• VRanges -“高”

函数和方法

• I/O和过滤:readVcf ()，readGeno ()，readInfo ()，readGT ()，writeVcf ()，filterVcf ()
• 注释:locateVariants ()(变量重叠范围)，predictCoding ()，summarizeVariants ()
• 单核苷酸多态性:genotypeToSnpMatrix ()，snpSummary ()

例子

• 从VCF文件中读取变体，并对已知的基因模型进行注释

  ##输入变量库(VariantAnnotation) fl <-系统。file("extdata"， "chr22. vcfens .gz"， package="VariantAnnotation") vcf <- readVcf(fl， "hg19") seqlevels(vcf) <- "chr22" ##已知基因模型库(TxDb.Hsapiens.UCSC.hg19. knowngene) coding <- locateVariants(rowRanges(vcf)， TxDb.Hsapiens.UCSC.hg19. hg19. hg19. hg19. vcf)knownGene, codingvariations())头(编码)

  GRanges对象有6个范围和9个元数据列:# # seqnames范围链|位置LOCSTART LOCEND QUERYID TXID # # < Rle > < IRanges > < Rle > | <因素> <整数> <整数> <整数> <人物> # # 1 chr22(50301422、50301422)- |编码939 939 24 75253 # # 2 chr22(50301476、50301476)25 - |编码885 885 75253 # # 3 chr22(50301488、50301488)26 - |编码873 873 75253 # # 4 chr22(50301494、50301494)27 - |编码867 867 75253 # # 5 chr22(50301584、50301584)28 - |编码777 777 75253 # # 6 chr22 (50302962,50302962] - |编码698 698 57 75253 ## CDSID GENEID preferdeid FOLLOWID ##     ## 1 218562 79087 ## 2 218562 79087 ## 3 218562 79087 ## 4 218562 79087 ## 5 218562 79087 ## 6 218563 79087 ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths

相关的包

• ensemblVEP
  • 正向变异集合变异效应预测器
• VariantTools，h5vc
  • 叫变异

参考

• Obenchain, V, Lawrence, M, Carey, V, Gogarten, S, Shannon, P, and Morgan, M. VariantAnnotation:用于探索和注释遗传变异的生物导体包。生物信息学，2014年3月28日首次在线发布doi: 10.1093 /生物信息学/ btu168

rtracklayer(基因组注释)

• 进口(): BED, GTF, WIG, 2bit等
• export (): GRanges to BED, GTF, WIG，…
• 访问UCSC基因组浏览器

SummarizedExperiment

• 将实验数据与样本、特征和实验范围的注释集成
• 矩阵，其中行由基因组范围索引，列由数据帧。

函数和方法

• 访问器:分析()/分析()，rowData ()/rowRanges ()，colData ()，元数据()
• 特别是基于范围的操作subsetByOverlaps ()

标准文件格式的输入和表示

BAM文件的对齐读取-GenomicAlignments

召回:整体工作流程

1. 实验设计
2. 湿实验室制备
3. 高通量测序
4. 对齐
   • 全基因组vs.转录组
5. 总结
6. 统计分析
7. 理解

对齐读取的BAM文件

• 头

  @HD VN:1.0 SO:坐标@SQ SN:chr1 LN:249250621 @SQ SN:chr10 LN:135534747 @SQ SN:chr11 LN:135006516…@SQ号:chrY LN:59373566 @PG ID:TopHat VN:2.0.8b CL:/home/hpages/ TopHat -2.0.8b。Linux_x86_64/tophat——match -inner-dist 150——solexa-quals——max-multihits 5——no- disdant——no-mixed——covere -search——microexon-search——libraries -type fr- un——num-threads 2——output-dir tophat2_out/ERR127306 /home/hpages/ bowtib2 -2.1.0/indexes/hg19 fastq/ERR127306_1。fastq fastq / ERR127306_2.fastq

• 比对

  • ID，旗帜，对中和同伴

    ERR127306.7941162 403 chr14 19653689 3 72M = 19652348 -1413…ERR127306.22648137 145 chr14 19653692 1 72M = 19650044 -3720…

  • 顺序与质量

    ．.．GAATTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCC  *'%%%%%#&&%''#'&%%%)&&%%$%%'%%'&*****$))$)'')'%)))&)%%%%$'%%%%&"))'')%)) ...TTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAG  '**)****)*'*&*********('&)****&***(**')))())%)))&)))*')&***********)****

  • 标签

    ．.．AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:72 YT:Z:UU NH:i:2 CC:Z:chr22 CP:i:16189276 HI:i:0…AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:72 YT:Z:UU NH:i:3 CC:Z:= CP:i:19921600 HI:i:0

• 通常情况下，排序(按位置)和索引('。白的文件)

GenomicAlignments

• 使用一个示例BAM文件(fl可以是您自己的BAM文件的路径)

  ##示例BAM数据库(RNAseqData.HNRNPC.bam.chr14)让R知道这是一个BAM文件，而不仅仅是一个字符向量库(Rsamtools) bfl <- BamFile(fl)

• 将数据输入R

  aln <- readGAlignments(bfl) aln . aln

  与800484年# # GAlignments对象对齐和0元数据列:# # seqnames链雪茄qwidth开始结束宽度njunc # # < Rle > < Rle > <人物> <整数> <整数> <整数> <整数> <整数> # # [1]chr14 + 72 72 19069583 19069654 72 0 # # [2] chr14 + 72 72 19363738 19363809 72 0 # # [3] chr14 - 72 72 19363755 19363826 72 0 # # [4] chr14 + 72 72 19369799 19369870 72 0 # # [5] chr14 - 72 72 19369828 19369899 72 0  ## ... ... ... ... ... ... ... ... ...[800480] chr14 - 72M 72 106989780 106989851 72 0 ## [800481] chr14 + 72M 72 106994763 106994834 72 0 ## [800482] chr14 - 72M 72 106994819 106994890 72 0 ## [800483] chr14 + 72M 72 107003080 107003151 72 0 ## [800484] chr14 - 72M 72 107003171 107003242 72 0 ## ------- # seqinfo:来自未知基因组的93个序列

  • readGAlignmentPairs ()/readGAlignmentsList ()如果对端数据
  • 有很多事要做，包括农庄/GRangesList操作

  方法(类=类(aln))

  # # [1] ! = % % < <= ## [ 5 ] == > >= NROW # # [9] ROWNAMES[[< -总# # [13]anyNA追加。字符为。complex ## [17] as.data.frame as。Env as.integer as。列出##[21]为。逻辑。数字。[25] cigar coerce compare countOverlaps ## [29] coverage duplicate elementMetadata elementMetadata<- # [33] end eval export extractROWS ## [37] findCompatibleOverlaps finddoverlaps findSpliceOverlaps granges ## [41] grglist head high2low连接## [45]length mapCoords mapFromAlignments mapToAlignments ## [49] match mcols mcols<-元数据##[53]元数据<- mstack names names<- # [57] narrow njunc overlapsAny parallelSlotNames ## [61] pintersect pmapCoords pmapFromAlignmentspmapToAlignments ## [65] qnarrow qwidth ranges rank ## [69] relist relistToClass rename rep ## [73] rep.int replaceROWS rev rglist ## [77] rname rname<- seqinfo seqinfo<- ## [81] seqlevelsInUse seqnames seqnames<- shiftApply ## [85] show showAsCell sort split ## [89] split<- start strand strand<- ##[93]子集subsetByOverlaps summarizeOverlaps table ## [97] tail tapply unique update ## [101] updateObject values values<- width ## [105] window window<- with xtfrm ## see '?方法来访问帮助和源代码

• 购者自慎: BAM文件过大。通常你会限制输入特定的基因组范围，或迭代通过BAM文件。关键Bioconductor函数(例如,GenomicAlignments: summarizeOverlaps ()为您执行此数据管理步骤。请参阅下一节!

其他格式和包

〇大数据BiocParallel，GenomicFiles

限制

• 只输入必要的数据，例如:ScanBamParam ()
• 哪一个:感兴趣的基因组范围
• 什么: BAM文件的“列”，例如，“seq”，“flag”

迭代

• 读取整个文件，但要分块读取
• 块大小小到可以轻松放入内存，
• 块大小大到可以从中受益R的向量化操作-每次记录10k到1M条
• 例如,BamFile(…，收益率Size=100000)

迭代规划模型

• 收益率一大块数据
• 地图输入数据以方便表示，通常将输入总结为简化形式
  • 例如，已对齐的读取坐标计数感兴趣的重叠区域
  • 例如，对齐读序列到GC内容
• 减少跨映射块

• 使用GenomicFiles: reduceByYield ()

  library(GenomicFiles) yield <- function(bfl){##输入一组序列库(GenomicAlignments) readGAlignments(bfl, param=ScanBamParam(what="seq"))} map <- function(aln){##计数G或C核苷酸每个读取库(Biostrings) gc <- letterFrequency(mcols(aln)$seq， " gc ") ##汇总0,1，…G或C核苷酸表(1 + gc, 73) # max。读取长度:72}减少<- ' + '

• 例子

  library(RNAseqData.HNRNPC.bam.chr14) fls <- RNAseqData.HNRNPC.bam。chr14_BAMFILES bf <- BamFile(fls[1]， yieldSize=100000) gc <- reduceByYield(bf, yield, map, reduce) plot(gc, type="h"， xlab="每次对齐读取的gc内容"，ylab="读取数量")

  

平行的评价

• 核心、计算机、集群、云
• 通常，需要内存管理技术，如限制或迭代-并行进程竞争共享内存

• 许多问题是尴尬的平行- - - - - -拉普兰人()-特别是在生物信息学中，并行评估是跨文件的

• 示例:几个BAM文件中的GC内容

  library(BiocParallel) gc <- bplapply(BamFileList(fls)， reduceByYield, yield, map, reduce) library(ggplot2) df <- stack(as.data.frame(lapply(gc, cumsum)) df$ gc <- 0:72 ggplot(df, aes(x= gc, y=values)) + geom_line(aes(color =ind)) + xlab(“每次读取的gc核苷酸数量”)+ ylab(“读取的数量”)