Shankar merge fixes

COPYING

@@ -1,4 +1,4 @@
-Copyright (c) 2009-2012 Illumina, Inc.
+Copyright (c) 2009-2016 Illumina, Inc.
 
 Permission is hereby granted, free of charge, to any person obtaining
 a copy of this software and associated documentation files (the

Changelog

@@ -1,8 +1,9 @@
-v0.17.0
-- add new check_reference tool
-- add tty detection for all tools requiring input on stdin
-v0.16.1
-- merge: Fix vcf header from FILT to FILTER, doc method limitations
+v0.16.x
+- Fixed merge_variants bug associated with haploid chromosomes
+- Fixed merge_variants bug around indels (samples moving past locus would get
+"./." genotype). Implemented rewind function
+- Fixed merge_variants bug to merge snv and indel alleles at the same pos - write one
+record instead of multiple
 v0.16
 - Add new remove_region tool (contributed by S. Ajay)
 - Add 'invert' option to extract_variants

Makefile

@@ -13,15 +13,15 @@ export BOOST_ROOT := $(REDIST_DIR)/boost/stage
 .PHONY: all build clean install redist test
 
 
-all: install
+all: install test
 
 build: redist
 	$(MAKE) -C $(SRC_DIR)
 
 redist:
 	$(MAKE) -C $(REDIST_DIR)
 
-install: test
+install: build
 	mkdir -p $(BIN_DIR) && $(MAKE) -C $(SRC_DIR) $@ 
 
 clean: srcclean

scratch/make_release_tarball.bash

@@ -49,7 +49,7 @@ git archive --prefix=$pname_root/ HEAD | tar -x -C $outdir
 rm -rf $pname/scratch
 
 # make version number substitutions:
-for f in README.md; do
+for f in README.txt; do
     sed "s/\${VERSION}/$gitversion/" < $f >| $pname/$f
 done
 

src/Makefile

@@ -35,7 +35,7 @@ export LDLIBS = -ltabix_and_faidx -lz -lm
 
 GVCFTOOLS_HH := gvcftools.hh
 TRIOPROGS := trio twins merge_variants
-BLOCKPROGS := break_blocks check_reference extract_variants gatk_to_gvcf get_called_regions set_haploid_region remove_region
+BLOCKPROGS := break_blocks extract_variants gatk_to_gvcf get_called_regions set_haploid_region remove_region
 PROGS = $(TRIOPROGS) $(BLOCKPROGS) 
 PROG_OBJS = $(PROGS:%=%.o)
 

src/break_blocks.cpp

@@ -42,7 +42,7 @@
 
 #include "boost/program_options.hpp"
 
-#include <unistd.h>
+//#include <ctime>
 
 #include <iostream>
 #include <string>
@@ -161,9 +161,7 @@ try_main(int argc,char* argv[]) {
         po_parse_fail=true;
     }
 
-    const bool isStdinTerminal(isatty(fileno(stdin)));
-
-    if ((argc<=1) || (vm.count("help")) || po_parse_fail || isStdinTerminal) {
+    if ((argc<=1) || (vm.count("help")) || po_parse_fail) {
         log_os << "\n" << progname << " converts non-reference blocks to individual positions in specified regions\n\n";
         log_os << "version: " << gvcftools_version() << "\n\n";
         log_os << "usage: " << progname << " [options] < (g)VCF > unblocked_(g)VCF\n\n";

src/extract_variants.cpp

@@ -34,10 +34,10 @@
 #include "VcfHeaderHandler.hh"
 #include "vcf_util.hh"
 
+
 #include "boost/program_options.hpp"
 
-#include <cstdio>
-#include <unistd.h>
+//#include <ctime>
 
 #include <iostream>
 #include <string>
@@ -159,9 +159,7 @@ try_main(int argc,char* argv[]) {
         po_parse_fail=true;
     }
 
-    const bool isStdinTerminal(isatty(fileno(stdin)));
-
-    if ((vm.count("help")) || po_parse_fail || isStdinTerminal) {
+    if ((vm.count("help")) || po_parse_fail) {
         log_os << "\n" << progname << " extracts variants from a VCF file\n\n";
         log_os << "version: " << gvcftools_version() << "\n\n";
         log_os << "usage: " << progname << " [options] < (g)VCF > variants_only_VCF\n\n";

src/gatk_to_gvcf.cpp

@@ -24,6 +24,7 @@
 //
 //
 
+/// \file
 ///
 /// \author Chris Saunders
 ///

src/get_called_regions.cpp

@@ -33,10 +33,10 @@
 #include "VcfHeaderHandler.hh"
 #include "vcf_util.hh"
 
+
 #include "boost/program_options.hpp"
 
-#include <cstdio>
-#include <unistd.h>
+//#include <ctime>
 
 #include <iostream>
 #include <string>
@@ -219,9 +219,7 @@ try_main(int argc,char* argv[]) {
         po_parse_fail=true;
     }
 
-    const bool isStdinTerminal(isatty(fileno(stdin)));
-
-    if ((vm.count("help")) || po_parse_fail || isStdinTerminal) {
+    if ((vm.count("help")) || po_parse_fail) {
         log_os << "\n" << progname << " creates a bed file of called regions from a gVCF\n\n";
         log_os << "version: " << gvcftools_version() << "\n\n";
         log_os << "usage: " << progname << " [options] < gVCF > called.bed\n\n";

src/libblock/GatkVcfRecord.hh

@@ -46,11 +46,9 @@
 
 
 ///
-/// A value which might exist, and might be convertible to an integer
-/// if it does.
-///
-/// 'Convertible to an integer' in this case are ints or floats which
-/// are rounded to the nearest int
+/// A value which might exist, and might be convertable to an integer
+/// if it does 'convertable to an integer' in this case are ints or
+/// floats which are rounded to the nearest int
 ///
 struct MaybeInt {
 

src/libtrio/related_sample_util.cpp

@@ -1,6 +1,6 @@
 // -*- mode: c++; indent-tabs-mode: nil; -*-
 //
-// Copyright (c) 2009-2012 Illumina, Inc.
+// Copyright (c) 2009-2015 Illumina, Inc.
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -26,7 +26,7 @@
 
 /// \file
 ///
-/// \author Chris Saunders
+/// \author Chris Saunders and Subramanian Shankar Ajay
 ///
 
 
@@ -118,7 +118,8 @@ get_indel_allele(
         }
     }
 
-    return (is_standard_diploid && _gtcode.size()==2);
+    //return (is_standard_diploid && _gtcode.size()==2);
+    return ((is_standard_diploid && _gtcode.size() == 2) || (_gtcode.size() == 1 && _gtcode[0] >= 0));
 }
 
 
@@ -159,10 +160,24 @@ get_site_allele(
         }
     }
 
-    return (is_standard_diploid && _gtcode.size()==2);
+    //return (is_standard_diploid && _gtcode.size()==2);
+    return ((is_standard_diploid && _gtcode.size() == 2) || (_gtcode.size() == 1 && _gtcode[0] >= 0));
 }
 
 
+bool
+snp_type_info::is_ref() const{
+
+    BOOST_FOREACH(const int gt, _gtcode)
+    {
+        if(gt == 0)
+            continue;
+        else
+            return false;
+    }
+    return true;
+}                       
+
 
 // extract unsigned value for key from the vcf info field
 bool
@@ -256,6 +271,7 @@ site_crawler(const sample_info& si,
     , _tabs(NULL)
     , _is_sample_begin_state(true)
     , _is_sample_end_state(false)
+    , _is_reference(false)
     , _next_file(0)
     , _ref_seg(ref_seg)
     , _n_word(0)
@@ -266,6 +282,9 @@ site_crawler(const sample_info& si,
     , _is_site_allele_current(false)
     , _is_indel_allele_current(false)
 {
+    enum {MAX_WORD=50};
+    memset(_prev_word,0,sizeof(char *)*(MAX_WORD));
+    memset(_next_word,0,sizeof(char *)*(MAX_WORD));
     update(is_store_header);
 }
 
@@ -274,19 +293,27 @@ site_crawler(const sample_info& si,
 site_crawler::
 ~site_crawler() {
     if (NULL != _tabs) delete _tabs;
+    if (NULL != _prev_word[0]){
+        unsigned idx(0);
+        while (idx < _n_word)
+            free(_prev_word[idx++]);
+    }
 }
 
 
 
 
-
 void
 site_crawler::
 dump_state(std::ostream& os) const {
     const std::string& afile(_si.file);
     os << "LOCUS_CRAWLER STATE:\n";
     os << "\tchrom: " << chrom();
     os << "\tposition: " << pos() << " offset: " << _locus_offset << "\n";
+    os << "\t_is_call: " << _is_call << "\n";
+    os << "\t_is_site_allele_current: " << _is_site_allele_current;
+    os << "\t_is_sample_begin_state: " << _is_sample_begin_state;
+    os << "\t_is_sample_end_state: " << _is_sample_end_state;
     os << "\tis_indel: " << is_indel() << "\n";
     os << "\tfile: '" << afile << "'\n";
     os << "\tline: '";
@@ -296,6 +323,8 @@ dump_state(std::ostream& os) const {
 
 
 
+
+
 bool
 site_crawler::
 update_site_allele() const
@@ -437,7 +466,8 @@ process_record_line(char* line)
     //}
 
     if (! _is_sample_begin_state) {
-        if (! (last_vpos < vpos()) ) {
+        //if (! (last_vpos < vpos()) ) {
+        if (vpos() < last_vpos ) {
             if (_opt.is_murdock_mode) {
                 _vpos=last_vpos;
                 _locus_size=0;
@@ -474,6 +504,46 @@ dump_header(std::ostream& os) const {
 }
 
 
+bool
+site_crawler::
+rewind_site(const vcf_pos& lpos) {
+
+    if(_locus_offset > 0){
+        while(_vpos.pos > lpos.pos){
+            _locus_offset--;
+            _vpos.pos--;
+        }
+        _is_call = update_allele();
+        _is_reference = _opt.sti().is_ref();
+    }
+    else{
+        unsigned idx(0);
+        while(idx < _n_word){
+            _next_word[idx] = _word[idx];
+            _word[idx] = _prev_word[idx];
+            idx++;
+        }
+        _vpos.pos = _opt.sti().pos(_word);
+        _is_site_allele_current = false;
+        _is_indel_allele_current = false;
+        _vpos.is_indel = (_opt.sti().get_is_indel(_word));
+        _opt.sti().get_nonindel_ref_length(pos(),is_indel(),_word,_locus_size);
+
+        if(is_indel())
+            _locus_size = 1;
+
+        _locus_offset = _locus_size - 1;
+        _vpos.pos += _locus_offset;
+        while(_vpos.pos > lpos.pos){
+            _locus_offset--;
+            _vpos.pos--;
+        }
+        _is_call = update_allele();
+        _is_reference = _opt.sti().is_ref();
+    }
+
+    return true;
+}
 
 
 void
@@ -502,6 +572,7 @@ update(bool is_store_header) {
             }
 
             _vpos.pos++;
+            _is_site_allele_current = false;
 
             if (_opt.is_region()) {
                 // check for pos preceding the start of region of interest in a multi-base record:
@@ -564,6 +635,38 @@ update(bool is_store_header) {
             _next_file++;
         }
 
+        if (_next_word[0]){
+            unsigned idx(0);
+            while(idx < _n_word){
+                _word[idx] = _next_word[idx];
+                _next_word[idx] = NULL;
+                idx++;
+            }
+            _vpos.pos = _opt.sti().pos(_word);
+            _locus_offset = 0;
+            _is_site_allele_current = false;
+            _is_indel_allele_current = false;
+            _vpos.is_indel=(_opt.sti().get_is_indel(_word));
+            _opt.sti().get_nonindel_ref_length(pos(),is_indel(),_word,_locus_size);
+            _is_call=update_allele();
+
+            return;
+        }
+
+        unsigned idx(0);
+        if ( ! _is_sample_begin_state){
+            while (idx < _n_word){
+                size_t size = strlen(_word[idx]) + sizeof(char) ;
+
+                if (NULL != _prev_word[idx])
+                    free(_prev_word[idx]);
+                _prev_word[idx] = (char *)malloc(size);
+                strncpy(_prev_word[idx], _word[idx], size);
+                idx++;
+            }
+        }
+
+
         // read through file to get to a data line:
         bool is_eof(true);
         char* line(NULL);