call: Trio constrained calling (no sex chromosomes yet)

junaruga · Oct 1, 2013 · c4172c5 · c4172c5
1 parent 62da94b
commit c4172c5
Show file tree

Hide file tree

Showing 5 changed files with 659 additions and 279 deletions.
diff --git a/Makefile b/Makefile
@@ -13,7 +13,8 @@ DFLAGS=
 OBJS=		main.o vcfview.o bcfidx.o tabix.o \
 			vcfstats.o vcfisec.o vcfmerge.o vcfquery.o vcffilter.o filter.o vcfsom.o \
             vcfnorm.o vcfgtcheck.o vcfsubset.o \
-            vcfcall.o mcall.o
+            vcfcall.o mcall.o \
+            ccall.o em.o prob1.o kmin.o # the original samtools calling
 INCLUDES=	-I. -I$(HTSDIR)
 SUBDIRS=    .
 
@@ -60,6 +61,7 @@ test:
 
 main.o: version.h $(HTSDIR)/version.h bcftools.h
 vcfcall.o: vcfcall.c call.h mcall.c prob1.h $(HTSDIR)/htslib/kfunc.h $(HTSDIR)/htslib/vcf.h
+mcall.o ccall.o: call.h
 vcffilter.o: filter.h
 vcfsubset.o: filter.h
 

diff --git a/call.h b/call.h
@@ -3,37 +3,59 @@
 
 #include <htslib/vcf.h>
 
-#define CALL_KEEPALT 1
-#define CALL_VARONLY (1<<1)
+#define CALL_KEEPALT        1
+#define CALL_VARONLY        (1<<1)
+#define CALL_CONSTR_ALLELES (1<<2)
+#define CALL_CONSTR_TRIO    (1<<3)
 
+#define FATHER 0
+#define MOTHER 1
+#define CHILD  2
+typedef struct
+{
+    char *name;
+    int sample[3];
+}
+family_t;
+
+typedef struct _ccall_t ccall_t;
 typedef struct
 {
     // mcall only
     double min_ma_lrt;  // variant accepted if P(chi^2)>=FLOAT [0.99]
-    int *PLs, nPLs, *gts;   // VCF PL likelihoods (rw); GTs (w)
-    float anno16[16];       // see anno[16] in bam2bcf.h
-    double *pl2p, *pdg;     // PL to 10^(-PL/10) table; PLs converted to P(D|G)
     float *qsum;            // QS(sum) values
     int nqsum, npdg;
     int *als_map, nals_map; // mapping from full set of alleles to trimmed set of alleles (old -> new)
     int *pl_map, npl_map;   // same as above for PLs, but reverse (new -> old)
     char **als;             // array to hold the trimmed set of alleles to appear on output
     int nals;               // size of the als array
+    family_t *fams;         // list of families and samples for trio calling
+    int nfams, mfams;
+    int ntrio[5];           // possible trio genotype combinations
+    uint16_t *trio[5];
+    double *GLs, *sumGLs;
+    int *GQs;               // VCF FORMAT genotype qualities
 
     // ccall only
     double indel_frac, theta, min_lrt, min_perm_p; 
     double prior_type, pref;
+    double ref_lk, lk_sum;
     int ngrp1_samples, n_perm;
+    int ac[4], nhets, ndiploid;
     char *prior_file;
+    ccall_t *cdat;
 
     // shared
     bcf1_t *rec;
     bcf_hdr_t *hdr;
-
-    uint32_t flag;      // One or more of the CALL_* flags defined above
+    uint32_t flag;          // One or more of the CALL_* flags defined above
     uint8_t *ploidy;
-    int nsamples;
-    uint32_t trio;
+
+    double pl2p[256];       // PL to 10^(-PL/10) table
+    int *PLs, nPLs, mPLs;   // VCF PL likelihoods (rw)
+    int *gts;               // GTs (w)
+    double *pdg;            // PLs converted to P(D|G)
+    float *anno16; int n16; // see anno[16] in bam2bcf.h
 }
 call_t;
 
@@ -55,5 +77,7 @@ void mcall_destroy(call_t *call);
 void ccall_destroy(call_t *call);
 void qcall_destroy(call_t *call);
 
+void call_init_pl2p(call_t *call);
+uint32_t *call_trio_prep(int is_x, int is_son);
 
 #endif