乐趣区

关于程序员:ChIPseq-分析数据比对3

jiezi
  • 读取 = reads（二者含意雷同，下文不做辨别）

1. ChIPseq reads 比对

在评估读取品质和咱们利用的任何读取过滤之后，咱们将心愿将咱们的读取与基因组对齐，以便辨认任何基因组地位显示比对读取高于背景的富集。

因为 ChIPseq 读数将与咱们的参考基因组间断比对，咱们能够应用咱们在之前中看到的基因组比对器。生成的 BAM 文件将蕴含用于进一步剖析的对齐序列读取。

2. 参考基因组生成

首先，咱们须要以 FASTA 格局检索感兴趣的基因组的序列信息。咱们能够应用 BSgenome 库来检索残缺的序列信息。对于小鼠 mm10 基因组，咱们加载包 BSgenome.Mmusculus.UCSC.mm10。

library(BSgenome.Mmusculus.UCSC.mm10)
BSgenome.Mmusculus.UCSC.mm10

咱们将仅应用次要染色体进行剖析，因而咱们可能会排除随机和未搁置的重叠群。在这里，咱们循环遍历次要染色体，并依据检索到的序列创立一个 DNAStringSet 对象。

mainChromosomes <- paste0("chr", c(1:19, "X", "Y", "M"))
mainChrSeq <- lapply(mainChromosomes, function(x) BSgenome.Mmusculus.UCSC.mm10[[x]])
names(mainChrSeq) <- mainChromosomes
mainChrSeqSet <- DNAStringSet(mainChrSeq)
mainChrSeqSet

当初咱们有了一个 DNAStringSet 对象，咱们能够应用 writeXStringSet 来创立咱们的 FASTA 序列文件来比对。

writeXStringSet(mainChrSeqSet, "BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa")

3. 索引创立

咱们将应用 subread 背地的 subjunc 算法进行对齐。因而，咱们能够应用 Rsubread 包。在咱们尝试比对咱们的 FASTQ 文件之前，咱们须要首先应用 buildindex() 函数从咱们的参考基因组构建一个索引。

buildindex() 函数仅采纳咱们所需的索引名称和要从中构建索引的 FASTA 文件的参数。

library(Rsubread)
buildindex("mm10_mainchrs", "BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa", memory = 8000,
    indexSplit = TRUE)
  • 请记住：建设索引会占用大量内存，默认状况下设置为 8GB。这对于您的笔记本电脑或台式机来说可能太大了。

4. 比对

4.1. Rsubread

咱们能够应用 Rsubread 包将 FASTQ 格局的原始序列数据与 mm10 基因组序列的新 FASTA 文件进行比对。具体来说，咱们将应用 align 函数，因为它利用了 subread 基因组比对算法。

myMapped <- align("mm10_mainchrs", "filtered_ENCFF001NQP.fastq.gz", output_format = "BAM",
    output_file = "Myc_Mel_1.bam", type = "dna", phredOffset = 64, nthreads = 4)

4.2. Rbowtie2

Bowtie 家族是最驰名的对齐算法之一。咱们能够应用 Rbowtie2 包拜访 Bowtie2。QuasR 包容许拜访原始的 Bowtie 对准器，但它有点慢并且须要内存。

library(Rbowtie2)

与 Rsubread 一样，Rbowtie2 包要求咱们首先创立一个要对齐的索引。咱们能够应用 bowtie2_build() 函数来实现此操作，指定咱们的 FASTA 文件和所需的索引名称。

bowtie2_build(references = "BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa", bt2Index = file.path("BSgenome.Mmusculus.UCSC.mm10.mainChrs"))

而后咱们能够应用 bowtie2() 函数对齐咱们的 FASTQ 数据，指定咱们新创建的索引、SAM 输入的所需名称和未压缩的 FASTQ。

咱们须要先解压缩咱们的 FASTQ。这里咱们应用 remove is FALSE 设置来放弃原始压缩的 FASTQ。

library(R.utils)
gunzip("filtered_ENCFF001NQP.fastq.gz", remove = FALSE)

bowtie2(bt2Index = "BSgenome.Mmusculus.UCSC.mm10.mainChrs", samOutput = "ENCFF001NQP.sam",
    seq1 = "filtered_ENCFF001NQP.fastq")

因为 Rbowtie2 还输入 SAM 文件，咱们须要将其转换为 BAM 文件。咱们能够应用 RSamtools 的 asBam() 函数来做到这一点。

bowtieBam <- asBam("ENCFF001NQP.sam")

应用 Rbowtie2 时的一个重要思考因素是其未压缩文件的输出和输入。在命令行上，咱们能够将输出流式传输到 Rbowtie2，但在 R 中这不是一个选项。咱们须要确保删除任何创立的临时文件（SAM 和 / 或未压缩的 FASTQ）以防止填满咱们的硬盘。咱们能够应用 unlink() 函数删除 R 中的文件。

unlink("ENCFF001NQP.sam")

4.3. 排序

和以前一样，咱们别离应用 Rsamtools 包 sortBam() 和 indexBam() 函数对文件进行排序和索引。生成的排序和索引 BAM 文件当初能够用于内部程序，例如 IGV，也能够用于 R 中的进一步上游剖析。

library(Rsamtools)
sortBam("Myc_Mel_1.bam", "SR_Myc_Mel_rep1")
indexBam("SR_Myc_Mel_rep1.bam")

广告

本文由 mdnice 多平台公布

退出移动版