Introduction to Bioconductor for Sequence Data
在院士组轮转的时候，由于没人安排我去做什么，我也不懂怎么和别人交流，于是那些日子就在翻译Biocondutor的教程中度过。本文为我学习Bioconductor所写的第一篇学习笔记。说是学习笔记，差不多等于全文翻译了bioconductor的 sequence Analysis一章。基本上高通量数据分析到了最后都会用到Bioconductor的包，甚至是前期都可能会用到，
什么是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, 1296 software packages, and an active user community. Bioconductor is also available as an AMI (Amazon Machine Image) and a series of Docker images.
简单的说Bioconductor就是以R语言为平台的一个高通量基因组数据的分析工具。那么下面就算正文环节
简介
Biconductor为高通量基因组数据提供了分析和理解方法。我们拥有大量包可对大量数据进行严格的统计分析，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，其他则在其他条目和代表性的包下，包括
ChIPSeq, e.g.,DiffBind, csaw, ChIPseeker, ChIPQC.
SNPs and other variants, e.g., VariantAnnotation, VariantFiltering, h5vc.
CopyNumberVariation e.g., DNAcopy, crlmm, fastseg.
Microbiome and metagenome sequencing, e.g., metagenomeSeq, phyloseq, DirichletMultinomial.
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.
RNAseqData.HNRNPC.bam.chr14_BAMFILES内有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
rtracklayer能够导入和导出许多常见的文件类型，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
AnnotationHub 同样包括多种基因组注释文件（BED,GTF,BigWig)，使用rtracklayer import()。更多有关AnnotationHub的介绍需要看 Annotation workflow and AnnotationHub HOW TO vignette。
顺便放一下自己的知识星球，如果你觉得我对你有帮助的话。
