1 /* The MIT License
2 
3    Copyright (c) 2016 Hyun Min Kang <hmkang.umich.edu> and Adrian Tan <atks@umich.edu>
4 
5    Permission is hereby granted, free of charge, to any person obtaining a copy
6    of this software and associated documentation files (the "Software"), to deal
7    in the Software without restriction, including without limitation the rights
8    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9    copies of the Software, and to permit persons to whom the Software is
10    furnished to do so, subject to the following conditions:
11 
12    The above copyright notice and this permission notice shall be included in
13    all copies or substantial portions of the Software.
14 
15    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21    THE SOFTWARE.
22 */
23 
24 #ifndef COMPLEX_GENOTYPING_RECORD_H
25 #define COMPLEX_GENOTYPING_RECORD_H
26 
27 #include "htslib/vcf.h"
28 #include "htslib/faidx.h"
29 #include "bcf_ordered_writer.h"
30 #include "variant.h"
31 #include "hts_utils.h"
32 #include "augmented_bam_record.h"
33 #include "estimator.h"
34 #include "genotyping_record.h"
35 
36 /**
37  * A generic record that holds information for genotyping a
38  * variant across multiple samples.
39  *
40  * Maintains read information and allows for additional reads
41  * till VCF record can be printed out.
42  */
43 class ComplexGenotypingRecord : public GenotypingRecord
44 {
45     public:
46     bcf_hdr_t *h;
47     //bcf1_t *v;
48     int32_t rid;
49     int32_t pos1; //position of variant
50     //[beg1,end1] is the required overlapping of the variant against the aligned read necessary to make a genotype call.
51     //for SNPs, beg1=end1=pos1
52     //
53     //for Indels, this refers to the flanking positions
54     //   insertion
55     //        if T/TG - beg1=pos1, end1=pos1+1
56     //        if T/GT - beg1=pos1-1, end1=pos1
57     //   deletion
58     //        if TG/T - beg1=pos1, end1=pos1+length(REF)
59     //        if TG/G - beg1=pos1-1, end1=pos1+length(REF)-1
60     int32_t beg1;
61     int32_t end1;
62     int32_t vtype;
63 
64     //indel specific record
65     int32_t dlen;
66     int32_t len;
67     std::string indel;
68 
69     //vntr specific record
70     std::string motif;
71 
72     //vntr specific record
73     //std::vector<float> counts;
74 
75     // sample level information
76     int32_t nsamples;
77     kstring_t alleles;
78     std::vector<std::string> v_alleles;
79     uint32_t n_filter;
80 
81     uint8_t* pls;
82     uint8_t* ads;
83 
84     // sufficient statistics for computing INFO field
85     float bqr_num, bqr_den;
86     float mqr_num, mqr_den;
87     float cyr_num, cyr_den;
88     float str_num, str_den;
89     float nmr_num, nmr_den;
90     float ior_num, ior_den;
91     float nm0_num, nm0_den;
92     float nm1_num, nm1_den;
93     float abe_num, abe_den;
94     float abz_num, abz_den;
95     float ns_nref, dp_sum, max_gq;
96 
97     int32_t tmp_dp_q20;
98     int32_t tmp_dp_ra;
99     int32_t tmp_bq_s1, tmp_bq_s2;
100     int32_t tmp_mq_s1, tmp_mq_s2;
101     float tmp_cy_s1, tmp_cy_s2;
102     int32_t tmp_st_s1, tmp_st_s2;
103     int32_t tmp_al_s1, tmp_bq_al, tmp_mq_al;
104     float  tmp_cy_al;
105     int32_t tmp_st_al, tmp_nm_al;
106     int32_t tmp_nm_s1, tmp_nm_s2;
107     double tmp_oth_exp_q20, tmp_oth_obs_q20;
108     double tmp_pls[3];
109     double tmp_ads[3];
110 
111     // temporary information to be cleared out per-sample basis
112 
113     /**
114      * Constructor.
115      * @v - VCF record.
116      */
ComplexGenotypingRecord()117     ComplexGenotypingRecord() {};
118 
119     /**
120      * Constructor.
121      * @v - VCF record.
122      */
123     ComplexGenotypingRecord(bcf_hdr_t *h, bcf1_t *v, int32_t vtype, int32_t nsamples);
124 
125     /**
126      * Clears this record.
127      */
128     void clear();
129     void clearTemp();
130     bcf1_t* flush_variant(bcf_hdr_t* hdr);
131     void flush_sample( int32_t sampleIndex );
132     void add_allele( double contam, int32_t allele, uint8_t mapq, bool fwd, uint32_t q, int32_t cycle, uint32_t nm );
133     void process_read(AugmentedBAMRecord& as, int32_t sampleIndex, double contam);
134 
135     /**
136      * Destructor.
137      */
138     ~ComplexGenotypingRecord();
139 };
140 
141 #endif
142