Merge pull request #52 from J35P312/find_sv_from_seq

Find sv from seq

Author: J35P312

setup.py

@@ -19,7 +19,7 @@
     ext_modules = []
 
 setup(name='svdb',
-      version='2.6.1',
+      version='2.6.2',
       url="https://github.com/J35P312/SVDB",
       author="Jesper Eisfeldt",
       author_email="jesper.eisfeldt@scilifelab.se",

svdb/__main__.py

@@ -36,7 +36,7 @@ def make_query_calls (args, queries, keyword):
         query_module.main(args)
 
 def main():
-    version = "2.6.1"
+    version = "2.6.2"
     parser = argparse.ArgumentParser(
         """SVDB-{}, use the build module to construct databases, use the query module to query the database usign vcf files, or use the hist module to generate histograms""".format(version), add_help=False)
     parser.add_argument('--build', help="create a db",

svdb/merge_vcf_module.py

@@ -109,6 +109,7 @@ def print_header(vcf_list, vcf_dictionary, args, command_line):
         print(subheader[entry].strip())
     if not args.notag:
         print("##INFO=<ID=VARID,Number=1,Type=String,Description=\"The variant ID of merged samples\">")
+        print("##INFO=<ID=FOUNDBY,Number=1,Type=Integer,Description=\"The number of files containing the variant\">")
         print("##INFO=<ID=set,Number=1,Type=String,Description=\"Source VCF for the merged record in SVDB\">")
         print("##INFO=<ID=svdb_origin,Number=1,Type=String,Description=\"pipe separated list of the VCF for the merged record in SVDB\">")
 

svdb/merge_vcf_module_cython.py

@@ -331,6 +331,7 @@ def merge(variants, samples, sample_order, sample_print_order, priority_order, a
             if not args.notag:
                 set_tag = determine_set_tag(priority_order, files)
                 line[7] += ";set={}".format(set_tag)
+                line[7] += ";FOUNDBY={}".format( len(set(filters_tag.keys())) )
 
             #add chrom information of all merged variants
             callers=[]

svdb/readVCF.py

@@ -52,33 +52,43 @@ def readVCFLine(line):
     elif "]" not in variation[4] and "[" not in variation[4]:
 
         chrB = chrA
+        posB = posA
+
         if "END" in description:
             posB = int(description["END"])
-
         elif "SVLEN" in description:
             posB = posA + abs(int(description["SVLEN"]))
-        else:
-            posB = posA
+
         # sometimes the fermikit intra chromosomal events are inverted i.e the end pos is a lower position than the start pos
         if posB < posA:
             tmp = posB
             posB = posA
             posA = tmp
-        event_type = variation[4].strip("<").rstrip(">")
 
+        nucleotides=set(["A","T","C","G","N"])
+        #SVtype is given in alt column
         if "<" in variation[4] and ">" in variation[4]:
+            event_type = variation[4].strip("<").rstrip(">")
             if "DUP" in event_type:
                 event_type = "DUP"
-        else:
-            if "SVTYPE" in description:
-                event_type = description["SVTYPE"]
+
+        #SVtype present in INFO
+        elif "SVTYPE" in description:
+            event_type = description["SVTYPE"]
+
+        #no SVTYPE, actual sequence is written in alt
+        elif set(list(variation[4])).union(nucleotides) == nucleotides:
+            if len(variation[4]) > len(variation[3]):
+                event_type="INS"
+            else:
+                event_type="DEL"
+                posB=posA+len(variation[3])-1
+
         #treat the insertion as single points
         if "INS" in event_type:
             posA=int(variation[1])
             posB=int(variation[1])
 
-
-
     # if the variant is given as a breakpoint, it is stored as a precise variant in the db
     else:
         B = variation[4]