glsearch36 e-value threshold

zorino commented 6 years ago

The option (-E) to filter the evalue results doesn't seem to work with glsearch36. I'm using the last release v36.3.8d_13Apr16

wrpearson commented 6 years ago

I need a bit more information before I can understand what the problem with glsearch36 is.

(1) The current version (now available at GitHub) is version 36.3.8f Sep, 2017, could you try your search with this version?

(2) sometimes, the statistical estimates fail, which produces E()-values of -inf, is this what you are seeing?

If you could send me your query sequence, command line, and the output through the list of high scores (I do not need to see the alignments), I should be able to help more.

Bill Pearson

Begin forwarded message:

From: Maxime Déraspe notifications@github.com<mailto:notifications@github.com>

Subject: [wrpearson/fasta36] glsearch36 e-value threshold (#2)

Date: December 21, 2017 at 2:51:19 PM EST

To: wrpearson/fasta36 fasta36@noreply.github.com<mailto:fasta36@noreply.github.com>

Cc: Subscribed subscribed@noreply.github.com<mailto:subscribed@noreply.github.com>

Reply-To: wrpearson/fasta36 reply@reply.github.com<mailto:reply@reply.github.com>

The option (-E) to filter the evalue results doesn't seem to work with glsearch36. I'm using the last release v36.3.8d_13Apr16

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/wrpearson/fasta36/issues/2, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AGGGvwKpfAR_WrqkWFSY7aZwlC83Kui6ks5tCra2gaJpZM4RKRoQ.

zorino commented 6 years ago

Dear Dr. Pearson,

1) OK it's because 36.3.8f Sep, 2017 is not part of an official release (https://github.com/wrpearson/fasta36/releases) I've clone the latest github version (at commit cecd77af3832d5b5d2bb45f2d7cf5b58cc7e3e9f)

2) Yes that's exactly what I'm seeing. For instance I filter with -E 1e-40 but I get results with high evalues, 97 for instance. It is also the case with the latest github version.

I'm also wondering how could I speed up a global-local search ? Does filtering have an impact at all ? maybe by quickly dropping highly divergent alignment ?

Thank you.

This is the command and output result

$ glsearch36 -b 12 -E 1e-40 -m 8 U00096.rna ../GCA_000794425.1_AZPAE14550_genomic.fna

b0201|rRNA|16S  JTXW01000042.1  86.38   1543    209     8       1       1542    355     1890    0       439.5
b0201|rRNA|16S  JTXW01000042.1  56.41   1654    629     211     1       1542    3099    4653    1.1e+02 26.3
b0204|rRNA|23S  JTXW01000042.1  85.52   2909    418     23      1       2904    2361    5251    0       712.0
b0205|rRNA|5S   JTXW01000042.1  80.00   120     24      0       1       120     5394    5513    9.8e-259        84.1
b2588|rRNA|5S   JTXW01000042.1  80.00   120     24      0       1       120     5394    5513    6.8e-256        83.7
b2589|rRNA|23S  JTXW01000042.1  85.41   2910    421     25      1       2904    2361    5251    0       682.9
b2591|rRNA|16S  JTXW01000042.1  85.93   1543    216     8       1       1542    355     1890    0       430.7
b2591|rRNA|16S  JTXW01000042.1  56.27   1654    631     211     1       1542    3099    4653    1.2e+02 26.2
b3272|rRNA|5S   JTXW01000042.1  77.50   120     27      0       1       120     5394    5513    1e-200  76.5
b3272|rRNA|5S   JTXW01000042.1  66.97   131     36      22      1       120     530     649     1.1e+03 12.0
b3272|rRNA|5S   JTXW01000042.1  70.00   137     33      27      1       120     3217    3343    1.1e+03 12.0
b3274|rRNA|5S   JTXW01000042.1  80.00   120     24      0       1       120     5394    5513    4e-256  83.8
b3275|rRNA|23S  JTXW01000042.1  85.34   2909    423     23      1       2904    2361    5251    0       686.0
b3278|rRNA|16S  JTXW01000042.1  85.61   1542    221     6       1       1542    355     1890    0       444.6
b3278|rRNA|16S  JTXW01000042.1  56.20   1654    632     211     1       1542    3099    4653    71      26.7
b3756|rRNA|16S  JTXW01000042.1  85.68   1542    220     6       1       1542    355     1890    0       440.7
b3756|rRNA|16S  JTXW01000042.1  56.09   1653    634     209     1       1542    3099    4653    91      26.5
b3758|rRNA|23S  JTXW01000042.1  85.38   2909    422     23      1       2904    2361    5251    0       697.3
b3759|rRNA|5S   JTXW01000042.1  80.00   120     24      0       1       120     5394    5513    1.3e-255        83.7
b3851|rRNA|16S  JTXW01000042.1  85.80   1543    218     8       1       1542    355     1890    0       438.3
b3851|rRNA|16S  JTXW01000042.1  56.30   1653    631     209     1       1542    3099    4653    67      26.8
b3854|rRNA|23S  JTXW01000042.1  85.68   2911    413     26      1       2905    2361    5251    0       696.6
b3855|rRNA|5S   JTXW01000042.1  78.33   120     26      0       1       120     5394    5513    1.2e-222        79.5
b3855|rRNA|5S   JTXW01000042.1  66.37   132     38      19      1       120     3465    3589    1.1e+03 11.6
b3968|rRNA|16S  JTXW01000042.1  85.86   1543    217     8       1       1542    355     1890    0       431.9
b3968|rRNA|16S  JTXW01000042.1  56.23   1653    632     209     1       1542    3099    4653    1.3e+02 26.0
b3970|rRNA|23S  JTXW01000042.1  85.41   2910    421     25      1       2904    2361    5251    0       684.6
b3971|rRNA|5S   JTXW01000042.1  79.17   120     25      0       1       120     5394    5513    2.9e-223        79.6
b3971|rRNA|5S   JTXW01000042.1  58.04   122     47      10      1       120     3465    3578    1.1e+03 11.1
b4007|rRNA|16S  JTXW01000042.1  85.86   1543    217     8       1       1542    355     1890    0       431.8
b4007|rRNA|16S  JTXW01000042.1  56.23   1653    632     209     1       1542    3099    4653    1.3e+02 26.0
b4009|rRNA|23S  JTXW01000042.1  85.48   2910    419     25      1       2904    2361    5251    0       684.0
b4010|rRNA|5S   JTXW01000042.1  79.17   120     25      0       1       120     5394    5513    1.4e-220        79.2
b4010|rRNA|5S   JTXW01000042.1  58.04   122     47      10      1       120     3465    3578    1.1e+03 11.2

wrpearson commented 6 years ago

Could you send me the ‘normal’ output without -m 8?

Thanks, Bill Pearson

On Dec 28, 2017, at 9:42 PM, Maxime Déraspe notifications@github.com<mailto:notifications@github.com> wrote:

Dear Dr. Pearson,

OK it's because 36.3.8f Sep, 2017 is not part of an official release (https://github.com/wrpearson/fasta36/releases) I've clone the latest github version (at commit cecd77ahttps://github.com/wrpearson/fasta36/commit/cecd77af3832d5b5d2bb45f2d7cf5b58cc7e3e9f)
Yes that's exactly what I'm seeing. For instance I filter with -E 1e-40 but I get results with high evalues, 97 for instance. It is also the case with the latest github version.

I'm also wondering how could I speed up a global-local search ? Does filtering have an impact at all ? maybe by quickly dropping highly divergent alignment ?

Thank you.

This is the command and output result

$ glsearch36 -b 12 -E 1e-40 -m 8 U00096.rna ../GCA_000794425.1_AZPAE14550_genomic.fna

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/wrpearson/fasta36/issues/2#issuecomment-354389865, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AGGGv308oPrSdkOvDXwderlmE1gn8dGOks5tFFF-gaJpZM4RKRoQ.

wrpearson commented 6 years ago

I believe I understand what is happening here, which is confusing but not completely unexpected.

It looks like you are getting a large number of hits to the same sequence, JTXW01000042.1.

This very much confuses the counting/ordering process used to limit the display of output.

Basically, glsearch36 is finding a large number of secondary alignments after finding the primary alignment. This would be more clear in the conventional output (-m 0), because the secondary alignments would start with a line starting with "+ " in the description line, and a '>-' to separate the primary (>>) from the secondary (>-) alignments. In the -m 8 output, this information is obscured, but can be inferred from the fact that the subject sequence is the same for all the alignments.

The threshold for secondary alignments is set by the second value in the -E option, e.g. -E "1e-40 -1" would prevent any secondary alignments (which would cause you to miss many homologies), -E "1e-40 1e-40" would show all secondary alignments with E() < 1e-40. The fact that you have an alignment with E() < 1.1e+2 suggests that there is a bug in the thresholding process, as you suggested.

The -b 12 option restricts the number of primary alignments to 12, but does not restrict the number of secondary alignments.

So you have found a bug in the thresholding process for secondary alignments, which I will look into.

Bill Pearson

Begin forwarded message:

From: Maxime Déraspe notifications@github.com<mailto:notifications@github.com>

Subject: Re: [wrpearson/fasta36] glsearch36 e-value threshold (#2)

Date: December 28, 2017 at 9:42:06 PM EST

To: wrpearson/fasta36 fasta36@noreply.github.com<mailto:fasta36@noreply.github.com>

Cc: William Pearson wrp@virginia.edu<mailto:wrp@virginia.edu>, Comment comment@noreply.github.com<mailto:comment@noreply.github.com>

Reply-To: wrpearson/fasta36 reply@reply.github.com<mailto:reply@reply.github.com>

Dear Dr. Pearson,

OK it's because 36.3.8f Sep, 2017 is not part of an official release (https://github.com/wrpearson/fasta36/releases) I've clone the latest github version (at commit cecd77ahttps://github.com/wrpearson/fasta36/commit/cecd77af3832d5b5d2bb45f2d7cf5b58cc7e3e9f)
Yes that's exactly what I'm seeing. For instance I filter with -E 1e-40 but I get results with high evalues, 97 for instance. It is also the case with the latest github version.

I'm also wondering how could I speed up a global-local search ? Does filtering have an impact at all ? maybe by quickly dropping highly divergent alignment ?

Thank you.

This is the command and output result

$ glsearch36 -b 12 -E 1e-40 -m 8 U00096.rna ../GCA_000794425.1_AZPAE14550_genomic.fna

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/wrpearson/fasta36/issues/2#issuecomment-354389865, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AGGGv308oPrSdkOvDXwderlmE1gn8dGOks5tFFF-gaJpZM4RKRoQ.

wrpearson commented 6 years ago

I believe the problem with E-value thresholds with glsearh36 (and ggsearh36) has been corrected.

The problem was occurring with very low E()-value thresholds (< 1e-20) which were not converted to score thresholds properly. I have modified the function used to convert from E()-values to z-scores, so that it should work better for very low thresholds.

The GitHub version has been updated and is now 36.3.8f Dec.2017.

Bill Pearson

Begin forwarded message:

From: Maxime Déraspe notifications@github.com<mailto:notifications@github.com>

Subject: Re: [wrpearson/fasta36] glsearch36 e-value threshold (#2)

Date: December 28, 2017 at 9:42:06 PM EST

To: wrpearson/fasta36 fasta36@noreply.github.com<mailto:fasta36@noreply.github.com>

Cc: William Pearson wrp@virginia.edu<mailto:wrp@virginia.edu>, Comment comment@noreply.github.com<mailto:comment@noreply.github.com>

Reply-To: wrpearson/fasta36 reply@reply.github.com<mailto:reply@reply.github.com>

Dear Dr. Pearson,

OK it's because 36.3.8f Sep, 2017 is not part of an official release (https://github.com/wrpearson/fasta36/releases) I've clone the latest github version (at commit cecd77ahttps://github.com/wrpearson/fasta36/commit/cecd77af3832d5b5d2bb45f2d7cf5b58cc7e3e9f)
Yes that's exactly what I'm seeing. For instance I filter with -E 1e-40 but I get results with high evalues, 97 for instance. It is also the case with the latest github version.

I'm also wondering how could I speed up a global-local search ? Does filtering have an impact at all ? maybe by quickly dropping highly divergent alignment ?

Thank you.

This is the command and output result

$ glsearch36 -b 12 -E 1e-40 -m 8 U00096.rna ../GCA_000794425.1_AZPAE14550_genomic.fna

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/wrpearson/fasta36/issues/2#issuecomment-354389865, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AGGGv308oPrSdkOvDXwderlmE1gn8dGOks5tFFF-gaJpZM4RKRoQ.

zorino commented 6 years ago

Thank you, I confirm that the last commit fixed the problem with my dataset.

Would it be possible for you to make another release with the fix ?

I'm closing this issue.

wrpearson / fasta36

glsearch36 e-value threshold #2