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
