Making TSS plots for some genes for a particular histone mark

Question

[Varun Gupta]

Hi Everyone,

I have 2 bam files which I produced by mapping the reads with bowtie using hg19. This is a CHIP seq data for one of the histone marks, where one file is input(control) and other file is the chip sample. I want to show for certain genes which I am interested in, read plots around the Transcription Start Sites(TSS). Let's say I have 6 genes of interest, so in one TSS plot(for that histone mark) I want to see how histone marks are represented for different genes.

How can I do it?

Thanks

Answer 1

[James Ashmore]

I've had great success following the BEDTools protocol for plotting occupancy. Starts on page 18.

[1769mkc]

The pdf is not there can the link is dead. Can you share the pdf again?

[Arup Ghosh]

I guess you can find that protocol in the following URL.

https://github.com/arq5x/lab_website/blob/master/pdf/bedtools.protocols.pdf

The same can be achieved using deepTools.

Ref: https://deeptools.readthedocs.io/en/latest/content/example_gallery.html?highlight=tss#tata-box-enrichments-around-the-tss-of-mouse-genes

[1769mkc]

Found the pdf I was looking that in bedtools site couldn't find it.

Answer 2

[dariober]

See if this Software for drawing peak shapes or this Visualizing average read coverage across features with varying length thread helps.

Answer 3

[Fidel]

DeepTools will do exactly what you need quickly and nicely. You can use our Galaxy instance or install deepTools using python pip command (see our github page for that).

First you need to transform you bam file to a bigwig file (using the tool called bamCoverage from deepTools). The following code will probably work:

bamCoverage -b file.bam -o coverage.bw
# compute the coverage for a region 3000 bp downstream and upstream of the TSS
computeMatrix reference-point -R my_genes.bed -S coverage.bw \
  -o gene_cov --afterRegionStartLength 3000 \
  --beforeRegionStartLength 3000 --referencePoint TSS
heatmapper -m gene.cov -o heatmap_and_profile.pdf
profiler -m gene.cov -o profile.pdf

As a recommendation, I suggest to use a log2 ratio of ChIP-seq / input instead of the coverage as this can be biased.

[Varun Gupta]

HI Fidel,

I was trying to use Deeptools, but don't you think for the chip sample, I should use bamCompare and normalize and get a normalized bigwig file

Something like this

bamCompare --bamfile1 chip1.bam --bamfile2 Input.bam --ratio log2 --scaleFactorsMethod SES --binSize 50 -o chip_vs_input.bw

And then should I use computeMatrix command? I think we should always normalize the chip with input(as you say log2 normalization)

Also what does my_genes.bed file consists of?

I think my bed files should have the TSS site of the gense I am interested in : something like this:

chr1    34565    34566    gene_b    100    +
chr17    43300230    43300231    gene_a    100    -

Will the heatmap show both these genes for the particular histone mark I am interested? I might increase the list of genes in my bed file

Let me know

Thanks for the help

[Fidel]

Yes, use bamCompare. Be careful with the scale factor. The SES method is good for sharp peaks but not for broad marks, so choose based on the mark that you have.

To get the TSS correctly the only requirement is that the bed file contains the strand information, you don't need to specify the TSS. The usual bed for genes like from UCSC is what is needed. The heatmap will show the two genes but I think it works best when you have many more genes.

[Varun Gupta]

Hi Fidel,

Thanks for the explanation. The default scale method is SES and it is good for sharp peaks as you say but, for peaks with H3K27ac , those are usually broader, what scale factor should I use?

About the bed file: What information should I include in it? The one I am using above, is it wrong??

[Fidel]

I recommend you to take a look at the bamCompare manual to understand the scaling methods. The bed file that you are using will also work.

[Varun Gupta]

Hi Fidel,

I am trying to view a heatmap for around 70 of my genes at TSS. Is there a way I can label those genes on the heatmap. All it labels is genes. Also is there a way to label them horizontally and not vertically? That would be a great help. Thanks for the wonderful tool.

Answer 4

[Sander]

Also NGS-Plot is a nice tool!

https://github.com/shenlab-sinai/ngsplot

Answer 5

[Ming Tommy Tang]

I want to add Genomation and my post here