various changes to simplify overlap matching with features of interest

Makefile

@@ -54,7 +54,9 @@ SOURCES = program\
 		estimator\
 		profile_afs\
 		profile_hwe\
-		profile_len
+		profile_len\
+		overlap_region_matcher\
+		bed
 
 SOURCESONLY = main.cpp
 

annotate_variants.cpp

@@ -40,14 +40,15 @@ class Igor : Program
     std::string output_vcf_file;
     std::vector<GenomeInterval> intervals;
     std::string interval_list;
-    std::string gencode_gtf_file;    
+    std::string gencode_gtf_file;
     bool annotate_coding;
 
     ///////
     //i/o//
     ///////
     BCFOrderedReader *odr;
     BCFOrderedWriter *odw;
+    OverlapRegionMatcher *orm;
 
     /////////
     //stats//
@@ -114,19 +115,20 @@ class Igor : Program
 //        bcf_hdr_append(odw->hdr, "##INFO=<ID=LFLANK,Number=1,Type=String,Description=\"Right Flank\">");
 //        bcf_hdr_append(odw->hdr, "##INFO=<ID=RFLANK,Number=1,Type=String,Description=\"Left Flank\">");
 //        bcf_hdr_append(odw->hdr, "##INFO=<ID=NS,Number=0,Type=Flag,Description=\"Near to STR\">");
+        bcf_hdr_append(odw->hdr, "##INFO=<ID=LOW_COMPLEXITY,Number=0,Type=Flag,Description=\"INDEL is in low complexity region\">");
 
         ///////////////////////
         //tool initialization//
         ///////////////////////
         vm = new VariantManip(ref_fasta_file);
-        
+
         if (annotate_coding)
-        {    
+        {
             bcf_hdr_append(odw->hdr, "##INFO=<ID=GENCODE_FS,Number=0,Type=Flag,Description=\"Frameshift INDEL\">");
             bcf_hdr_append(odw->hdr, "##INFO=<ID=GENCODE_NFS,Number=0,Type=Flag,Description=\"Non Frameshift INDEL\">");
             gc = new GENCODE(gencode_gtf_file, ref_fasta_file);
         }
-        
+
         ////////////////////////
         //stats initialization//
         ////////////////////////
@@ -155,6 +157,9 @@ class Igor : Program
     {
         odw->write_hdr();
 
+        std::string str = "/net/fantasia/home/atks/ref/vt/grch37/mdust.bed.gz";
+        orm = new OverlapRegionMatcher(str);
+
         bcf1_t *v = odw->get_bcf1_from_pool();
         std::vector<Interval*> overlaps;
         Variant variant;
@@ -166,40 +171,33 @@ class Igor : Program
             std::string chrom = bcf_get_chrom(odr->hdr,v);
             int32_t start1 = bcf_get_pos1(v);
             int32_t end1 = bcf_get_end_pos1(v);
-            
+
             vm->vtype2string(vtype, &s);
             if (s.l)
-            {    
+            {
                 bcf_update_info_string(odr->hdr, v, "VT", s.s);
             }
-            
+
             if (vtype==VT_SNP)
             {
                 //synonymous and non synonymous annotation
-                
-            }    
+
+            }
             else if (vtype&VT_INDEL)
             {
+                if (orm->overlaps_with(chrom, start1+1, end1))
+                {
+                    bcf_update_info_flag(odr->hdr, v, "LOW_COMPLEXITY", "", 1);
+                }
+
                 //frame shift annotation
                 if (annotate_coding)
                 {
-//                    if (gc->overlaps_with(chrom, start1+1, end1, GC_FT_CDS))
-//                    {
-//                        if (variant.exists_frame_shift())
-//                        {
-//                            bcf_update_info_flag(odr->hdr, v, "GENCODE_FS", "", 1);
-//                        }
-//                        else
-//                        {
-//                            bcf_update_info_flag(odr->hdr, v, "GENCODE_NFS", "", 1);
-//                        }
-//                    }    
-
                     gc->search(chrom, start1+1, end1, overlaps);
-    
+
                     bool cds_found = false;
                     bool is_fs = false;
-    
+
                     for (int32_t i=0; i<overlaps.size(); ++i)
                     {
                         GENCODERecord *rec = (GENCODERecord *) overlaps[i];
@@ -213,7 +211,7 @@ class Igor : Program
                             }
                         }
                     }
-                    
+
                     if (cds_found)
                     {
                         if (is_fs)
@@ -225,23 +223,23 @@ class Igor : Program
                             bcf_update_info_flag(odr->hdr, v, "GENCODE_NFS", "", 1);
                         }
                     }
-                    
-                    //classify STR 
+
+                    //classify STR
                     std::string ru = "ACGT";
                     int32_t rl = 4;
                 }
     //            bcf_update_info_string(odr->hdr, v, "RU", ru.c_str());
-    //            bcf_update_info_int32(odr->hdr, v, "RL", &rl, 1); 
+    //            bcf_update_info_int32(odr->hdr, v, "RL", &rl, 1);
             }
-            
+
             ++no_variants_annotated;
             odw->write(v);
             v = odw->get_bcf1_from_pool();
         }
-        
+
         odw->close();
     };
-    
+
     private:
 
 };

annotate_variants.h

@@ -48,6 +48,7 @@
 #include "variant_manip.h"
 #include "log_tool.h"
 #include "gencode.h"
+#include "overlap_region_matcher.h"
 
 void annotate_variants(int argc, char ** argv);
 

bed.cpp

@@ -0,0 +1,24 @@
+/* The MIT License
+
+Copyright (c) 2014 Adrian Tan <[email protected]>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+*/
+
+#include "bed.h"

bed.h

@@ -0,0 +1,65 @@
+/* The MIT License
+
+   Copyright (c) 2014 Adrian Tan <[email protected]>
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+   copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in
+   all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+   THE SOFTWARE.
+*/
+
+#ifndef BED_H
+#define BED_H
+
+#include "htslib/kstring.h"
+#include "utils.h"
+#include "interval_tree.h"
+
+class BEDRecord : public Interval
+{
+    public:
+    std::string chrom;
+    int32_t start1, end1;
+
+    BEDRecord(kstring_t *s)
+    {
+        std::vector<std::string> fields;
+        split(fields, "\t", s->s);
+
+        chrom = fields[0];
+        str2int32(fields[3], start1);
+        str2int32(fields[4], end1);
+    };
+
+    BEDRecord(std::string& chrom, int32_t start1, int32_t end1)
+    {
+        this->chrom = chrom;
+        this->start1 = start1;
+        this->end1 = end1;
+    };
+
+    /**
+     * Prints this BED record to STDERR.
+     */
+    void print()
+    {
+        std::cerr << this->chrom << ":" << this->start1 << "-" <<this->end1 << "\n";
+    };
+
+    private:
+};
+
+#endif

genome_interval.cpp

@@ -126,4 +126,12 @@ void GenomeInterval::to_string(kstring_t *interval)
         kputc('-', interval);
         kputw(end1, interval);
     }
+};
+
+/**
+ * Checks if this interval overlap.
+ */
+bool GenomeInterval::overlaps_with(std::string& chrom, int32_t start1, int32_t end1)
+{
+    return (seq==chrom && this->start1<=end1 && this->end1>=start1);
 };

genome_interval.h

@@ -72,7 +72,7 @@ class GenomeInterval
      * Sets an interval.
      */
     void set(std::string interval);
-        
+
     /**
      * Converts genome interval into the entire chromosome.
      */
@@ -87,6 +87,11 @@ class GenomeInterval
      * Returns a string representation of this Genome Interval.
      */
     void to_string(kstring_t *interval);
+
+    /**
+     * Checks if this interval overlap.
+     */
+    bool overlaps_with(std::string& chrom, int32_t start1, int32_t end1);
 };
 
 #endif

overlap_region_matcher.cpp

@@ -0,0 +1,25 @@
+/* The MIT License
+
+   Copyright (c) 2014 Adrian Tan <[email protected]>
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+   copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in
+   all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+   THE SOFTWARE.
+*/
+
+#include "overlap_region_matcher.h"
+