Bioconductor序列数据介绍
原文作者是: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,一名不太合格的翻译工。
推荐阅读: