PageSourceSearch

https://wresch.github.io/2014/11/18/process-vcf-file-with-htslib.html

html wresch.github.io collected 2026-10-03 10:10:56 UTC 32,539 bytes, 328 lines download raw bytes

1<!DOCTYPE html>
2<html lang="en">
3  <head>
4    <meta charset="utf-8">
5    <title>Process a VCF file with htslib</title>
6    <link rel="stylesheet" href="/css/syntax.css">
7    <link rel="stylesheet" href="/css/main.css">
8    <link href="http://fonts.googleapis.com/css?family=Inconsolata" rel="stylesheet" type="text/css">
9    <!--[if lt IE 9]>
10	
10<script src="js/html5shiv.js"></script>
10
11    <![endif]-->
12    
12<script type="text/javascript"
13	    src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML">
14    </script>
14
15    
15<script type="text/x-mathjax-config">
16      MathJax.Hub.Config({
17        tex2jax: {
18          skipTags: ["script", "noscript", "style", "textarea", "pre"]
19        }
20      });
21    </script>
21
22    
22<script>
vendor: 339 bytes, lines 22-28
22
23  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
24  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
25  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
26  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');
27
28  ga('create', '
28UA-52346637-1
vendor: 39 bytes, lines 28-30
28', 'auto');
29  ga('send', 'pageview');
30
31</script>
31
32
33  </head>
34
35<body>
36  <a id="page-top"></a>
37    <header class="site-hdr">
38    <h1><a href="/">Wolfgang Resch - Notes</a></h1>
39    <nav>
40      <a href="/">Home</a>
41      <a href="/publications.html">Publications</a>
42    </nav>
43  </header>
44
45  <article>
46    <header class="article-hdr">
47      <h1>Process a VCF file with htslib</h1>
48      <span class="date">November 18, 2014</span>
49    </header>
50    <hr class="light">
51    <p><a href="http://www.htslib.org/">Samtools and Bcftools</a> are migrating to a single underlying
52library for dealing with the various high-throughput sequencing data formats (SAM,
53BAM, CRAM, VCF, BCF, and tabix). The library is called
54<a href="https://github.com/samtools/htslib">htslib</a>.  As of now, documentation is limited
55to comments in the code and examples can be found in the newer revisions of
56samtools and bcftools using htslib.  Htslib has minimal dependencies (zlib) and
57can easily be downloaded and installed.  The makefile creates a static and a
58dynamic library.</p>
59
60<p>I wanted to filter the <a href="ftp://ftp-mouse.sanger.ac.uk/REL-1410-SNPs_Indels/">VCF file</a> from the
61<a href="http://www.sanger.ac.uk/resources/mouse/genomes/">Sanger mouse genome project</a>. The
62goal was to extract high quality SNPs for a single mouse strain in a tabular
63format for further processing downstream (in my case create a mm10 genome incorporating
64all the SNPs from one strain).</p>
65
66<h3 id="open-vcf-or-bcf-file-and-extract-the-number-of-samples-present-and-the-sequence-names">Open VCF (or BCF) file and extract the number of samples present and the sequence names</h3>
67
68<div class="language-c highlighter-rouge"><pre class="highlight"><code><span class="cp">#include &lt;stdio.h&gt;  //puts and printf
69#include &lt;stdlib.h&gt; //EXIT_FAILURE 
70#include "vcf.h"
71</span>
72<span class="kt">int</span> <span class="nf">main</span><span class="p">(</span><span class="kt">int</span> <span class="n">argc</span><span class="p">,</span> <span class="kt">char</span> <span class="o">**</span><span class="n">argv</span><span class="p">)</span> <span class="p">{</span>
73        <span class="k">if</span> <span class="p">(</span><span class="n">argc</span> <span class="o">!=</span> <span class="mi">2</span><span class="p">)</span> <span class="p">{</span>
74                <span class="k">return</span> <span class="n">EXIT_FAILURE</span><span class="p">;</span>
75        <span class="p">}</span>
76
77        <span class="c1">// counters
78</span>        <span class="kt">int</span> <span class="n">nseq</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
79
80        <span class="c1">// open VCF/BCF file
81</span>        <span class="c1">//    * use '-' for stdin
82</span>        <span class="c1">//    * bcf_open will open bcf and vcf files
83</span>        <span class="c1">//    * bcf_open is a macro that expands to hts_open
84</span>        <span class="c1">//    * returns NULL when file could not be opened
85</span>        <span class="c1">//    * by default also writes message to stderr if file could not be found
86</span>        <span class="n">htsFile</span> <span class="o">*</span> <span class="n">inf</span> <span class="o">=</span> <span class="n">bcf_open</span><span class="p">(</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="s">"r"</span><span class="p">);</span>
87        <span class="k">if</span> <span class="p">(</span><span class="n">inf</span> <span class="o">==</span> <span class="nb">NULL</span><span class="p">)</span> <span class="p">{</span>
88                <span class="k">return</span> <span class="n">EXIT_FAILURE</span><span class="p">;</span>
89        <span class="p">}</span>
90
91        <span class="c1">// read header
92</span>        <span class="n">bcf_hdr_t</span> <span class="o">*</span><span class="n">hdr</span> <span class="o">=</span> <span class="n">bcf_hdr_read</span><span class="p">(</span><span class="n">inf</span><span class="p">);</span>
93        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"File %s contains %i samples</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">bcf_hdr_nsamples</span><span class="p">(</span><span class="n">hdr</span><span class="p">));</span>
94
95        <span class="c1">// report names of all the sequences in the VCF file
96</span>        <span class="k">const</span> <span class="kt">char</span> <span class="o">**</span><span class="n">seqnames</span> <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
97        <span class="c1">// bcf_hdr_seqnames returns a newly allocated array of pointers to the seq names
98</span>        <span class="c1">// caller has to deallocate the array, but not the seqnames themselves; the number
99</span>        <span class="c1">// of sequences is stored in the int pointer passed in as the second argument.
100</span>        <span class="c1">// The id in each record can be used to index into the array to obtain the sequence
101</span>        <span class="c1">// name
102</span>        <span class="n">seqnames</span> <span class="o">=</span> <span class="n">bcf_hdr_seqnames</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">nseq</span><span class="p">);</span>
103        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"Sequence names:</span><span class="se">\n</span><span class="s">"</span><span class="p">);</span>
104        <span class="k">for</span> <span class="p">(</span><span class="kt">int</span> <span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o">&lt;</span> <span class="n">nseq</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
105                <span class="c1">// bcf_hdr_id2name is another way to get the name of a sequence
106</span>                <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"  [%2i] %s (bcf_hdr_id2name -&gt; %s)</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">seqnames</span><span class="p">[</span><span class="n">i</span><span class="p">],</span>
107                       <span class="n">bcf_hdr_id2name</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">i</span><span class="p">));</span>
108        <span class="p">}</span>
109
110
111        <span class="c1">// clean up memory
112</span>        <span class="k">if</span> <span class="p">(</span><span class="n">seqnames</span> <span class="o">!=</span> <span class="nb">NULL</span><span class="p">)</span>
113                <span class="n">free</span><span class="p">(</span><span class="n">seqnames</span><span class="p">);</span>
114        <span class="n">bcf_hdr_destroy</span><span class="p">(</span><span class="n">hdr</span><span class="p">);</span>
115        <span class="n">bcf_close</span><span class="p">(</span><span class="n">inf</span><span class="p">);</span>
116        <span class="k">return</span> <span class="mi">0</span><span class="p">;</span>
117<span class="p">}</span>
118
119</code></pre>
120</div>
121
122<p>With the comments this should be pretty self explanatory.  Internally, htslib represents
123records as <code class="highlighter-rouge">bcf1_t</code> structures, which is why most the functions use the <code class="highlighter-rouge">bcf_</code> prefix.
124The <code class="highlighter-rouge">bcf_</code> functions appear to be general and work with vcf and bcf format files. For
125example, <code class="highlighter-rouge">bcf_open</code> will open vcf and bcf files (and actually any of the other formats
126since it simply points to <code class="highlighter-rouge">hts_open</code>, the function that is used to open all file types.
127Even though vcf uses 1-based indexing (i.e. first base is base 1), htslib internally uses
1280-based indexing (i.e. bcf1_t::pos is 0 based).</p>
129
130<h3 id="iterate-through-all-snps-in-file-and-extract-data">Iterate through all SNPs in file and extract data</h3>
131
132<p>Next, I want to iterate through all the SNPs for one particular mouse strain, filter
133to high quality SNPs and output the data in tabular form.  It seems that this can
134be accomplished with the <code class="highlighter-rouge">bcf_read</code> function, which will read the next record and
135return 0 on success. When reading vcf files (as done here), it is necessary to
136call <code class="highlighter-rouge">bcf_unpack</code> after read to populate the <code class="highlighter-rouge">bcf1_t::d</code> field.  However, each
137of the function for fetching values from samples calles <code class="highlighter-rouge">bcf_unpack</code> if it hasn’t
138already been called, so no explicit call is required here.  Note that unpacking
139the data for each gt for each sample is time comsuming for vcf data. Therefore,
140if only a subset of samples is going to be considered, <code class="highlighter-rouge">bcf_hdr_set_samples</code>
141can be used to limit which samples are parsed. It can be given a single sample or
142a list of comma separated samples to include or exclude (prefix with <code class="highlighter-rouge">^</code>). It
143can also take a filename for a file containing the information.</p>
144
145<p>The <code class="highlighter-rouge">bcf_get_format_*()</code> functions extract data from the genotypes
146for each record. This is subject to the restrictions imposed with
147<code class="highlighter-rouge">bcf_hdr_set_samples</code>.  These functions allocate new memory if necessary. In
148the case of the code here, they will only allocate on the first call and then
149re-use the memory passed in, since all records contain the same number of
150samples.</p>
151
152<p>The actual filters applied here check for high quality, homozygous ALT calls
153(FI == 1) with a genotype quality score &gt; 20. See the
154<a href="ftp://ftp-mouse.sanger.ac.uk/REL-1410-SNPs_Indels/README">sanger site</a> for
155more details.</p>
156
157<div class="language-c highlighter-rouge"><pre class="highlight"><code><span class="cp">#include &lt;stdio.h&gt;
158#include "vcf.h"
159#include "vcfutils.h"
160</span>
161<span class="kt">void</span> <span class="nf">usage</span><span class="p">()</span> <span class="p">{</span>
162        <span class="n">puts</span><span class="p">(</span>
163        <span class="s">"NAME</span><span class="se">\n</span><span class="s">"</span>
164        <span class="s">"    03_vcf - High quality calls for single sample</span><span class="se">\n</span><span class="s">"</span>
165        <span class="s">"SYNOPSIS</span><span class="se">\n</span><span class="s">"</span>
166        <span class="s">"    03_vcf vcf_file sample</span><span class="se">\n</span><span class="s">"</span>
167        <span class="s">"DESCRIPTION</span><span class="se">\n</span><span class="s">"</span>
168        <span class="s">"    Given a &lt;vcf file&gt;, extract all calls for &lt;sample&gt; and filter for</span><span class="se">\n</span><span class="s">"</span>
169        <span class="s">"    high quality, homozygous SNPs. This will omit any positions</span><span class="se">\n</span><span class="s">"</span>
170        <span class="s">"    that are homozygous ref (0/0) or heterozygous.  The exact filter</span><span class="se">\n</span><span class="s">"</span>
171        <span class="s">"    used is</span><span class="se">\n</span><span class="s">"</span>
172        <span class="s">"        FI == 1 &amp; GQ &gt; 20 &amp; GT != '0/0'</span><span class="se">\n</span><span class="s">"</span>
173        <span class="s">"        [NOTE: FI == 1 implies homozygous call]</span><span class="se">\n</span><span class="s">"</span>
174        <span class="s">"        </span><span class="se">\n</span><span class="s">"</span>
175        <span class="s">"    The returned format is</span><span class="se">\n</span><span class="s">"</span>
176        <span class="s">"        chrom pos[0-based]  REF  ALT GQ|DP</span><span class="se">\n</span><span class="s">"</span>
177        <span class="s">"    and can be used for Marei's personalizer.py</span><span class="se">\n</span><span class="s">"</span>
178                <span class="p">);</span>
179<span class="p">}</span>
180
181
182
183<span class="kt">int</span> <span class="nf">main</span><span class="p">(</span><span class="kt">int</span> <span class="n">argc</span><span class="p">,</span> <span class="kt">char</span> <span class="o">**</span><span class="n">argv</span><span class="p">)</span> <span class="p">{</span>
184        <span class="k">if</span> <span class="p">(</span><span class="n">argc</span> <span class="o">!=</span> <span class="mi">3</span><span class="p">)</span> <span class="p">{</span>
185                <span class="n">usage</span><span class="p">();</span>
186                <span class="k">return</span> <span class="mi">1</span><span class="p">;</span>
187        <span class="p">}</span>
188        <span class="c1">// counters
189</span>        <span class="kt">int</span> <span class="n">n</span>    <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>  <span class="c1">// total number of records in file
190</span>        <span class="kt">int</span> <span class="n">nsnp</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>  <span class="c1">// number of SNP records in file
191</span>        <span class="kt">int</span> <span class="n">nhq</span>  <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>  <span class="c1">// number of SNPs for the single sample passing filters
192</span>        <span class="kt">int</span> <span class="n">nseq</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>  <span class="c1">// number of sequences
193</span>        <span class="c1">// filter data for each call
194</span>        <span class="kt">int</span> <span class="n">nfi_arr</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
195        <span class="kt">int</span> <span class="n">nfi</span>     <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
196        <span class="kt">int</span> <span class="o">*</span><span class="n">fi</span>     <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
197        <span class="c1">// quality data for each call
198</span>        <span class="kt">int</span> <span class="n">ngq_arr</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
199        <span class="kt">int</span> <span class="n">ngq</span>     <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
200        <span class="kt">int</span> <span class="o">*</span><span class="n">gq</span>     <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
201        <span class="c1">// coverage data for each call
202</span>        <span class="kt">int</span> <span class="n">ndp_arr</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
203        <span class="kt">int</span> <span class="n">ndp</span>     <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
204        <span class="kt">int</span> <span class="o">*</span><span class="n">dp</span>     <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
205        <span class="c1">// genotype data for each call
206</span>        <span class="c1">// genotype arrays are twice as large as
207</span>        <span class="c1">// the other arrays as there are two values for each sample
208</span>        <span class="kt">int</span> <span class="n">ngt_arr</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
209        <span class="kt">int</span> <span class="n">ngt</span>     <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
210        <span class="kt">int</span> <span class="o">*</span><span class="n">gt</span>     <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
211        
212        <span class="c1">// open VCF/BCF file
213</span>        <span class="c1">//    * use '-' for stdin
214</span>        <span class="c1">//    * bcf_open will open bcf and vcf files
215</span>        <span class="c1">//    * bcf_open is a macro that expands to hts_open
216</span>        <span class="c1">//    * returns NULL when file could not be opened
217</span>        <span class="c1">//    * by default also writes message to stderr if file could not be found
218</span>        <span class="n">htsFile</span> <span class="o">*</span> <span class="n">inf</span> <span class="o">=</span> <span class="n">bcf_open</span><span class="p">(</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="s">"r"</span><span class="p">);</span>
219        <span class="k">if</span> <span class="p">(</span><span class="n">inf</span> <span class="o">==</span> <span class="nb">NULL</span><span class="p">)</span> <span class="p">{</span>
220                <span class="k">return</span> <span class="n">EXIT_FAILURE</span><span class="p">;</span>
221        <span class="p">}</span>
222        
223        <span class="c1">// read header
224</span>        <span class="n">bcf_hdr_t</span> <span class="o">*</span><span class="n">hdr</span> <span class="o">=</span> <span class="n">bcf_hdr_read</span><span class="p">(</span><span class="n">inf</span><span class="p">);</span>
225        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"File %s contains %i samples</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">bcf_hdr_nsamples</span><span class="p">(</span><span class="n">hdr</span><span class="p">));</span>
226        <span class="c1">// report names of all the sequences in the VCF file
227</span>        <span class="k">const</span> <span class="kt">char</span> <span class="o">**</span><span class="n">seqnames</span> <span class="o">=</span> <span class="nb">NULL</span><span class="p">;</span>
228        <span class="c1">// bcf_hdr_seqnames returns a newly allocated array of pointers to the seq names
229</span>
229        <span class="c1">// caller has to deallocate the array, but not the seqnames themselves; the number
230</span>        <span class="c1">// of sequences is stored in the int pointer passed in as the second argument.
231</span>        <span class="c1">// The id in each record can be used to index into the array to obtain the sequence
232</span>        <span class="c1">// name
233</span>        <span class="n">seqnames</span> <span class="o">=</span> <span class="n">bcf_hdr_seqnames</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">nseq</span><span class="p">);</span>
234        <span class="k">if</span> <span class="p">(</span><span class="n">seqnames</span> <span class="o">==</span> <span class="nb">NULL</span><span class="p">)</span> <span class="p">{</span>
235                <span class="k">goto</span> <span class="n">error1</span><span class="p">;</span>
236        <span class="p">}</span>
237        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"Sequence names:</span><span class="se">\n</span><span class="s">"</span><span class="p">);</span>
238        <span class="k">for</span> <span class="p">(</span><span class="kt">int</span> <span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o">&lt;</span> <span class="n">nseq</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
239                <span class="c1">// bcf_hdr_id2name is another way to get the name of a sequence
240</span>                <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"  [%2i] %s (bcf_hdr_id2name -&gt; %s)</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">seqnames</span><span class="p">[</span><span class="n">i</span><span class="p">],</span>
241                       <span class="n">bcf_hdr_id2name</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">i</span><span class="p">));</span>
242        <span class="p">}</span>
243
244        <span class="c1">// limit the VCF data to the sample name passed in
245</span>        <span class="n">bcf_hdr_set_samples</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">argv</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="mi">0</span><span class="p">);</span>
246        <span class="k">if</span> <span class="p">(</span><span class="n">bcf_hdr_nsamples</span><span class="p">(</span><span class="n">hdr</span><span class="p">)</span> <span class="o">!=</span> <span class="mi">1</span><span class="p">)</span> <span class="p">{</span>
247                <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"ERROR: please limit to a single sample</span><span class="se">\n</span><span class="s">"</span><span class="p">);</span>
248                <span class="k">goto</span> <span class="n">error2</span><span class="p">;</span>
249        <span class="p">}</span>
250
251        <span class="c1">// struc for storing each record
252</span>        <span class="n">bcf1_t</span> <span class="o">*</span><span class="n">rec</span> <span class="o">=</span> <span class="n">bcf_init</span><span class="p">();</span>
253        <span class="k">if</span> <span class="p">(</span><span class="n">rec</span> <span class="o">==</span> <span class="nb">NULL</span><span class="p">)</span> <span class="p">{</span>
254                <span class="k">goto</span> <span class="n">error2</span><span class="p">;</span>
255        <span class="p">}</span>
256        
257        <span class="k">while</span> <span class="p">(</span><span class="n">bcf_read</span><span class="p">(</span><span class="n">inf</span><span class="p">,</span> <span class="n">hdr</span><span class="p">,</span> <span class="n">rec</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="p">{</span>
258                <span class="n">n</span><span class="o">++</span><span class="p">;</span>
259                <span class="k">if</span> <span class="p">(</span><span class="n">bcf_is_snp</span><span class="p">(</span><span class="n">rec</span><span class="p">))</span> <span class="p">{</span>
260                        <span class="n">nsnp</span><span class="o">++</span><span class="p">;</span>
261                <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
262                        <span class="k">continue</span><span class="p">;</span>
263                <span class="p">}</span>
264                <span class="c1">// the bcf_get_format_int32 function does not appear to reallocate
265</span>                <span class="c1">// the array it returns for each of the samples on each call. Just
266</span>                <span class="c1">// needs to be freed in the end. First call to bcf_get_format_*
267</span>                <span class="c1">// takes care of calling bcf_unpack, which fills the `d` member
268</span>                <span class="c1">// of bcf1_t
269</span>                <span class="n">nfi</span> <span class="o">=</span> <span class="n">bcf_get_format_int32</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">rec</span><span class="p">,</span> <span class="s">"FI"</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">fi</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">nfi_arr</span><span class="p">);</span>
270                <span class="c1">// GQ can be missing (".") in this VCF file; The htslib version
271</span>                <span class="c1">// used right now does not return a negative value in that case,
272</span>                <span class="c1">// so we can't check for it.  As it turns out, all homozygous
273</span>                <span class="c1">// good quality calls for alt allele have GQ values, so it doesn't matter
274</span>                <span class="c1">// here, but it's important to keep in mind.
275</span>                <span class="n">ngq</span> <span class="o">=</span> <span class="n">bcf_get_format_int32</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">rec</span><span class="p">,</span> <span class="s">"GQ"</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">gq</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">ngq_arr</span><span class="p">);</span>
276                <span class="n">ndp</span> <span class="o">=</span> <span class="n">bcf_get_format_int32</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">rec</span><span class="p">,</span> <span class="s">"DP"</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">dp</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">ndp_arr</span><span class="p">);</span>
277                <span class="n">ngt</span> <span class="o">=</span> <span class="n">bcf_get_format_int32</span><span class="p">(</span><span class="n">hdr</span><span class="p">,</span> <span class="n">rec</span><span class="p">,</span> <span class="s">"GT"</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">gt</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">ngt_arr</span><span class="p">);</span>
278                <span class="k">if</span> <span class="p">(</span><span class="n">fi</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">==</span> <span class="mi">1</span> <span class="o">&amp;&amp;</span> <span class="n">gq</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mi">20</span> <span class="o">&amp;&amp;</span> <span class="n">gt</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">0</span> <span class="o">&amp;&amp;</span> <span class="n">gt</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">)</span> <span class="p">{</span>
279                        <span class="n">nhq</span><span class="o">++</span><span class="p">;</span>
280                        <span class="n">printf</span><span class="p">(</span><span class="s">"chr%s</span><span class="se">\t</span><span class="s">%i</span><span class="se">\t</span><span class="s">%s</span><span class="se">\t</span><span class="s">%s</span><span class="se">\t</span><span class="s">%i|%i</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">seqnames</span><span class="p">[</span><span class="n">rec</span><span class="o">-&gt;</span><span class="n">rid</span><span class="p">],</span>
281                               <span class="n">rec</span><span class="o">-&gt;</span><span class="n">pos</span><span class="p">,</span>
282                               <span class="n">rec</span><span class="o">-&gt;</span><span class="n">d</span><span class="p">.</span><span class="n">allele</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span>
283                               <span class="n">rec</span><span class="o">-&gt;</span><span class="n">d</span><span class="p">.</span><span class="n">allele</span><span class="p">[</span><span class="n">bcf_gt_allele</span><span class="p">(</span><span class="n">gt</span><span class="p">[</span><span class="mi">0</span><span class="p">])],</span>
284                               <span class="n">gq</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">dp</span><span class="p">[</span><span class="mi">0</span><span class="p">]);</span>
285                <span class="p">}</span>
286                
287        <span class="p">}</span>
288        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"Read %i records %i of which were SNPs</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">nsnp</span><span class="p">);</span>
289        <span class="n">fprintf</span><span class="p">(</span><span class="n">stderr</span><span class="p">,</span> <span class="s">"%i records for the selected sample were high quality homozygous ALT SNPs in sample %s</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">nhq</span><span class="p">,</span> <span class="n">argv</span><span class="p">[</span><span class="mi">2</span><span class="p">]);</span>
290        <span class="n">free</span><span class="p">(</span><span class="n">fi</span><span class="p">);</span>
291        <span class="n">free</span><span class="p">(</span><span class="n">gq</span><span class="p">);</span>
292        <span class="n">free</span><span class="p">(</span><span class="n">gt</span><span class="p">);</span>
293        <span class="n">free</span><span class="p">(</span><span class="n">dp</span><span class="p">);</span>
294        <span class="n">free</span><span class="p">(</span><span class="n">seqnames</span><span class="p">);</span>
295        <span class="n">bcf_hdr_destroy</span><span class="p">(</span><span class="n">hdr</span><span class="p">);</span>
296        <span class="n">bcf_close</span><span class="p">(</span><span class="n">inf</span><span class="p">);</span>
297        <span class="n">bcf_destroy</span><span class="p">(</span><span class="n">rec</span><span class="p">);</span>
298        <span class="k">return</span> <span class="n">EXIT_SUCCESS</span><span class="p">;</span>
299<span class="nl">error2:</span>
300        <span class="n">free</span><span class="p">(</span><span class="n">seqnames</span><span class="p">);</span>
301<span class="nl">error1:</span>
302        <span class="n">bcf_close</span><span class="p">(</span><span class="n">inf</span><span class="p">);</span>
303        <span class="n">bcf_hdr_destroy</span><span class="p">(</span><span class="n">hdr</span><span class="p">);</span>
304        <span class="k">return</span> <span class="n">EXIT_FAILURE</span><span class="p">;</span>
305<span class="p">}</span>
306
307</code></pre>
308</div>
309
310  </article>
311  <footer>
312    
313      <span class="disabled" href="">&laquo; Newer post</span>
314    
315       | <a href="#page-top">top</a> | 
316     
317      <a href="/2014/07/19/model-protein-dna-dwell-time-part1.html">Older post &raquo;</a>
318     
319  </footer>
320
321
322  
323</body>
324
325
326</html>
327
328

Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.