Thanks a lot !
I now understand better how these two p-values are computed.
Indeed, the second is far more conservative than the first, as stated in your
paper.
However, the one that is output in the vcf does not correspond to either of
these two p-values.
In fact, it is still computed using the contrast2() function.
My question is therefore: which p-value would you advise me to use ?
Should I follow the advice given in your paper and use the first p-value
(computed with D_alpha1) and filter out the cases with low D_alpha2 ?
Thanks in advance !
Vincent
PS: I am using samtools v0.1.19, which is the latest version in which you can
specify that you want to compare two groups.
________________________________________
De : Heng Li [hen...@broadinstitute.org]
Date d'envoi : mercredi 11 février 2015 14:05
À : LACROIX VINCENT
Cc: Samtools-help List
Objet : Re: Association Test with mpileup
I forgot which tag gives the P-value, but you should not look at this part of
the code. It is conservative and should be actually deprecated. You should look
at this function:
// x[0]: ref frequency
// x[1..3]: alt-alt, alt-ref, ref-ref frequenc
// x[4]: HWE P-value
// x[5..6]: group1 freq, group2 freq.
// x[7]: 1-degree P-value
// x[8]: 2-degree P-value
int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
The P-values are returned at x[7] and x[8]. See the equation (11) and (12) in
the samtools paper for the computation.
Heng
On Feb 11, 2015, at 6:55, Vincent Lacroix <vincent.lacr...@univ-lyon1.fr> wrote:
> Dear Heng Li,
>
> I am using mpileup to call SNPs from pooled RNAseq data.
> I have two conditions and I wanted to find SNPs that are condition specific.
> I have tried to use the association test available in mpileup, but I am
> having a hard time understanding why it is so conservative.
> I see in the documentation that what you output is the posterior weighted
> chi2 P-value, but I am not familiar with this concept.
> I think I have isolated the chunk of code which corresponds to this
> computation (reproduced below), but I am having a hard time understanding
> what "p" corresponds to...
> Could you give me an explanation ? Or point me to a reference where this is
> explained ?
>
> Thanks a lot in advance,
>
> Vincent Lacroix
>
>
> // chi2=(a+b+c+d)(ad-bc)^2/[(a+b)(c+d)(a+c)(b+d)]
> static inline double contrast2_aux(const bcf_p1aux_t *p1, double sum, int k1,
> int k2, double x[3])
> {
> double p = p1->phi[k1+k2] * p1->z1[k1] * p1->z2[k2] / sum * p1->hg[k1][k2];
> int n1 = p1->n1, n2 = p1->n - p1->n1;
> if (p < CONTRAST_TINY) return -1;
> if (.5*k1/n1 < .5*k2/n2) x[1] += p;
> else if (.5*k1/n1 > .5*k2/n2) x[2] += p;
> else x[0] += p;
> return p * chi2_test(k1, k2, (n1<<1) - k1, (n2<<1) - k2);
> }
>
> --
> Vincent Lacroix
> Université Claude Bernard - Lyon I
> UMR CNRS 5558 Biometrie et Biologie Évolutive
> INRIA Projet BAMBOO
> http://lbbe.univ-lyon1.fr/-Lacroix-Vincent-.html
>
> ************** ATTENTION NOUVELLE ADRESSE ELECTRONIQUE **************
> * *
> * vincent.lacr...@univ-lyon1.fr *
> * *
> *********************************************************************
>
------------------------------------------------------------------------------
Dive into the World of Parallel Programming. The Go Parallel Website,
sponsored by Intel and developed in partnership with Slashdot Media, is your
hub for all things parallel software development, from weekly thought
leadership blogs to news, videos, case studies, tutorials and more. Take a
look and join the conversation now. http://goparallel.sourceforge.net/
_______________________________________________
Samtools-help mailing list
Samtools-help@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/samtools-help
