Merge pull request #93 from miRTop/fix_mirgenedb_parsing_name_error

fix mirgenedb parser

Author: lpantano

data/examples/annotate/mirtop.db

@@ -4,4 +4,5 @@ dependencies:
   - bioconda::pybedtools
   - bioconda::samtools=1.21
   - conda-forge::pandas
+  - conda-forge::sqlite
   - conda-forge::biopython=1.83

environment.yml

@@ -28,7 +28,7 @@ def reader(args):
     args.database = database
     precursors = fasta.read_precursor(args.hairpin, args.sps)
     args.precursors = precursors
-    matures = mapper.read_gtf_to_precursor(args.gtf,database)
+    matures = mapper.read_gtf_to_precursor(args.gtf, database)
     args.matures = matures
     # TODO check numbers of miRNA and precursors read
     # TODO print message if numbers mismatch
@@ -78,13 +78,14 @@ def reader(args):
 def _write(lines, header, fn, args = None):
     out_handle = open(fn, 'w')
     print(header, file=out_handle)
-    mapper = read_gtf_to_mirna(args.gtf)
+    database = mapper.guess_database(args)
+    dbmapper = read_gtf_to_mirna(args.gtf, database)
     for m in lines:
         for s in sorted(lines[m].keys()):
             for hit in lines[m][s]:
                 # TODO: convert to genomic if args.out_genomic
                 if args and args.out_genomic:
-                    lifted = body.lift_to_genome(hit[4], mapper)
+                    lifted = body.lift_to_genome(hit[4], dbmapper)
                     print(lifted, file=out_handle)
                 else:
                     print(hit[4], file=out_handle)

mirtop/gff/__init__.py

@@ -1,6 +1,7 @@
 """Read database information"""
 
 from collections import defaultdict
+import re
 
 import mirtop.libs.logger as mylog
 
@@ -26,6 +27,8 @@ def guess_database(args):
 
 
 def _guess_database_file(gff, database=None):
+    if database:
+        return database
     with open(gff) as in_handle:
         for line in in_handle:
             if not line.startswith("#"):
@@ -56,7 +59,7 @@ def get_primary_transcript(database):
         raise ValueError("Only miRBase is supported for this action.")
 
 
-def read_gtf_to_mirna(gtf):
+def read_gtf_to_mirna(gtf, database=None):
     """
     Load GTF file with precursor positions on genome.
 
@@ -70,9 +73,11 @@ def read_gtf_to_mirna(gtf):
     """
     if not gtf:
         return gtf
-    if _guess_database_file(gtf).find("miRBase") > -1:
+    if not database:
+        database = _guess_database_file(gtf)
+    if database.find("miRBase") > -1:
         mapped = read_gtf_to_precursor_mirbase(gtf, format="genomic")
-    elif _guess_database_file(gtf).find("MirGeneDB") > -1:
+    elif database.find("MirGeneDB") > -1:
         mapped = read_gtf_to_precursor_mirgenedb(gtf, format="genomic")
     else:
         logger.info("Database different than miRBase or MirGeneDB")
@@ -145,7 +150,7 @@ def read_gtf_chr2mirna2(gtf):  # to remove
     return db_mir
 
 
-def read_gtf_to_precursor(gtf,database):
+def read_gtf_to_precursor(gtf, database=None):
     """
     Load GTF file with precursor positions on genome
     Return dict with key being precursor name and
@@ -163,9 +168,11 @@ def read_gtf_to_precursor(gtf,database):
     """
     if not gtf:
         return gtf
-    if _guess_database_file(gtf,database).find("miRBase") > -1:
+    if not database:
+        database = _guess_database_file(gtf)
+    if database.find("miRBase") > -1:
         mapped = read_gtf_to_precursor_mirbase(gtf)
-    elif _guess_database_file(gtf,database).find("MirGeneDB") > -1:
+    elif database.find("MirGeneDB") > -1:
         mapped = read_gtf_to_precursor_mirgenedb(gtf)
     else:
         logger.info("Database different than miRBase or MirGeneDB")
@@ -284,6 +291,7 @@ def read_gtf_to_precursor_mirgenedb(gtf, format="precursor"):
     db = defaultdict(list)
     db_mir = defaultdict(list)
     id_dict = dict()
+    pattern = r'(_3p\*?|_5p\*?)'
     with open(gtf) as in_handle:
         for line in in_handle:
             if line.startswith("#"):
@@ -299,7 +307,10 @@ def read_gtf_to_precursor_mirgenedb(gtf, format="precursor"):
             if cols[2] == "miRNA":
                 idname_mi = [n.split("=")[1] for n in cols[-1].split(";")
                              if n.startswith("ID")][0]
-                parent = "%s_pre" % idname_mi.split("_")[0]
+                # parent = "%s_pre" % idname_mi.replace("_3p.*", "").replace("_5p.*", "")
+                parent = re.sub(pattern, '', idname_mi)
+                parent = "%s_pre" % parent
+                # import pdb; pdb.set_trace()
                 db_mir[(parent, name)] = [chrom,
                                           int(start), int(end),
                                           strand, parent]

mirtop/mirna/mapper.py

@@ -3,7 +3,7 @@
 import os
 from setuptools import setup, find_packages
 
-version = '0.4.29'
+version = '0.4.30'
 url = 'http://github.com/mirtop/mirtop'