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 <stdio.h> //puts and printf 69#include <stdlib.h> //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">&</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"><</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 -> %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 > 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 <stdio.h> 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 <vcf file>, extract all calls for <sample> 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 & GQ > 20 & 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">&</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"><</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 -> %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">&</span><span class="n">fi</span><span class="p">,</span> <span class="o">&</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">&</span><span class="n">gq</span><span class="p">,</span> <span class="o">&</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">&</span><span class="n">dp</span><span class="p">,</span> <span class="o">&</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">&</span><span class="n">gt</span><span class="p">,</span> <span class="o">&</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">&&</span> <span class="n">gq</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">></span> <span class="mi">20</span> <span class="o">&&</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">&&</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">-></span><span class="n">rid</span><span class="p">],</span> 281 <span class="n">rec</span><span class="o">-></span><span class="n">pos</span><span class="p">,</span> 282 <span class="n">rec</span><span class="o">-></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">-></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="">« 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 »</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.