序列分析生物导体

马丁•摩根
2015年2月2日

序列分析工作流程

  1. 实验设计
  2. 湿式实验室序列准备
  3. 测序(Bentley et al.， 2008，doi: 10.1038 / nature07517）
  4. 对齐
  5. 减少
  6. 分析
  7. 理解

数据移动

Alt测序生态系统

序列数据表示

DNA /氨基酸序列:FASTA文件

输入和操作:Biostrings

>NM_078863_up_2000_chr2L_16764737_f chr2L:16764737-16766736 gttggggcccaccagtgccaaaatacacaagaagaacagcatctt gacactaaaatgaaaattgctttgcgtcaatgactcaaaacgaaaatg…atgggtatcaagttgccccgtataaaaggcaagtttaccggttgcacggt >NM_001201794_up_2000_chr2L_8382455_f chr2L:8382455-8384454 ttatttatgtaggcgcccgttcccgcagccagagcactcagaattccggg cgtgtagcgcaacgaccatctacaaggcaatattttgatcgcttgttaggttagg…

整个基因组:2位而且.fa格式:rtracklayer，Rsamtools；BSgenome

读取:FASTQ文件

输入和操作:ShortReadreadFastq ()，FastqStreamer ()，FastqSampler ()

@err127302.1703 hwi - eas350_041:1:1460:19184#0/1 cctgagtgaagctgatcttcttagagagagagagatcttgatcgtcgaggaggagatgctgaccttgacct + hhghhghhhhhhdgg < gdgge@gdggd > ce ?=896=: @err127302.1704 hwi - eas350_041:1:1460:16861#0/1 gcggtatgctggaaggtgctcgaatggagagcgccagcgccccggcgctgagccgccccccc>ed4 > eee > de8eeede8b ? eb <@3; ba79 ?， 881b ?@73;########################

对齐读取:BAM文件(例如ERR127306_chr14.bam)

输入和操作:“低级”Rsamtools，scanBam ()，BamFile ()；“高级”GenomicAlignments

被称为变量:VCF文件

输入和操作:VariantAnnotationreadVcf ()，readInfo ()，readGeno ()选择性地与ScanVcfParam ()．

基因组注释:BED, WIG, GTF等文件

输入:rtracklayer进口()

序列数据R/Bioconductor

Bioconductor的作用

Alt测序生态系统

序列

BiostringsDNA或氨基酸序列类

方法

相关的包

例子

suppressPackageStartupMessages({library(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

范围

范围是:

许多常见的生物学问题都是基于范围的

的GenomicRanges包定义基本类和方法

业务范围

Alt范围代数

范围

Intra-range方法

Inter-range方法

Between-range方法

例子

suppressPackageStartupMessages({library(GenomicRanges)}) gr <- GRanges("A"， IRanges(c(10,20,22)， width=5)， "+") shift(gr, 1) # 1-based坐标!
## seqnames ranges strand ##    ## [1] A [11,15] + ## [2] A [21,25] + ## [3] A [23,27] + ## ------- ## seqinfo:来自未指定基因组的1个序列;没有seqlengths
Range (gr) # intra-range
## seqnames ranges strand ##    ## [1] A [10,26] + ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths
减少(gr) # inter-range
## seqnames ranges strand ##    ## [1] A [10,14] + ## [2] A [20,26] + ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths
覆盖(gr)
##长度为1的RleList ## $A ##整数-长度为26的rllist，运行6次##长度:9 5 5 2 3 2 ##值:0 1 0 1 2 1
Setdiff (range(gr)， gr) # '内含子'
## seqnames ranges strand ##    ## [1] A [15,19] + ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths

IRangesList, GRangesList

参考

GenomicAlignments(对齐读取)

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

方法

例子

suppressPackageStartupMessages({library(GenomicRanges) library(GenomicAlignments) library(Rsamtools)}) ## our 'region of interest' roi <- GRanges("chr14"， IRanges(19653773, width=1)) ##样本数据suppressPackageStartupMessages({library('RNAseqData.HNRNPC.bam.chr14')}) bf <- BamFile(RNAseqData.HNRNPC.bam. bam. GRanges)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(称为变种)

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

函数和方法

例子

##输入变量suppressPackageStartupMessages({library(VariantAnnotation)}) fl <- system. txtfile("extdata"， "chr22. hsapiens .gz"， package="VariantAnnotation") vcf <- readVcf(fl， "hg19") seqlevels(vcf) <- "chr22" ##已知基因模型suppressPackageStartupMessages({library(TxDb.Hsapiens.UCSC.hg19. knowngene)}) coding <- locateVariants(rowData(vcf)， TxDb.Hsapiens.UCSC.hg19. hg19. hg19. rowData)。knownGene, codingvariations())头(编码)
GRanges对象有6个范围和9个元数据列:## seqnames range strand | LOCATION LOCSTART LOCEND ## #    |    ## 1 chr22[50301422, 50301476] * |编码939 939 ## 2 chr22[50301476, 50301476] * |编码885 885 ## 3 chr22[50301488, 50301488] * |编码873 873 ## 4 chr22[50301494, 50301494] * |编码867 867 ## 5 chr22[50301584, 50301584] * |编码777 777 ## 6 chr22 [50302962，698] * |编码698 50302962 # # QUERYID TXID CDSID GENEID PRECEDEID # # <整数> <人物> < IntegerList > <人物> < CharacterList > 24 75253 218562 79087 # # 1 # # 2 25 75253 218562 79087 # # 3 26 # # 4 27 75253 218562 79087 75253 218562 79087 5 # # 28 75253 218562 79087 # 6 # 57 75253 218563 79087 # # FOLLOWID # # < CharacterList > # # 1 # # 2 # 3 # 4 # 5 # # 6  ## ------- ## seqinfo: 1从一个未指明的基因组序列;没有seqlengths

相关的包

参考

rtracklayer(基因组注释)

序列分析旅行团

这个非常开放的主题指向一些最突出的Bioconductor包序列分析。利用这个实验的机会来探索下面突出显示的软件包小插图和帮助页面;欧洲杯2021体育彩票许多材料将在后续的实验和讲座中更详细地介绍。

基础知识

特定于领域的分析——探索以下两个或三个包的登录页、小插图和参考手册。

工作序列，对齐，常见的web文件格式，和原始数据;这些包在很大程度上依赖于IRanges/GenomicRanges我们稍后会讲到的基础设施。

可视化

实验室

短读质量评估

fastqc是一个Java程序，通常用于总结fastq文件的质量。它有一个直观的图形用户界面。这里我们将使用命令行版本。

  1. 从内部Rstudio，选择“Tools - > Shell…”，或使用Mac / linux终端或Windows上的PuTTY程序登录到您的亚马逊机器实例。

  2. 在示例fastq文件上运行fastqc，将输出发送到~ / fastqc_report目录中。

    Fastqc fastq/*fastq——threads 8——outdir=fastqc_reports

  3. 在线研究质量报告和结果文档:在“文件”页签中，单击fastqc_reports．点击那里的HTML文件，然后点击“在Web浏览器中查看”。

ShortRead提供类似的功能，但是从内部R．下面显示了R可以处理大数据，并说明了一些基本的方式，其中一个可能与功能的交互提供Bioconductor包中。

# # 1。附加ShortRead和BiocParallel suppressPackageStartupMessages({library(ShortRead) library(BiocParallel)}) ## 2。创建文件路径fls <- dir("~/fastq"， pattern="*fastq"， full=TRUE)
# # 3。收集统计信息stats0 <- qa(fls)
# # 4。生成并浏览报告if (interactive())

检查所有车道的质量报告

数据(qa_all) if (interactive())

比对(和基因组注释)

此数据来自气道Bioconductor注释包;看到装饰图案有关详细信息,

Bioconductor:我们将探索如何在不同类型的标识符之间进行映射，如何导航基因组坐标，以及如何查询BAM文件以进行对齐读取。

  1. 附上包含基因符号信息的“注释”包org.Hs.eg.db基因组坐标(如基因、外显子、cd、转录本)TxDb.Hsapiens.UCSC.hg19.knownGene．安排TxDb包中的“seqlevels”(染色体名称)与BAM文件中的名称相匹配。

  2. 使用适当的org . *包映射从基因符号到Entrez基因id，并相应TxDb。*包检索SPARCL1基因的基因坐标。注意:下面使用单个基因符号，但我们可以使用1、2或a中的所有基因符号矢量化时尚。

  3. 附加GenomicAlignments用于处理对齐读取的包。使用range ()以获得跨越SPARCL1第一个和最后一个外显子的基因组坐标。输入重叠SPARCL1的配对读。

  4. 关于这些排列，你可以很容易地回答哪些问题?例如，有多少读取重叠这个感兴趣的区域?

    # # 1。a“Annotation”包suppressPackageStartupMessages({library(txdb . hspapiens . ucsc .hg19. knowngene) library(org.Hs.eg.db)}) ## 1。b——将TxDb文件中记录的“seqlevels”映射到## BAM文件fl <- "~/igv/genomes/hg19_alias。tab" map <- with(read.delim(fl, header=FALSE, stringsAsFactors=FALSE)， setNames(V1, V2)) seqlevels(TxDb.Hsapiens.UCSC.hg19。knownGene, force=TRUE) <- map ## 2。Symbol -> Entrez ID ->基因坐标sym2eg <- select(org.Hs.eg.db， "SPARCL1"， "ENTREZID"， " Symbol ") exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19. exonby)knownGene， "gene") sparcl1exons <- exByGn[[sym2eg$ENTREZID]] ## 3。Aligned reads suppressPackageStartupMessages({library(GenomicAlignments)}) fl <- "~/bam/SRR1039508_sorted. "bam" sparcl1gene <- range(sparcl1exons) param <- ScanBamParam(which=sparcl1gene) aln <- readGAlignmentPairs(fl, param=param)

  5. 作为另一个练习，我们问我们输入的读取数据中有多少与已知的基因模型兼容。我们必须找到属于我们基因的转录本，然后按转录本分组的外显子

    5. # #。txids <- select(TxDb.Hsapiens.UCSC.hg19)。knownGene, sym2eg$ENTREZID， "TXID"， "GENEID")$TXID exByTx <- exonsBy(TxDb.Hsapiens.UCSC.hg19。knownGene，“tx”)[txids] ## 5;b compatible alignment hits <- findCompatibleOverlaps(query=aln, subject=exByTx) good <- seq_along(aln) %in% queryHits(hits) table(good)
    ##好##假真## 14

  6. 最后，让我们从基因模型转到蛋白质编码序列。(a)提取按转录本分组的CDS区域，只选择我们感兴趣的转录本，(b)从相应的参考基因组中附加并提取编码序列。将编码序列转化为蛋白质。

    rerestoreseqlevels (TxDb.Hsapiens.UCSC.hg19.knownGene)
    # # TxDb对象:# # # Db型:TxDb支持包:# # # # # # GenomicFeatures数据来源:UCSC基因组:# # # # # # hg19生物:智人# # # UCSC的表:knownGene # # #资源URL: http://genome.ucsc.edu/ # # #的基因类型ID: Entrez基因ID # # #完整数据集:是的# # # miRBase构建ID: GRCh37 # # # transcript_nrow: 82960 # # # exon_nrow: 289969 # # # cds_nrow: 237533 # # # Db由:GenomicFeatures包从Bioconductor # # #创建时间:2014-09-26 11:16:12 -0700(2014年9月26日星期五)## #基因组特征创建时的版本:1.17.17 ## RSQLite创建时的版本:0.11.4 ## DBSCHEMAVERSION: 1.0
    ## a. cds坐标，按转录txids分组<- select(TxDb.Hsapiens.UCSC.hg19)。knownGene, sym2eg$ENTREZID， "TXID"， "GENEID")$TXID cdsByTx <- cdsBy(TxDb.Hsapiens.UCSC.hg19。knownGene， "tx")[txids] ## b.来自相关参考基因组的编码序列suppressPackageStartupMessages({library(BSgenome.Hsapiens.UCSC.hg19)}) dna <- extractTranscriptSeqs(BSgenome.Hsapiens.UCSC. ucsc . c)。hg19, cdsByTx)蛋白<-翻译(dna)

研究基因组范围

访问“GenomicRanges HOWTOs”小插图。

browseVignettes(“GenomicRanges”)

阅读第1节，做练习2.2、2.4、2.5、2.8、2.12和2.13。也许可以选择你特别感兴趣的其他主题。

资源

R/Bioconductor

刊物及简报