1. Mapped reads
当初咱们有了 BAM 文件的索引,咱们能够应用 idxstatsBam() 函数检索和绘制映射读取的数量。
mappedReads <- idxstatsBam("SR_Myc_Mel_rep1.bam")
TotalMapped <- sum(mappedReads[, "mapped"])
ggplot(mappedReads, aes(x = seqnames, y = mapped)) + geom_bar(stat = "identity") +
coord_flip()
2. bigWig 创立
咱们还能够从咱们排序的、索引的 BAM 文件中创立一个 bigWig,以容许咱们疾速查看 IGV 中的数据。
首先,咱们应用 coverage() 函数创立一个蕴含咱们的覆盖率分数的 RLElist 对象。
forBigWig <- coverage("SR_Myc_Mel_rep1.bam")
forBigWig
咱们当初能够应用 rtracklayer 包的 export.bw() 函数将 RLElist 对象导出为 bigWig。
library(rtracklayer)
export.bw(forBigWig, con = "SR_Myc_Mel_rep1.bw")
咱们可能心愿标准化咱们的覆盖范围,以便咱们可能比拟样本之间的富集。
咱们能够应用 coverage() 中的权重参数将咱们的读取缩放到映射读取数乘以一百万(每百万读取数)。
forBigWig <- coverage("SR_Myc_Mel_rep1.bam", weight = (10^6)/TotalMapped)
forBigWig
export.bw(forBigWig, con = "SR_Myc_Mel_rep1_weighted.bw")
本文由 mdnice 多平台公布