This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Same data, very different heatmaps

Me and a ex lab member generated two heatmaps using the same Pol II dataset using HOMER and deeptools.

I used deeptools, and the ex lab member used HOMER. The heatmaps look 'very' different, though their patterns are similar. My PI however seems to think that because of this, my heatmap is un-reliable so I'm here hoping that someone could shed some light on the discrepancies.

HOMER heatmap was generated using protocol from this site: http://homer.salk.edu/homer/chipseq/heatmap.html

Command:

annotatePeaks.pl refseqTSS.bed hg19 -gtf reference_data.gtf -size 10000 -hist 50 -ghist -bedGraph bedgraphfiles.bdg

Deeptools heatmap was generated using the following command:

computeMatrix reference-point \
  -S /Users/Carlos/Dropbox/ChIP-Seq/hg19/old_data/Sample_Pol-II-Chip-N20/Tracks/Pol-II-Chip.normalized.bigWig \
  -R /Users/Carlos/Desktop/Projects/GenomicDataPublication_Project/Gencode/genes_proteincoding_gencodev19.bed \
  -a 5000 \
  -b 5000 \
  -out computeMatrix_pol \
  --referencePoint TSS \
  --averageTypeBins max \
  --missingDataAsZero \
  --skipZeros \
  --numberOfProcessors "max" \
  -bs 50

plotHeatmap \
  -m /Users/Carlos/Desktop/Projects/GenomicDataPublication_Project/Gencode/computeMatrix_pol \
  -out pol_heatmap2.pdf \
  --sortUsing mean \
  --averageTypeSummaryPlot mean \
  --colorList white black \
  --heatmapHeight 15 \
  --samplesLabel "Pol II" \
  --legendLocation none \
  --heatmapWidth 3

The images generated:

HOMER:

Home image

DeepTools:

DeepTools image

Any ideas?

chip-seq javatreeviewer homer deeptools
• 720 views • link
updated by Ram Ram • written by Sam Sam

I think you may want to first merge the data from each output in the same image so you can be confident how each gene is affected. Also, just some thoughts, you may want to look at the output values from each pipeline. See why they are different, maybe one is normalized in one way in which the other isn't. Also there are default settings you may be using in one but not the other, so flags that you are using or not using which may affect the results.

• 0 views • link
updated by Ram Ram • written by Jason Jason

1 answer

Maybe a difference in expected inputs and types of data? HOMER appears to use peak bed calls, DeepTools uses BigWig so probably coverage info (e.g. not limiting to only peak call regions). Scaling issue? The problem with the HOMER output is there isn't a scale associated with it as there is with DeepTools.

I tend to trust the DeepTools version myself.

EDIT: and what Jason said above :)

• 0 views • link
updated by Ram Ram • written by Chris Fields Chris Fields

Would the results vary that much between peak files and bigwigs? I figured that might have something to do with it, but the difference is just so big.

Scaling could definitely be a problem, JavaTreeViewer messes with scaling a lot.

I've tried to make sure most default settings are tuned so that they are at least similar in both runs.

I guess the real question is how do I justify to my PI that what he is seeing is 'correct'? In his eyes, the difference representation is too noticeable to justify publishing the figure in a manuscript.

• 0 views • link
updated by Ram Ram • written by Sam Sam

Actually, I made on key assumption that likely isn't correct, namely that everything in both runs is centered around the TSS. Key question: would the Pol II data in the BED file given to HOMER (pol.bed) represent the peak calls (e.g. regions of enrichment), or regions around the TSS? If the former (peak calls) then what you are seeing is a heat map generated only around the peak call regions (so would be enriched for signal for obvious reasons), whereas in the DeepTools output the data are generated around all TSS regions (which may have reduced signal overall).

If this is the case, one way to check how comparable they are is to run HOMER as shown on this page but look for regions around TSS instead of peak calls. Conversely, you could try using the Pol II peak calls from the HOMER run in DeepTools.

• 0 views • link
updated by Ram Ram • written by Chris Fields Chris Fields

You could be on to something. The bedfile given to homer is indeed the peak calls from MACS2 and is not the regions around the TSS. However, according to the ex lab member, what the heatmap is showing is each line represents a gene (~30k genes) and approximately 16.5k genes have Pol II in a -5k/5k window, so would your assumption still be correct based on this?

I realize it's hard to say for sure because you can only make guesses without having the data in front of you, but it'll give me something to look into tomorrow morning.

• 0 views • link
written by Sam Sam

Acc. to the HOMER docs for annotatePeaks.pl the -ghist option should output counts for each gene (so each TSS in the provided GTF file), but this appears to require read coverage counts from file(s) passed via -bedGraph (which aren't in the given command line call example). If those aren't present I'm not sure what is going on, since acc. to the docs you don't need peak calls for generating heat maps around TSS.

Frankly, the most comparable run to me would be something like (based blindly off the HOMER docs):

annotatePeaks.pl tss hg19 -size 10000 -hist 50 -ghist -wig <WIG file(s)>

or

annotatePeaks.pl hg19 -gtf reference_data.gtf -size 10000 -hist 50 -ghist -wig <WIG file(s)>

Where the .wig files are generated from the same BigWig information used for DeepTools, and/or the GTF is comparable to the BED file used for DeepTools.

• 0 views • link
updated by Ram Ram • written by Chris Fields Chris Fields

Sorry yes you're right. Bedgraph files were indeed passed to the original command. That's my fault. The peakfile is actually a TSS file of RefSeq TSS. I've updated the example. I will attempt to run the command with a wig file instead of a bedgraph and see if there's a difference.

• 0 views • link
written by Sam Sam

I'm having trouble getting the heatmaps .cdt files to be sorted by Pol II density. The script that was left over by a former 'actual' bioinformatic person does not work because it relies on a very specific MySQL database that is based on RefSeq TSS annotation (and I am attempting to redo this using gencode). Could you point me in the direction to be able to sort my .cdt files by density? Is there a tool for this? Would I need to script something by myself?

• 0 views • link
written by Sam Sam

Log in to answer this question.

More posts like this