Hello,
I'd like to identify the peak summit for each peak in a list of ChIPseq
peaks. I've been using IRanges/Genomic Ranges & Rsamtools to read in the BAM
file with mapped reads, seperated reads according to '+' and '-'
strand-mappings ( & dropped unmapped reads).
I have a predefined list of peak regions that I'm interested in, so I only
take the reads that map to these regions.
x=BAM_data[BAM_data %in% Peak_regions]
x_views=views(coverage(x), Peak_regions)
which gives me an Rle looking like this
start end width
[1] 3114398 3114591 194 [0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 ...]
[2] 3198876 3199070 195 [0 0 0 0 0 0 0 0 0 0 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 0 0 0 0 0 0 ...]
Does anybody have some advice on what is the best way to access the coverage
data in this object and find the maximum value/summit?
Many thanks!!
Daniela
[[alternative HTML version deleted]]
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
