समृद्धि का आकलन

R में Bioconductor के साथ ChIP-seq

Peter Humburg

Statistician, Macquarie University

रीड्स को बढ़ाएँ

डेटा लोड करें:

reads <- readGAlignments(bam)
reads_gr <- granges(reads[[1]])

औसत फ्रैगमेंट लंबाई प्राप्त करें:

frag_length <- fragmentlength(qc_report)["GSM1598218"]

रीड्स बढ़ाएँ और कवरेज निकालें:

reads_ext <- resize(reads_gr, width=frag_length)
cover_ext <- coverage(reads_ext)

पीक्स के लिए कवरेज

जीनोम पर 200 bp के बिन बनाएँ।

bins <- tileGenome(seqinfo(reads), tilewidth=200, 
                   cut.last.tile.in.chrom=TRUE)

पीक से ओवरलैप होते सभी बिन खोजें।

peak_bins_overlap <- findOverlaps(bins, peaks)
peak_bins <- bins[from(peak_bins_overlap), ]

हर पीक बिन से ओवरलैप रीड्स की गिनती करें।

peak_bins$score <- countOverlaps(peak_bins, reads)

बिन्ड कवरेज फ़ंक्शन

count_bins <- function(reads, target, bins){
  # Find all bins overlapping peaks
  overlap <- from(findOverlaps(bins, target))
  target_bins <- bins[overlap, ]

  # Count the number of reads overlapping each peak bin
  target_bins$score <- countOverlaps(target_bins, reads)
  target_bins
}

ब्लैकलिस्टेड क्षेत्रों का कवरेज

peak_bins <- count_bins(reads_ext, peaks, bins)
bl_bins <- count_bins(reads_ext, blacklist.hg19, bins)

बैकग्राउंड कवरेज

पहले से गिने गए सभी बिन हटा दें।

bkg_bins <- subset(bins, !bins %in% peak_bins & !bins %in% bl_bins)

शेष हर बिन से ओवरलैप रीड्स की गिनती करें।

bkg_bins$score <- countOverlaps(bkg_bins, reads_ext)

अभ्यास करते हैं!

R में Bioconductor के साथ ChIP-seq