search

nboley / idr

IDR
GNU General Public License v2.0
168 stars 46 forks source link

Using --rank option no output #13

Closed nataliesilmon closed 8 years ago

nataliesilmon commented 8 years ago

When --rank is specified (p.value, signal.value or q.value), no peaks are called. If it is not specified, signal.value is used but for MACS2 called peaks p.value needs to be used.

nboley commented 8 years ago

It seems to work on the test data, e.g.

idr --samples peak1 peak2 --verbose --rank q.value

returns output. What is the exact command that isnt working?

nataliesilmon commented 8 years ago

Hi,

Perhaps I am doing something wrong, sorry - I tried the test files and it does work. I'm using the command below (and have tried --rank p.value and --rank signal.value) and get no peaks, but if I do not specify --rank I get peaks ranked by signal.value (I can tell this from the output file because columns 8 and 9 are set to -1).

_idr --samples H3K4me3_rep1_callpeak_p1e-2_peaks.narrowPeak H3K4me3_rep2_callpeak_p1e-2_peaks.narrowPeak -o H3K4me3_IDR_rep1vsrep2qvalue --idr-threshold 0.05 --rank q.value --plot

Here is the start of my dataset (peaks called in MACS2):

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1a 98 . 2.45547 9.83350 8.70161 225

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1b 174 . 3.16719 17.41792 15.97919 480

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1c 106 . 2.57598 10.66525 9.50756 895

ABPA02001268 251 793 H3K4me3_rep1_callpeak_p1e-2_peak_2 46 . 1.67656 4.65375 3.66030 328

ABPA02001273 332 620 H3K4me3_rep1_callpeak_p1e-2_peak_3 98 . 2.49152 9.80494 8.67419 164

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4a 111 . 2.01017 11.19162 10.01817 524

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4b 117 . 2.06141 11.76029 10.56994 854

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5a 202 . 2.46462 20.20859 18.61673 182

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5b 137 . 2.19984 13.73618 12.47353 504

ABPA02001287 102 993 H3K4me3_rep1_callpeak_p1e-2_peak_6a 103 . 2.13005 10.33780 9.19190 184

On 10 December 2015 at 13:04, Nathan Boley notifications@github.com wrote:

It seems to work on the test data, e.g.

idr --samples peak1 peak2 --verbose --rank q.value

returns output. What is the exact command that isnt working?

— Reply to this email directly or view it on GitHub https://github.com/nboley/idr/issues/13#issuecomment-163704489.

Natalie Clare Silmon de Monerri

email: nsilmondemonerri@gmail.com

nboley commented 8 years ago

Can you please send me a self contained test case?

On Thu, Dec 10, 2015 at 11:33 AM, Natalie Silmon de Monerri < notifications@github.com> wrote:

Hi,

Perhaps I am doing something wrong, sorry - I tried the test files and it does work. I'm using the command below (and have tried --rank p.value and --rank signal.value) and get no peaks, but if I do not specify --rank I get peaks ranked by signal.value (I can tell this from the output file because columns 8 and 9 are set to -1).

_idr --samples H3K4me3_rep1_callpeak_p1e-2_peaks.narrowPeak H3K4me3_rep2_callpeak_p1e-2_peaks.narrowPeak -o H3K4me3_IDR_rep1vsrep2qvalue --idr-threshold 0.05 --rank q.value --plot

Here is the start of my dataset (peaks called in MACS2):

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1a 98 . 2.45547 9.83350 8.70161 225

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1b 174 . 3.16719 17.41792 15.97919 480

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1c 106 . 2.57598 10.66525 9.50756 895

ABPA02001268 251 793 H3K4me3_rep1_callpeak_p1e-2_peak_2 46 . 1.67656 4.65375 3.66030 328

ABPA02001273 332 620 H3K4me3_rep1_callpeak_p1e-2_peak_3 98 . 2.49152 9.80494 8.67419 164

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4a 111 . 2.01017 11.19162 10.01817 524

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4b 117 . 2.06141 11.76029 10.56994 854

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5a 202 . 2.46462 20.20859 18.61673 182

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5b 137 . 2.19984 13.73618 12.47353 504

ABPA02001287 102 993 H3K4me3_rep1_callpeak_p1e-2_peak_6a 103 . 2.13005 10.33780 9.19190 184

On 10 December 2015 at 13:04, Nathan Boley notifications@github.com wrote:

It seems to work on the test data, e.g.

idr --samples peak1 peak2 --verbose --rank q.value

returns output. What is the exact command that isnt working?

— Reply to this email directly or view it on GitHub https://github.com/nboley/idr/issues/13#issuecomment-163704489.

Natalie Clare Silmon de Monerri

email: nsilmondemonerri@gmail.com

— Reply to this email directly or view it on GitHub https://github.com/nboley/idr/issues/13#issuecomment-163727270.

nataliesilmon commented 8 years ago

Hi, Please find script and files attached - I think this is what you meant by a self contained test? If not let me know and I'll send anything else. Thank you!

On 10 December 2015 at 14:33, Natalie Silmon de Monerri < nsilmondemonerri@gmail.com> wrote:

Hi,

Perhaps I am doing something wrong, sorry - I tried the test files and it does work. I'm using the command below (and have tried --rank p.value and --rank signal.value) and get no peaks, but if I do not specify --rank I get peaks ranked by signal.value (I can tell this from the output file because columns 8 and 9 are set to -1).

_idr --samples H3K4me3_rep1_callpeak_p1e-2_peaks.narrowPeak H3K4me3_rep2_callpeak_p1e-2_peaks.narrowPeak -o H3K4me3_IDR_rep1vsrep2qvalue --idr-threshold 0.05 --rank q.value --plot

Here is the start of my dataset (peaks called in MACS2):

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1a 98 . 2.45547 9.83350 8.70161 225

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1b 174 . 3.16719 17.41792 15.97919 480

ABPA02001267 1279 2300 H3K4me3_rep1_callpeak_p1e-2_peak_1c 106 . 2.57598 10.66525 9.50756 895

ABPA02001268 251 793 H3K4me3_rep1_callpeak_p1e-2_peak_2 46 . 1.67656 4.65375 3.66030 328

ABPA02001273 332 620 H3K4me3_rep1_callpeak_p1e-2_peak_3 98 . 2.49152 9.80494 8.67419 164

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4a 111 . 2.01017 11.19162 10.01817 524

ABPA02001282 72 1337 H3K4me3_rep1_callpeak_p1e-2_peak_4b 117 . 2.06141 11.76029 10.56994 854

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5a 202 . 2.46462 20.20859 18.61673 182

ABPA02001286 16 946 H3K4me3_rep1_callpeak_p1e-2_peak_5b 137 . 2.19984 13.73618 12.47353 504

ABPA02001287 102 993 H3K4me3_rep1_callpeak_p1e-2_peak_6a 103 . 2.13005 10.33780 9.19190 184

On 10 December 2015 at 13:04, Nathan Boley notifications@github.com wrote:

It seems to work on the test data, e.g.

idr --samples peak1 peak2 --verbose --rank q.value

returns output. What is the exact command that isnt working?

— Reply to this email directly or view it on GitHub https://github.com/nboley/idr/issues/13#issuecomment-163704489.

Natalie Clare Silmon de Monerri

email: nsilmondemonerri@gmail.com

Natalie Clare Silmon de Monerri

email: nsilmondemonerri@gmail.com

nboley commented 8 years ago

I dont see an attachment.

nataliesilmon commented 8 years ago

Sorry here: test.zip

nboley commented 8 years ago

It looks like the code is working - it's just that no peaks pass the 5% IDR threshold when using the q value ranking. Try using --soft-idr-threshold instead of --idr-threshold and then you can see all of the peaks and their idr values.

nataliesilmon commented 8 years ago

ok, I'll try that. Thank you!