1 #include "bcftools.pysam.h"
2
3 /* gvcf.c -- support for gVCF files.
4
5 Copyright (C) 2014-2015 Genome Research Ltd.
6
7 Author: Petr Danecek <pd3@sanger.ac.uk>
8
9 Permission is hereby granted, free of charge, to any person obtaining a copy
10 of this software and associated documentation files (the "Software"), to deal
11 in the Software without restriction, including without limitation the rights
12 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 copies of the Software, and to permit persons to whom the Software is
14 furnished to do so, subject to the following conditions:
15
16 The above copyright notice and this permission notice shall be included in
17 all copies or substantial portions of the Software.
18
19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 DEALINGS IN THE SOFTWARE. */
26
27 #include "gvcf.h"
28 #include "bcftools.h"
29
30 struct _gvcf_t
31 {
32 int *dp_range, ndp_range; // per-sample DP ranges
33 int prev_range; // 0 if not in a block
34 int32_t *dp, mdp, *pl, mpl, npl;
35 int32_t *tmp, mtmp, *gts, ngts,mgts, nqsum,mqsum;
36 float *qsum;
37 int32_t rid, start, end, min_dp;
38 kstring_t als;
39 bcf1_t *line;
40 };
41
gvcf_update_header(gvcf_t * gvcf,bcf_hdr_t * hdr)42 void gvcf_update_header(gvcf_t *gvcf, bcf_hdr_t *hdr)
43 {
44 bcf_hdr_append(hdr,"##INFO=<ID=END,Number=1,Type=Integer,Description=\"End position of the variant described in this record\">");
45 bcf_hdr_append(hdr,"##INFO=<ID=MinDP,Number=1,Type=Integer,Description=\"Minimum per-sample depth in this gVCF block\">");
46 }
47
gvcf_init(const char * dp_ranges)48 gvcf_t *gvcf_init(const char *dp_ranges)
49 {
50 gvcf_t *gvcf = (gvcf_t*) calloc(1,sizeof(gvcf_t));
51 gvcf->line = bcf_init();
52
53 int n = 1;
54 const char *ss = dp_ranges;
55 while ( *ss )
56 {
57 if ( *ss==',' ) n++;
58 ss++;
59 }
60 gvcf->ndp_range = n;
61 gvcf->dp_range = (int*) malloc(sizeof(int)*gvcf->ndp_range);
62
63 n = 0;
64 ss = dp_ranges;
65 while ( *ss )
66 {
67 char *se = (char*) ss;
68 gvcf->dp_range[n++] = strtol(ss,&se,10);
69 if ( se==ss ) return NULL;
70 if ( *se==',' && se[1] ) { ss = se+1; continue; }
71 else if ( !*se ) break;
72 return NULL;
73 }
74 return gvcf;
75 }
76
gvcf_destroy(gvcf_t * gvcf)77 void gvcf_destroy(gvcf_t *gvcf)
78 {
79 free(gvcf->dp_range);
80 free(gvcf->dp);
81 free(gvcf->pl);
82 free(gvcf->tmp);
83 free(gvcf->qsum);
84 free(gvcf->gts);
85 free(gvcf->als.s);
86 if ( gvcf->line ) bcf_destroy(gvcf->line);
87 free(gvcf);
88 }
89
gvcf_write(gvcf_t * gvcf,htsFile * fh,bcf_hdr_t * hdr,bcf1_t * rec,int is_ref)90 bcf1_t *gvcf_write(gvcf_t *gvcf, htsFile *fh, bcf_hdr_t *hdr, bcf1_t *rec, int is_ref)
91 {
92 int i, ret, nsmpl = bcf_hdr_nsamples(hdr);
93 int can_collapse = is_ref ? 1 : 0;
94 int32_t dp_range = 0, min_dp = 0;
95
96 // No record and nothing to flush?
97 if ( !rec && !gvcf->prev_range ) return NULL;
98
99 // Flush gVCF block if there are no more records, chr changed, a gap
100 // encountered, or other conditions not met (block broken by a non-ref or DP too low).
101 int needs_flush = can_collapse ? 0 : 1;
102
103
104 // Can the record be included in a gVCF block? That is, is this a ref-only site?
105 if ( rec && can_collapse )
106 {
107 bcf_unpack(rec, BCF_UN_ALL);
108
109 // per-sample depth
110 ret = bcf_get_format_int32(hdr, rec, "DP", &gvcf->tmp, &gvcf->mtmp);
111 if ( ret==nsmpl )
112 {
113 min_dp = gvcf->tmp[0];
114 for (i=1; i<nsmpl; i++)
115 if ( min_dp > gvcf->tmp[i] ) min_dp = gvcf->tmp[i];
116
117 for (i=0; i<gvcf->ndp_range; i++)
118 if ( min_dp < gvcf->dp_range[i] ) break;
119
120 dp_range = i;
121 if ( !dp_range )
122 {
123 // leave the record unchanged, DP is too small. Alternatively, return NULL here
124 // to skip these sites
125 needs_flush = 1;
126 can_collapse = 0;
127 }
128 }
129 else
130 needs_flush = 1; // DP field not present
131 }
132
133 if ( gvcf->prev_range && gvcf->prev_range!=dp_range ) needs_flush = 1;
134 if ( !rec || gvcf->rid!=rec->rid || rec->pos > gvcf->end+1 ) needs_flush = 1;
135
136 // If prev_range is set, something can be flushed
137 if ( gvcf->prev_range && needs_flush )
138 {
139 // mpileup can output two records with the same position, SNP and
140 // indel. Make sure the end position does not include the non-variant
141 // SNP position just before the indel.
142 if ( rec && rec->rid==gvcf->rid && rec->pos==gvcf->end ) gvcf->end--;
143
144 gvcf->end++; // from 0-based to 1-based coordinate
145
146 bcf_clear1(gvcf->line);
147 gvcf->line->rid = gvcf->rid;
148 gvcf->line->pos = gvcf->start;
149 gvcf->line->rlen = gvcf->end - gvcf->start;
150 bcf_update_alleles_str(hdr, gvcf->line, gvcf->als.s);
151 if ( gvcf->start+1 < gvcf->end ) // create gVCF record only if it spans at least two sites
152 bcf_update_info_int32(hdr, gvcf->line, "END", &gvcf->end, 1);
153 bcf_update_info_int32(hdr, gvcf->line, "MinDP", &gvcf->min_dp, 1);
154 if ( gvcf->nqsum>0 )
155 bcf_update_info_float(hdr, gvcf->line, "QS", gvcf->qsum, gvcf->nqsum);
156 if ( gvcf->ngts )
157 bcf_update_genotypes(hdr,gvcf->line,gvcf->gts,gvcf->ngts);
158 if ( gvcf->npl>0 )
159 bcf_update_format_int32(hdr, gvcf->line, "PL", gvcf->pl, gvcf->npl);
160 bcf_update_format_int32(hdr, gvcf->line, "DP", gvcf->dp, nsmpl);
161 if ( bcf_write1(fh, hdr, gvcf->line)!=0 ) error("[%s] Error: failed to write the record\n", __func__);
162 gvcf->prev_range = 0;
163 gvcf->rid = -1;
164 gvcf->npl = 0;
165 gvcf->nqsum = 0;
166 gvcf->ngts = 0;
167
168 if ( !rec ) return NULL; // just flushing the buffer, this was last record
169 }
170
171 if ( can_collapse )
172 {
173 if ( !gvcf->prev_range )
174 {
175 hts_expand(int32_t,nsmpl,gvcf->mdp,gvcf->dp);
176 memcpy(gvcf->dp,gvcf->tmp,nsmpl*sizeof(int32_t)); // tmp still contains DP from rec
177 gvcf->npl = bcf_get_format_int32(hdr, rec, "PL", &gvcf->pl, &gvcf->mpl);
178
179 gvcf->nqsum = bcf_get_info_float(hdr,rec,"QS",&gvcf->qsum,&gvcf->mqsum);
180 gvcf->ngts = bcf_get_genotypes(hdr,rec,&gvcf->gts,&gvcf->mgts);
181
182 gvcf->rid = rec->rid;
183 gvcf->start = rec->pos;
184 gvcf->als.l = 0;
185 kputs(rec->d.allele[0],&gvcf->als);
186 for (i=1; i<rec->n_allele; i++)
187 {
188 kputc(',',&gvcf->als);
189 kputs(rec->d.allele[i],&gvcf->als);
190 }
191 gvcf->min_dp = min_dp;
192 }
193 else
194 {
195 if ( gvcf->min_dp > min_dp ) gvcf->min_dp = min_dp;
196 for (i=0; i<nsmpl; i++)
197 if ( gvcf->dp[i] > gvcf->tmp[i] ) gvcf->dp[i] = gvcf->tmp[i];
198 ret = bcf_get_format_int32(hdr, rec, "PL", &gvcf->tmp, &gvcf->mtmp);
199 if ( ret>=0 )
200 {
201 if ( ret!=nsmpl*3 ) error("Unexpected number of PL fields\n");
202 for (i=0; i<nsmpl; i++)
203 {
204 if ( gvcf->pl[3*i+1] > gvcf->tmp[3*i+1] )
205 {
206 gvcf->pl[3*i+1] = gvcf->tmp[3*i+1];
207 gvcf->pl[3*i+2] = gvcf->tmp[3*i+2];
208 }
209 else if ( gvcf->pl[3*i+1]==gvcf->tmp[3*i+1] && gvcf->pl[3*i+2] > gvcf->tmp[3*i+2] )
210 gvcf->pl[3*i+2] = gvcf->tmp[3*i+2];
211 }
212 }
213 else
214 gvcf->npl = 0;
215 }
216 gvcf->prev_range = dp_range;
217 if ( bcf_get_info_int32(hdr,rec,"END",&gvcf->tmp,&gvcf->mtmp)==1 )
218 gvcf->end = gvcf->tmp[0] - 1; // from 1-based to 0-based
219 else
220 gvcf->end = rec->pos;
221 return NULL;
222 }
223
224 if ( is_ref && min_dp )
225 bcf_update_info_int32(hdr, rec, "MinDP", &min_dp, 1);
226
227 return rec;
228 }
229
230