Merge pull request #45 from ramprasadn/feat/multiple_dbs

J35P312 · web-flow · commit bcff1c94ba00 · 2022-03-11T09:15:01.000+01:00
Nicely done!
diff --git a/README.md b/README.md
@@ -18,7 +18,7 @@ The GNOMAD SV database:
 
 https://storage.googleapis.com/gnomad-public/papers/2019-sv/gnomad_v2_sv.sites.vcf.gz
 
-external databses are run like this:
+external databases are run like this:
 
 ```bash
 svdb --query \
@@ -67,7 +67,7 @@ This module is used to construct structural variant databases from vcf files. Th
 
 
 ## Export
-This module is used to export the variants of the SVDB sqlite database. The variants of the sqlite svdb database is clustered using one out of three algorihms, overlap or DBSCAN.
+This module is used to export the variants of the SVDB sqlite database. The variants of the sqlite svdb database is clustered using one out of three algorithms, overlap or DBSCAN.
  
     print a help message
         svdb  --export --help  
@@ -81,9 +81,9 @@ This module is used to export the variants of the SVDB sqlite database. The vari
  
          --overlap OVERLAP     the overlap required to merge two events(0 means anything that touches will be merged, 1 means that two events must be identical to be merged), default = 0.8
 
-         --DBSCAN              use dbscan to cluster the variants, overides the overlap based clustering algoritm
+         --DBSCAN              use dbscan to cluster the variants, overides the overlap based clustering algorithm
 
-         --epsilon EPSILON     used together with --DBSCAN; sets the epsilon paramter(default = 500bp)
+         --epsilon EPSILON     used together with --DBSCAN; sets the epsilon parameter(default = 500bp)
 
          --min_pts MIN_PTS     the min_pts parameter(default = 2
 
@@ -92,30 +92,34 @@ This module is used to export the variants of the SVDB sqlite database. The vari
           --memory              load the database into memory: increases the memory requirements, but lowers the time consumption
 
 ## Query
-The query module is used to query a structural variant database. Typically a database is constructed using the build module. However, since this module utilize the genotype field of the sructural variant database vcf to compute the frequency of structural variants, a wide range of files could be used as database. The query module requires a query vcf, as well as a database file(either multisample vcf or SVDB sqlite database):
+The query module is used to query one or more structural variant databases. Typically a database is constructed using the build module. However, since this module utilize the genotype field of the structural variant database vcf to compute the frequency of structural variants, a wide range of files could be used as database. The query module requires a query vcf, as well as a database file(either multisample vcf or SVDB sqlite database):
 
     print a help message
        svdb --query --help
     Query a structural variant database, using a vcf file as query:
 
         svdb --query --query_vcf patient1.vcf --db control_db.vcf
 
+    Query multiple databases, using a vcf file as query:
+    
+        svdb --query --query_vcf patient1.vcf --db control_db1.vcf,control_db2.vcf --prefix test --in_occ default,Obs --in_frq FRQ,default --out_frq db1_AF,db2_Frq --out_occ db1_AC,db2_Obs
+
     optional arguments:
 
 	-h, --help            show this help message and exit
 
-	--db DB				path to a db vcf
-	--sqdb SQDB			path to a SVDB sqlite db
-	--bedpedb BEDPEDB		path to a SV database of the following format chrA-posA-chrB-posB-type-count-frequency
-	--in_occ IN_OCC			The allele count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AN or OCC)
-	--in_frq IN_FRQ			The frequency count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AF or FRQ)
-	--out_occ OUT_OCC		the allle count tag, as annotated by SVDB variant(defualt=OCC)
-	--out_frq OUT_FRQ		the tag used to describe the frequency of the variant(defualt=FRQ)
-	--prefix PREFIX			the prefix of the output file, default = print to stdout
+	--db DB				path to a db vcf, or a comma separated list of vcfs
+	--sqdb SQDB			path to a SVDB sqlite db, or a comma separated list of dbs
+	--bedpedb BEDPEDB		path to a SV database of the following format chrA-posA-chrB-posB-type-count-frequency, or a comma separated list of files
+	--in_occ IN_OCC			The allele count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AN or OCC). This parameter is required if multiple databases are queried. 
+	--in_frq IN_FRQ			The frequency count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AF or FRQ). This parameter is required if multiple databases are queried.
+	--out_occ OUT_OCC		the allele count tag, as annotated by SVDB variant(default=OCC). This parameter is required if multiple databases are queried.
+	--out_frq OUT_FRQ		the tag used to describe the frequency of the variant(default=FRQ). This parameter is required if multiple databases are queried.
+	--prefix PREFIX			the prefix of the output file, default = print to stdout. Required if multiple databases are queried. 
 	--bnd_distance BND_DISTANCE	the maximum distance between two similar breakpoints(default = 10000)
 	--overlap OVERLAP		the overlap required to merge two events(0 means anything that touches will be merged, 1 means that two events must be identical to be merged), default = 0.6
 	--memory 			load the database into memory: increases the memory requirements, but lowers the time consumption(may only be used with sqdb)
-	--no_var			count overlaping variants of different type as hits in the db
+	--no_var			count overlapping variants of different type as hits in the db
 
     
 ## Merge
diff --git a/setup.py b/setup.py
@@ -19,7 +19,7 @@
     ext_modules = []
 
 setup(name='svdb',
-      version='2.5.1',
+      version='2.5.2',
       url="https://github.com/J35P312/SVDB",
       author="Jesper Eisfeldt",
       author_email="jesper.eisfeldt@scilifelab.se",
diff --git a/svdb/__main__.py b/svdb/__main__.py
@@ -1,10 +1,42 @@
-import argparse
+import argparse, os
 
 from . import build_module, export_module, merge_vcf_module, query_module
 
+def make_query_calls (args, queries, keyword):
+    if len(queries) > 1 and args.prefix:
+        if all(variable is not None for variable in [args.in_occ, args.out_occ, args.in_frq, args.out_frq]):
+            in_occs     = args.in_occ.split(",")
+            in_frqs     = args.in_frq.split(",")
+            out_occs    = args.out_occ.split(",")
+            out_frqs    = args.out_frq.split(",")
+            orig_prefix = args.prefix
+            if (len(queries) == len(in_occs) == len(in_frqs) == len(out_occs) == len(out_frqs)):
+                for ind in range(len(queries)):
+                    if keyword   == "db":
+                        args.db      = queries[ind]
+                    elif keyword == "sqdb":
+                        args.sqdb    = queries[ind]
+                    elif keyword == "bedpedb":
+                        args.bedpedb = queries[ind]
+                    args.in_occ      = None if in_occs[ind] == "default" else in_occs[ind]
+                    args.in_frq      = None if in_frqs[ind] == "default" else in_frqs[ind]
+                    args.out_occ     = out_occs[ind]
+                    args.out_frq     = out_frqs[ind]
+                    if ind < len(queries)-1:
+                        args.prefix  = orig_prefix + "_" + str(ind)
+                    else:
+                        args.prefix = orig_prefix
+                    query_module.main(args)
+                    if ind > 0:
+                        os.remove(args.query_vcf)
+                    args.query_vcf = args.prefix + "_query.vcf"
+        else:
+            print("please ensure that both count and frequency tags are specified for all samples")
+    else:
+        query_module.main(args)
 
 def main():
-    version = "2.5.1"
+    version = "2.5.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",
@@ -22,22 +54,22 @@ def main():
             """SVDB.{}: query module""".format(version))
         parser.add_argument('--query', help="query a db", required=False, action="store_true")
         parser.add_argument('--query_vcf', type=str, help="a vcf used to query the db", required=True)
-        parser.add_argument('--db', type=str, help="path to a SVDB db vcf ")
-        parser.add_argument('--sqdb', type=str, help="path to a SVDB sqlite db")
+        parser.add_argument('--db', type=str, help="path to a SVDB db vcf or a comma separated list of vcfs")
+        parser.add_argument('--sqdb', type=str, help="path to a SVDB sqlite db or a comma separated list of dbs")
         parser.add_argument('--bedpedb', type=str,
-                            help="path to a SV database of the following format chrA-posA-chrB-posB-type-count-frequency")
+                            help="path to a SV database of the following format chrA-posA-chrB-posB-type-count-frequency, or a or a comma separated list of dbs")
         parser.add_argument('--in_occ', type=str,
-                            help="The allele count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AC or OCC)")
+                            help="The allele count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AC or OCC), required if multiple databases are queried. Use default (as shown in the example in README) if you'd like to use default tag for a specific database")
         parser.add_argument('--in_frq', type=str,
-                            help="The frequency count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AF or FRQ)")
+                            help="The frequency count tag, if used, this tag must be present in the INFO column of the input DB(usually set to AF or FRQ), required if multiple databases are queried. Use default (as shown in the example in README) if you'd like to use default tag for a specific database")
         parser.add_argument('--out_occ', type=str, default="OCC",
-                            help="the allle count tag, as annotated by SVDBvariant(defualt=OCC)")
+                            help="the allele count tag, as annotated by SVDBvariant(default=OCC), required if multiple databases are queried.")
         parser.add_argument('--out_frq', type=str, default="FRQ",
-                            help="the tag used to describe the frequency of the variant(defualt=FRQ)")
+                            help="the tag used to describe the frequency of the variant(default=FRQ), required if multiple databases are queried.")
         parser.add_argument('--max_frq', type=float, default=1,
                             help='Only include variants with a higher frequency than given here between 0 and 1. All new variants are always included. (default: 1)')
         parser.add_argument('--prefix', type=str, default=None,
-                            help="the prefix of the output file, default = print to stdout")
+                            help="the prefix of the output file, default = print to stdout. Required, if multiple databases are queried")
         parser.add_argument('--bnd_distance', type=int, default=10000,
                             help="the maximum distance between two similar breakpoints(default = 10000)")
         parser.add_argument('--ins_distance', type=int, default=50,
@@ -52,7 +84,15 @@ def main():
         args.version = version
 
         if(args.db or args.sqdb or args.bedpedb):
-            query_module.main(args)
+            if(args.db):
+                queries = args.db.split(",")
+                make_query_calls(args, queries, "db")
+            if(args.sqdb):
+                queries = args.sqdb.split(",")
+                make_query_calls(args, queries, "sqdb")
+            if(args.bedpedb):
+                queries = args.bedpedb.split(",")
+                make_query_calls(args, queries, "bedpedb")
         else:
             print("invalid db option, choose --db to use the vcf db or sqdb to use the sqlite db")
 
