Getting peak heights from TF chIP-seq data (wig file)

Question

[arsala521]

Hello everyone,

I have TF ChIP seq data from NCBI GEO in wig format. I converted wig to bedgraph and then used MACS peak caller to get bed narrowpeak files.I further uploaded file on genome browser to get graphical map of peak regions. This figure is more of a binary (yes/no) indicating just presence of peaks at particular regions.

Is there anyway I can get figure showing peak heights, using the data I have. I don't have raw reads. I only have data in wig, bedgraph and narrowpeak files. I was wondering if there is some way I can get peak heights using signal values from wig/bedgraph or using score values in the fifth field of bed narrowpeak?

Thank you

Answer 1

[Alex Reynolds]

For WIG files:

wig2bed < signal.wig | bedmap --echo --max-element - > answer.bed

Bedgraph files have to be converted to BED:

awk -v FS="\t" -v OFS="\t" '{ print $1, $2, $3, ".", $4 }' in.bedgraph | sort-bed - | bedmap --echo --max-element - > answer.bed

If you have compressed files, use gunzip -c to extract them, e.g:

gunzip -c signal.wig.gz | wig2bed - | bedmap --echo --max-element - > answer.bed

[arsala521]

Thanks for the answer. Would you please clarify a bit more. I can't exactly follow what "bedmap --echo --max-element" is doing. I did run this on my wig file. Just pasting top few lines of my wig and output file:

[D_Melanogaster]$ wig2bed < R1.wig | bedmap --echo --max-element - > answer.bed
[D_Melanogaster]$ head -n5 R1.wig
chrX 5 40 8.01596910471638
chrX 105 140 9.84274069026994
chrX 205 240 10.3630144423714
chrX 305 340 9.76566981010343
chrX 405 440 8.91672795759724
[D_Melanogaster]$ head -n5 answer.bed
chr2L   192     227     id-539832       15.701525|chr2L 192     227     id-539832       15.701525
chr2L   292     327     id-539833       18.030737|chr2L 292     327     id-539833       18.030737
chr2L   392     427     id-539834       16.980151|chr2L 392     427     id-539834       16.980151
chr2L   492     527     id-539835       13.950963|chr2L 492     527     id-539835       13.950963
chr2L   592     627     id-539836       11.384652|chr2L 592     627     id-539836       11.384652

What is the difference between signal values shown in wig and output file? Also to further clarify my question, I have already called peaks from wig file using MACS. I already have peak regions in bed narrowpeak file. Following are top few lines of my bed narrowpeak file:

[D_Melanogaster]$ head -n5 R1_peak.bed
chr2L   193     5728    GoatAb_R1_peak.bed_narrowPeak1  323     .       0       0       0       2585
chr2L   5928    8328    GoatAb_R1_peak.bed_narrowPeak2  404     .       0       0       0       350
chr2L   10128   10728   GoatAb_R1_peak.bed_narrowPeak3  80      .       0       0       0       250
chr2L   11128   12228   GoatAb_R1_peak.bed_narrowPeak4  135     .       0       0       0       450

I need a figure showing height of the peaks. Can the score value in the 5th column of narrowpeak file be used for that?

Thank you

[seidel]

The 5th column is simply a re-scaled q-value from MACS, not really a peak height, but it does tend to correlate. However, if you've run MACS, then don't you have the .xls file which has a column for pileup - which is the height at the peak summit? In addition, you would also have the summits file, which is the base at the peak summit. What will your figure look like?

[seidel]

I want to have figure showing peaks with varying peak heights

I see - well isn't that the wig file? And by "figure" I think you mean to say that you want a track. The wig file would be the signal trace of read heights along the genome (aka read coverage), and the BED file shows the regions of signal that have been classified into a peak. Maybe you want to convert the WIG file to a bigWig (binary version of wig)? Is this correct, or am I misunderstanding you? In this figure top upper track is a bigWig of IP signal, and the lower track is a BED file of peak calls.

[arsala521]

Thank you so much for the very helping reply. Yes I want this upper track. Right now I only have the lower one (with bed file). How can I get this type of figure (track) showing both peaks and signals together? Should I upload both bigwig and bed file on the browser?

[seidel]

Yes, you would load the wig or bigwig file to the browser. Not sure what your set up is (there are various ways of doing this). The wig file is meant to display continuous data along the genome (data that wiggles).

[arsala521]

I want two tracks simultaneously like you showed, BED only demarcating peak positions and wig showing peak heights too, I am usIng Jbrowse genome browser. I think I need to upload both bed and wig files, right? Thank you again.

[seidel]

Yes, upload both files. You should be able to view multiple tracks simultaneously.

[arsala521]

Thank you so much

[arsala521]

Thank you. Yes, I used macs. I converted wig to bedgraph and used macs3 bdgpeakcall function to get bed narrowpeak file:

macs3 bdgpeakcall -i R1.bedgraph --outdir . -o R1_peak.bed

No, I don't get any .xls file which has height at the peak summit. That info would be useful. How can I get that? I also don't get any summit fiie. For figure, right now when I upload my bed file on browser, it demarcates the peak region. This figure is more of a binary (yes/no) indicating just presence of peaks at particular regions. I want a figure with quantitative information showing how strong is the signal for every peak. I want to have figure showing peaks with varying peak heights.