r - 使用 .bam 文件映射 CpG 坐标

Question

我找到了一种在鸡基因组上绘制 CpG 位置的方法，如下：

library(BSgenome.Ggallus.UCSC.galGal4)
library(Biostrings)

#Create the chromossome set

chr_all=paste("chr", c(1:28, 32, "Z", "W"), sep="")

#CpG Mapping on Genome
        CpGpos<-lapply(chr_all,function(x) { 
          RES<-matchPattern("CG",Ggallus[[x]])
          RES<-cbind(rep(x,length(RES)),RES@ranges@start)
            return(RES) 
        })

        CpGpos<-do.call("rbind",CpGpos)
        CpGstartend<-cbind(CpGpos,as.numeric(CpGpos[,2])+1)
        CpGstartend<-as.data.frame(CpGstartend,stringsAsFactors=F) 
write.table(CpGstartend, "CPGstartendWGChicken.txt", sep="\t", col.names=F, row.names=T, quote=F)

结果：

[fibio@milou2 Stacks]$ cat CPGstartendWGChicken.txt |head

    1       chr1    38      39
    2       chr1    84      85
    3       chr1    146     147

...

现在我想对我的 MEDIPseq 结果做同样的事情。我已经有 .bam 格式的对齐序列 (BSgenome.Ggallus.UCSC.galGal4) 可供使用。

超级有趣的是每个 CpG 的覆盖率值列

预期结果：

               Chr   Start    End  Depth
               chr1    38      39    10
               chr1    84      85    12
               chr1    146     147   20

我知道名为 MEDiPS 的 Bioconductor R 包可以做到这一点，但它使用 windows 大小。但是，我需要每个 CpG 的信息。定义 2bp 窗口大小时，脚本被终止。