Merge pull request #58 from miRTop/dev

release 0.4.22

HISTORY.md

@@ -1,3 +1,8 @@
+- 0.4.22
+
+* fix when reads map halfway on to the edge
+* fix edge case where limit==variant_size 
+
 - 0.4.21
 
 * Missing trimming events since 0.4.19

mirtop/bam/bam.py

@@ -110,6 +110,7 @@ def low_memory_genomic_bam(bam_fn, sample, out_handle, args):
     logger.info("Intersecting bed file.")
     intersect_fn = intersect(bed_fn, args.gtf)
     logger.info("Loading database.")
+    # TODO this'll return conn_reads and conn_counts
     conn = _read_lifted_bam_alpha(intersect_fn, bam_fn, args)
     rows = sql.select_all_reads(conn)
     lines = []
@@ -120,6 +121,8 @@ def low_memory_genomic_bam(bam_fn, sample, out_handle, args):
             lines.append(row)
             current = row[0]
         else:
+            # TODO counts of sequence = conn_counts.query UID
+            # it could be counts only same location UID+chrom+start, or counts all UID
             reads = _read_lifted_lines(lines, precursors, database)
             ann = annotate(reads, args.matures, args.precursors, quiet=True)
             gff_lines = body.create(ann, args.database, sample, args, quiet=True)
@@ -152,23 +155,22 @@ def _analyze_line(line, reads, precursors, handle, args):
     if query_name not in reads:
         reads[query_name].set_sequence(sequence)
         reads[query_name].counts = _get_freq(query_name)
+        # TODO if args.quant set to 0
+    # TODO if args.quant increase by 1
     if line.is_reverse and not args.genomic:
         logger.debug("READ::Sequence is reverse: %s" % line.query_name)
         return reads
     chrom = handle.getrname(line.reference_id)
     start = line.reference_start
-    # If genomic endcode, liftover to precursor position
- #   if not start:
- #       logger.debug(("READ::not start found %s" % line.reference_start))
- #       return reads
+
     cigar = line.cigartuples
     # if line.cigarstring.find("I") > -1:
     #     indels_skip += 1
     iso = isomir()
     iso.align = line
     iso.set_pos(start, len(reads[query_name].sequence))
     logger.debug("READ::From BAM start %s end %s at chrom %s" % (iso.start, iso.end, chrom))
-    if len(precursors[chrom]) < start + len(reads[query_name].sequence):
+    if len(precursors[chrom].replace("N","")) + 3 < start + len(reads[query_name].sequence):
         logger.debug("READ::%s start + %s sequence size are bigger than"
                      " size precursor %s" % (
                                              line.reference_id,
@@ -198,11 +200,13 @@ def _read_lifted_bam_alpha(bed_fn, bam_fn, args):
     conn = sql.create_connection()
     key = "name" if args.keep_name else "sequence"
     sql.create_reads_table(conn, key)
+    # TODO create counts table sequence and autoincrement or from read
     cur = conn.cursor()
     counts = 0
     seen = set()
     for line in bed_fn:
         fields = _parse_intersect(line, database, bed=True)
+        # TODO add sequence to count table args.quant on/off name=UID or name=UID+chrom+pos
         if fields:
             hit = ".".join(fields[:3])
             if hit not in seen:
@@ -211,9 +215,10 @@ def _read_lifted_bam_alpha(bed_fn, bam_fn, args):
                 seen.add(hit)
         # if counts == 1000:
         #     counts = 0
-    del(hit)
+    del(seen)
     logger.info("Read %s lines that intersected with miRNAs." % counts)
     conn.commit()
+    # TODO this'll return conn_reads and conn_counts
     return conn
 
 

mirtop/command_line.py

@@ -24,7 +24,7 @@ def main(**kwargs):
                       kwargs['args'].print_debug)
     logger = mylog.getLogger(__name__)
     start = time.time()
-
+    logger.warning("This is devel-live changes")
     if "gff" in kwargs:
         logger.info("Run annotation")
         reader(kwargs["args"])

mirtop/importer/manatee.py

@@ -139,7 +139,7 @@ def _bed(handle, bed_fn):
             query_sequence = cols[9]
             counts = cols[14]
             start = int(cols[3])
-            strand = cols[1]
+            strand = "-" if cols[1] == "16" else "+"
             chrom = cols[2]
             # if there no hits
             # if the sequence always matching the read, assuming YES now

mirtop/mirna/realign.py

@@ -464,13 +464,14 @@ def align_from_variants(sequence, mature, variants):
     if "iso_add3p" in k:
         sequence = sequence[:-1 * var_dict["iso_add3p"]]
     if "iso_3p" in k:
-        mature = mature[:-(7 + (-1 * var_dict["iso_3p"]))]
+        shift = 7 + (-1 * var_dict["iso_3p"])
+        if shift != 0:
+            mature = mature[:-(shift)]
     else:
         mature = mature[:-7]
     logger.debug("realign::align_from_variants::snp %s" % snp)
     logger.debug("realign::align_from_variants::sequence %s" % sequence)
     logger.debug("realign::align_from_variants::mature %s" % mature)
-
     if len(sequence) != len(mature):  # in case of indels, align again
         a = align(sequence, mature)
         sequence = a[0]

setup.py

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