查看原文
其他

Bioconductor序列数据介绍

2017-04-11 Sonali ,Martin 生信媛


原文作者是:Sonali Arora, Martin Morgan,本文已征得作者授权。


本文为我学习Bioconductor所写的学习笔记。说是学习笔记,差不多等于全文翻译了bioconductor的 sequence Analysis一章。基本上高通量数据分析到了最后都会用到Bioconductor的包,甚至是前期都可能会用到,对于英语好的同学而言,可以直接去官网看它的common work flow,对于哪些英文不错,就是懒得看英文的同学,你要是不嫌弃我的翻译水平太差的话,你就勉强看下去吧。欢迎在留言区指出错误,让我们共同进步。

什么是Bioconductor

官方的简介是:

Bioconductor provides tools for the analysis and comprehension of high-throughput genomic data. Bioconductor uses the R statistical programming language, and is open source and open development. It has two releases each year, , and an active user community. Bioconductor is also available as an (Amazon Machine Image) and a series of images.

简单的说Bioconductor就是以R语言为平台的一个高通量基因组数据的分析工具。

简介

Biconductor为高通量基因组数据提供了分析和理解方法。项目提供了许多R包用于对大量的数据进行严格的统计分析,while keeping techological aritiacts in mind. Bioconductor能进行多种类型分析:

  • 测序: RNASeq, ChIPSeq, variants, copy number…

  • 微矩阵: 表达,SNP

  • Domain specific analysis: Flow cytometry, Proteomics…

对于这些分析,通常需要先导入并使用多种序列相关的文件类型,包括fasta, fastq,BAM,gtf,bed,和wig文件等等.Bionconductor支持导入,常见和高级的序列操作,triming,transformation, and alignment(质量评估)

Sequencing Resources

不同阶段,不同数据格式会用到的包

  • IRange, GenomicRanges:基于range的计算,数据操作,常见数据表示。Biostring: DNA和氨基酸序列表示,比对和模式匹配和大规模生物学序列的数据操作。ShortRead处理FASTQ文件

  • Rsamtools, GenomicAlignments用于比对read(BAM文件)的I/O和数据操作 rtracklayer用于导入和导出多种数据格式(BED,WIG,bigWig,GTF,GFF)和UCSC genome browser tracks操作

  • BSgenomes: 用于获取和操作规划的全基因组。GenomicFeatures常见基因组之间序列特征注释,biomaRt 获取Biomart数据库

  • SRAdb用于从Sequnce Read Archive(SRA)查找和获取数据。

Bioconductor包通过biocViews进行管理,有些词条位于Sequencing,其他则在其他条目和代表性的包下,包括

  • , e.g., , , , , and .

  • , e.g.,, , , .

  • and other variants, e.g., , , .

  • e.g., , , .

  • and metagenome sequencing, e.g., , , .

Ranges Infrastructure

许多其他的Bioconductor包高度依赖于IRanges/GenomicRanges所提供的低层数据结构

library(GenomicRanges) GRanges(seqnames=Rle(c('chr1', 'chr2', 'chr3'), c(3, 3, 4)),      IRanges(1:10, width=5), strand='-',      score=101:110, GC = runif(10))

GenomicRanges使得我们可以将染色体坐标的一个范围和一段序列名(如chromosome)以及正负链关联起来。这对描述数据(aligned reads坐标,called ChIP peaks 或copy number variants)和注释(gene models, Roadmap Epigenomics regularoty elements, dbSNP的已知临床相关变异)

GRanes对象,用来储存基因组位置和关联注释的向量。
Genomic ranges基本都是通过导入数据 (e.g., viaGenomicAlignments::readGAlignments())或注释(e.g., via GenomicFeatures::select() or rtracklayer::import() of BED, WIG, GTF, and other common file formats).进行创建

help(package="GenomicRanges") vignette(package="GenomicRanges")

GRanges常用操作包括findOverlaps, nearest等

FASTA文件中的DNA/氨基酸序列

Biostrings类用于表示DNA或氨基酸序列

library(Biostrings) d <- DNAString("TTGAAAA-CTC-N") length(d)

使用AnnotationHub在Ensembl上下载Homo sapiens cDNA序列,文件名为‘Homo_sapiens.GRCh38.cdna.all.fa’

library(AnnotationHub) ah <- AnnotationHub() ah2 <- query(ah, c("fasta","homo sapiens","Ensembl")) fa <- ah2[["AH18522]] fa

使用Rsamtools打开fasta文件读取里面的序列和宽度记录

library(Rsamtools) idx <- scanFaIndex(fa) idx long <- idx[width(idx)>82000] getSeq(fa,param=long)

BSgenome提供全基因组序列,

library(BSgenome.Hsapiens.UCSC.hg19) chr14_range <- GRanges("chr14",IRanges(1,seqlengths(HSapiens)["chr14"])) chr12_dna <- getSeq(Hsapiens,chr14_range) letterFrequency(chr14_dna,‘GC",as.prob=TRUE)

FASTQ文件中的Reads

ShortRead用于fastq文件.BiocParallel用于加速任务进程

##1. attach ShortRead and BiocParallel library(ShortRead) library(BiocParallel) ## 2. create a vectior of file paths fls <- dir("~/fastq",pattern="*fastq",full=TRUE) ## 3. collect statistics stats0 <- qa(fls) ## 4. generate and browser the report if (interactive())    browseURL(report(stats))

BAM中的Aligned Reads

GenomicAlignments用于输入比对到参考基因组的reads
下一个案例中,我们会读入一个BAM文件,and specifically read in reads / supporting an apparent exon splice junction/ spanning position 19653773 of chromosome 14.

内有8个BAM文件。我们只是用第一个BAM文件。我们会载入软件包和数据包,为目标区域构建一个GRanges对象,然后使用summarizeJunctions()寻找目标区域的read

# 1. 载入软件包 library(GenomicRanges) library(GenomicAlignments) # 2. 载入样本数据 library('RNAseqDta.HNRNPC.bam.chr14') bf <- BamFile(RNAseqData.HNRNPC.bam.chr14_BAMFILES[[1),asMates=TRUE) # 3. 定义目标区域 roi <- GRanges("chr14", IRanges(19653773,width=1)) # 4. alignments, junctons, overlapping our roi paln <- readAlignmentsList(bf) j <- summarizeJunctions(paln,with.revmap=TRUE) j_overlap <- j[j %over% roi] # 5.supporting reads paln[j_overlap$revmap[[1]]]

Called Variants from VCF files

VCF(Variant Call Files)描述了SNP和其他变异。该文件包括元信息行,含有列名的header line和众多的数据行,每一行都包括基因组上位置信息和每个位置上可选基因型信息。
数据通过VariantAnnotation的readVcf()读入

library(VariantAnnotation) fl <- system.file("extdata","chr22.vcf.gz",package=“VariantAnnotation") vcf <- readVcf(fl,"hg19")

Genome Annotations from BED, WIG, GTF etc files

能够导入和导出许多常见的文件类型,BED,WIG,GTF,还能查询和导航UCSC genome Browser
rtracklayer包含测序所用BED文件

library(rtracklayer) test_path <- system.file("test",package = "rtracklayer") test_bed <- file.path(test_path,"test.bed") test <- import(test_bed,format="bed") test

同样包括多种基因组注释文件(BED,GTF,BigWig),使用rtracklayer import()。更多有关AnnotationHub的介绍需要看 and 。

我是hoptop,一名不太合格的翻译工。

推荐阅读:

R语言的数据管理

生信人修炼手册

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存