1. Mapped reads

現(xiàn)在我們有了 BAM 文件的索引，我們可以使用 idxstatsBam() 函數(shù)檢索和繪制映射讀取的數(shù)量。

mappedReads <- idxstatsBam("SR_Myc_Mel_rep1.bam")
TotalMapped <- sum(mappedReads[, "mapped"])
ggplot(mappedReads, aes(x = seqnames, y = mapped)) + geom_bar(stat = "identity") +
    coord_flip()

TotalMapped

2. bigWig 創(chuàng)建

我們還可以從我們排序的、索引的 BAM 文件中創(chuàng)建一個 bigWig，以允許我們快速查看 IGV 中的數(shù)據(jù)。

首先，我們使用 coverage() 函數(shù)創(chuàng)建一個包含我們的覆蓋率分?jǐn)?shù)的 RLElist 對象。

forBigWig <- coverage("SR_Myc_Mel_rep1.bam")
forBigWig

我們現(xiàn)在可以使用 rtracklayer 包的 export.bw() 函數(shù)將 RLElist 對象導(dǎo)出為 bigWig。

library(rtracklayer)
export.bw(forBigWig, con = "SR_Myc_Mel_rep1.bw")

我們可能希望標(biāo)準(zhǔn)化我們的覆蓋范圍，以便我們能夠比較樣本之間的富集。

我們可以使用 coverage() 中的權(quán)重參數(shù)將我們的讀取縮放到映射讀取數(shù)乘以一百萬（每百萬讀取數(shù)）。

forBigWig <- coverage("SR_Myc_Mel_rep1.bam", weight = (10^6)/TotalMapped)
forBigWig
export.bw(forBigWig, con = "SR_Myc_Mel_rep1_weighted.bw")

SR_Myc_Mel_rep1_weighted.bw

本文由mdnice多平臺發(fā)布