ChIP-seq 数据比对:技术实践指南
参考基因组的准备
进行基因组比对前,首先需要获取目标基因组的序列数据(FASTA格式)。R语言中的BSgenome包提供了多种物种的完整基因组序列。对于常见的研究物种,如小鼠的mm10版本,可以直接加载对应的BSgenome数据包。
library(BSgenome.Mmusculus.UCSC.mm10)
BSgenome.Mmusculus.UCSC.mm10考虑到分析效率和实际研究需求,通常只选择主要染色体进行后续分析。需要排除随机染色体和未定位的scaffolds。以下代码提取了小鼠的主要染色体序列并创建为DNAStringSet对象:
primaryChroms <- paste0("chr", c(1:19, "X", "Y", "M"))
primaryChrSeq <- lapply(primaryChroms, function(x) BSgenome.Mmusculus.UCSC.mm10[[x]])
names(primaryChrSeq) <- primaryChroms
primaryChrSeqSet <- DNAStringSet(primaryChrSeq)
primaryChrSeqSet随后使用writeXStringSet函数将序列导出为FASTA文件,供比对工具使用:
writeXStringSet(primaryChrSeqSet, "reference_mm10_chroms.fa")基因组索引的构建
大多数高效的序列比对算法都基于FM-index索引结构,因此在使用比对工具前需要先对参考基因组建立索引。本教程采用subread算法的实现版本subjunc进行比对操作,该功能封装在Rsubread包中。构建索引通过buildindex函数完成,仅需要指定索引名称和参考基因组FASTA文件路径:
library(Rsubread)
buildindex("mm10_reference", "reference_mm10_chroms.fa", memory = 8000,
indexSplit = TRUE)注意:索引构建过程会消耗大量内存,默认配置为8GB。如果计算资源有限,可能需要调整参数或使用内存更小的物种基因组。
序列比对操作
方法一:Rsubread
Rsubread包提供了align函数实现DNA序列与参考基因组的比对。该函数支持直接输出BAM格式文件,并允许用户自定义比对参数。
myAlignment <- align("mm10_reference", "filtered_ENCFF001NQP.fastq.gz", output_format = "BAM",
output_file = "Sample_Myc_1.bam", type = "dna", phredOffset = 64, nthreads = 4)方法二:Rbowtie2
Bowtie2是序列比对领域最广泛使用的工具之一,性能稳定且支持多种比对模式。通过Rbowtie2包可以在R环境中调用Bowtie2。首先需要构建Bowtie2索引:
library(Rbowtie2)
bowtie2_build(references = "reference_mm10_chroms.fa", bt2Index = file.path("reference_mm10_chroms"))由于Bowtie2需要输入未压缩的FASTQ文件,在比对前需要先解压原始数据:
library(R.utils)
gunzip("filtered_ENCFF001NQP.fastq.gz", remove = FALSE)
bowtie2(bt2Index = "reference_mm10_chroms", samOutput = "ENCFF001NQP.sam",
seq1 = "filtered_ENCFF001NQP.fastq")Bowtie2默认输出SAM格式,如需转换为BAM格式,可使用RSamtools包的asBam函数:
alignmentBam <- asBam("ENCFF001NQP.sam")使用Rbowtie2时需要特别注意磁盘空间管理。该包不支持流式处理压缩输入,需要确保有足够空间存储中间文件(SAM和未压缩的FASTQ)。完成比对后应及时清理临时文件:
unlink("ENCFF001NQP.sam")BAM文件排序与索引
完成比对生成的BAM文件通常需要按染色体坐标排序并建立索引,以便后续分析工具快速访问。Rsamtools包提供了sortBam和indexBam两个函数完成这些操作:
library(Rsamtools)
sortBam("Sample_Myc_1.bam", "sorted_Sample_Myc_rep1")
indexBam("sorted_Sample_Myc_rep1.bam")排序和索引后的BAM文件可直接用于IGV等基因组浏览器进行可视化展示,也可继续进行Peak Calling等下游分析流程。