ChIP-seq 分析:数据比对(3)
- 读取 = 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")
相关文章
- 关于Android大数据收集,埋点统计的详细讲解以及案例代码分析附github代码
- 智能视频分析系统的大数据应用
- python pandas对社保数据进行整理整合
- ChIP-seq 分析:数据比对(3)
- 【视频】Python用LSTM长短期记忆神经网络对不稳定降雨量时间序列进行预测分析|数据分享|附代码数据
- Redis数据倾斜与JD开源hotkey源码分析揭秘
- 数据分享|R语言分析上海空气质量指数数据:kmean聚类、层次聚类、时间序列分析:arima模型、指数平滑法|附代码数据
- R语言GARCH族模型:正态分布、t、GED分布EGARCH、TGARCH的VaR分析股票指数|附代码数据
- R语言中贝叶斯网络(BN)、动态贝叶斯网络、线性模型分析错颌畸形数据|附代码数据
- R语言nlme、nlmer、lme4用(非)线性混合模型non-linear mixed model分析藻类数据实例|附代码数据
- R语言使用马尔可夫链对营销中的渠道归因建模|附代码数据
- Oracle数据库升级或数据迁移方法研究
- 字段MySQL中使用Text字段存储大型文本数据(mysql中的text)
- 轻松拷贝Oracle数据迁移新方案(oracle数据拷贝)
- 恢复灾难:MySQL 数据误删案例分析(mysql数据误删)
- 利用 Redis 挖掘数据潜能: Redis 分析工具(redis分析工具)
- 使用Linux Mhash加密算法保护您的数据(linuxmhash)
- 架起数据云梦:搭建MSSQL云数据库(搭建云数据库mssql)
- 时间MSSQL分析一段时间内的数据查询(mssql 查询 一段)
- MySQL去重技巧如何避免数据中包含重复数据(mysql不包含重复)
- Redis集群助力数据深度分析(redis 集群qdas)
- 比较Oracle数据库两个字段数据对比分析(oracle两个字段内容)
- 利用Redis实现简单分页技术(redis返回分页数据)
- Google Nest Hub将增加空气质量数据 警告用户附近的污染和烟雾
- Repeater控件数据导出Excel(附演示动画)
- Android序列化XML数据