How to compare peak profiles for a subset of ChIPseq genes

Question

[bertb]

Hello,

My goal is to compare the peak profiles between WT and MT strains for a subset of genes in a yeast experiment I have completed. Specifically, I would like to:

• generate a ranked list of the genes with highest occupancy in the WT strain. (Though I know I can do this with the set of differentially bound genes I have generated, I want to find the top 500 (or so) genes, whether they are differentially expressed or not.)
• create a composite peak profile (will something like computeMatrix work for this?) of these genes to compare against my mutant

Any advice would be greatly appreciated!

Thank you, B

Answer 1

[Rory Stark]

Addressing the first question, if you do a differential analysis in DiffBind, you can get a report with all the sites by calling dba.report() and setting th=1. This will include the mean read count (in log2 form) for every peak across each sample group. By default these will be mean normalized read counts (you can get raw reads by setting bNormalized=FALSE).

Here's an example using the sample dataset with all the consensus peaks sorted by mean read density in the tamoxifen Resistant and Responsive samples:

data(tamoxifen_analysis)
consensus  <- dba.report(tamoxifen, th=1)
resistant  <- consensus[order(consensus$Conc_Resistant,  decreasing=TRUE),]
responsive <- consensus[order(consensus$Conc_Responsive, decreasing=TRUE),]

[bertb]

Great - thank you