Computing P-Value Or Qvalue Or Both From Peak File

Question

[Dataminer]

Hi!

I have a peak file and also its treatment and control BAM file. The peak file has peaks and I have computed the fold enrichment for these peaks by taking ration of tags in treatment to control.

Now, I want to compute p-value for these peaks, how can I do it?

Any scripts available or method that I should follow?

Thank you

[KCC]

Just to be clear. Are you trying to implement your own simple peak caller?

[Dataminer]

No, I have broad peak output file from MACS2 and I want to know how can I get p-value for these peaks.

Answer 1

[KCC]

MACS2 has a README file which has most of the formats described. Here is the link, https://github.com/taoliu/MACS/blob/master/README

NAMEbroadpeaks.bed is in BED12 format which contains both the broad region and narrow peaks. The 5th column is 100*-log10pvalue, to be more compatible to UCSC standard. Tht 7th is the start of the first narrow peak in the region, and the 8th column is the end. The 9th column should be RGB color key, however, we keep 0 here to use the default color, so change it if you want. The 10th column tells how many blocks including the starting 1bp and ending 1bp of broad regions. The 11th column shows the length of each blocks, and 12th for the starts of each blocks.

Note that the 5th column has the information you want. I would probably leave them in the -log10 (pvalue) format as they are are often quite small values. However, if y is the value in column 5, then you want (-y/100) for the log10(pvalue) and 10^(-y/100) for the pvalue. If you want I can write a script to do this, but essentially, you just want to pull out the 5th column and do the computation I described.

[Dataminer]

It worked, problem solved thank you.

[KCC]

Just a small edit. I forgot that the logarithm was base 10. So, to get the p-value, it should be 10^(-y/100)

Answer 2

[Puriney]

Maybe this post? Peak Finding Of Chip-Seq (Or Clip-Seq) And Poisson Distribution, Am I Right?