Thanks a lot !
I now understand better how these two p-values are computed.
However, the p-value that is output in the vcf does not correspond to
either of these two p-values... it actually still corresponds to the
contrast2() function.
I therefore still have the following questions:
1- Is there a paper you could point me to where the computation of the
p-value output in the vcf is described ? Or should it really be
deprecated and replaced by one of the two p-values you pointed to me ?
2- Is it possible to specify the number of individuals per sample when
running bcftools with pooled samples ?
For now, I followed the recommandation on
http://samtools.sourceforge.net/mpileup.shtml, and used the following
command:
samtools mpileup -uf ref.fa all-aln.bam | bcftools view -vcs xxx -1 yyy - > out.vcf
But I guess that if I do not specify the number of individuals per
sample, each sample will be treated as one diploid individual, which is
not what I need...
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.
On 02/11/2015 02:05 PM, Heng Li wrote:
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 *
* *
*********************************************************************
.
--
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 *
* *
*********************************************************************
------------------------------------------------------------------------------
Download BIRT iHub F-Type - The Free Enterprise-Grade BIRT Server
from Actuate! Instantly Supercharge Your Business Reports and Dashboards
with Interactivity, Sharing, Native Excel Exports, App Integration & more
Get technology previously reserved for billion-dollar corporations, FREE
http://pubads.g.doubleclick.net/gampad/clk?id=190641631&iu=/4140/ostg.clktrk
_______________________________________________
Samtools-help mailing list
Samtools-help@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/samtools-help
