How to visualize selected few ChIPseq peaks

Question

[Mike]

Hello everyone,

I have some selected most significant peaks from MACS2 output bed file, I want to visualize in any Genome Browse (IGV/UCSC etc..). How to visualize only few peaks?

Thanks,

Answer 1

[Devon Ryan]

1. Make a small BED file containing only those peaks.
2. Index that in IGV
3. Load the BED file

That's basically it. You might want to load bigWig tracks at the same time so you can see the underlying signal.

[Mike]

Thanks Devon, really appreciated for your smart response.

Answer 2

[dariober]

Somewhat of a shameless plug, but hopefully useful... I wrote (I'm writing) ASCIIGenome, a text based genome browser designed to work from terminal and optionally in non-interactive mode. In the OP's case it might be handy to loop through the peak regions and output a screenshot for each of them.

Just working from command line, a script like this could do:

while read line
do
  reg=`echo $line | awk '{print $1 ":" $2 "-" $3}'`
  ASCIIGenome -ni -r $reg -x 'zo 3 && ylim 0 na && save .png' my_peaks.narrowPeak coverage.tdf coverage2.tdf
done < my_peaks.narrowPeak

Explanation:

• my_peaks.narrowPeak is a file of selected peaks in bed or gff format (narrowPeak is effectively bed).

• For each line, i.e. peak, in this file prepare a target region to pass to ASCIIGenome in the format chrom:from-to (awk step).

• Run ASCIIGenome in non-interactive mode (-ni), load the peak file itself and some coverage files (tdf format here, can be bigwig or else), go to the target region (-r $reg).

• Then zoom out 3x, reset the y axis limits from 0 to local max, save the screenshot in png format (-x 'zo 3 && ylim 0 na && save .png').

Each output png should look like this:

Caveats: ASCIIGenome is restarted from scratch for each input region, this takes a few seconds to start the JVM itself plus some seconds to load the files. It's still ok-ish but non very smart. And of course the text output is not as pretty as IGV's.

If interested, feel free to send comments, issues etc. See also this post here on Biostars ASCIIGenome: Text Only Genome Viewer!

[Mike]

Thanks dariober, what are coverage.tdf , coverage2.tdf in your above script??

[dariober]

Hi- The tdf format has been devised by the IGV developer(s) to show quantitative data along the genome (e.g. RNA-Seq or ChIP-Seq coverage), see also the recommended formats section from igv. In my script coverage.tdf and coverage2.tdf are just dummy examples which in reality could be two chip samples or chip and input etc.

I like tdf format because you can very quickly show large intervals like entire chromosomes. However, bigWig are also very fast, it's your choice really. Of course you might want to show just the bed intervals regardless of the coverage profile, but usually I prefer to see how the peaks look like to check how credible they are.

Answer 3

[steve]

As mentioned, the same can be accomplished in the UCSC Genome Browser:

1. Navigate here: https://genome.ucsc.edu/cgi-bin/hgCustom

2. Select your reference genome, paste in your BED regions, and click 'Submit'

3. Select the track, and view it in Genome Browser

4. Navigate to your region of interest. You can paste the coordinates you want to see directly into the white box at the top (the one that says "enter position, gene symbol, ... etc.") and hit 'Go' to jump to a specific peak.

One advantage of UCSC is that you are able to share your sessions with collaborators a little easier than with IGV.

[Mike]

Thanks steve, according to your above step, I used my selected peaks files:

chr3    84775414        84775903        R1_peak_2335
chr3    84949410        84949871        R1_peak_2337
chr15   61955243        61955645        R1_peak_1439
chr9    102686656       102687105       R1_peak_3861
chr10   24554476        24554797        R1_peak_303
chr3    84827963        84828360        R1_peak_2336

But in step 3rd (manage custom tracks), its shows 6 items, but only chr3 , why not chr15 , ch9, ch10 ??

[steve]

It looks like it is using the first "chr3" it encounters as the default first position. But if you select that custom track and hit 'Go' for the Genome Browser, you can then paste in the other BED coordinates in the box at the top (see the red arrow in the pic for step 4) to jump to the other positions. So just paste chr10 24554476 24554797 into that box and hit 'Go' and you should see the peak at that region. Of course if the regions you're interested in are close together, you can just keep zooming out until they appear on the same screen.

Answer 4

[Charles Plessy]

You can visualise your data and select peaks of interest in the ZENBU browser. After registering an account, upload:

• The BED file containing MACS2 peaks,
• The BAM file(s) containing your ChIP-Seq alignments.

Then create one annotation track for the peaks and one numerical signal track for the alignments.

After this, you can browse any peak of interest (locating it for instance by its genomic coordinates), and see the shape of the raw data. For a rough impression of how the data can be visualised, you can see the example view "ZENBU - Severin et al. Supp. Fig. 4. ZENBU on-demand dynamic ChipSeq peak calling - SH3BP4". (Not using MACS2 peaks...)

Answer 5

[Persistent LABS]

+1 to Devon.

Also If you are looking for any free platform for carrying out end-to-end ChipSeq analysis along with data visualization, you can try SanGeniX. Its free. It does Peak calling and gives several visualization such as Region-wise distribution of peak (upstream/downstream/genic/intergenic), Chromosome-wise peak distribution, Tag density plot (plot with peaks plotted using tag densites in control and treatment), Fold enrichment plot (distribution of peaks based on fold-enrichment values given by MACS), GO terms distribution corresponding to the peaks in specific regions of the genome, Peak visualization using JBrowse and Motif analysis.

Thanks

Priyabrata

Persistent LABS