Dev (#794)

* Add --human_chr to gff.format()

* update gtf attributes

* fix tests and deprecation warnings

Author: tanghaibao

src/jcvi/apps/phylo.py

@@ -481,7 +481,7 @@ def SH_raxml(reftree, querytree, phy_file, shout="SH_out.txt"):
     except Exception:
         print("SH test failed.", file=sys.stderr)
     else:
-        pval = pval.strip().replace("\t", " ").replace("%", "\%")
+        pval = pval.strip().replace("\t", " ").replace("%", r"\%")
         print("{0}\t{1}".format(op.basename(querytree), pval), file=shout)
         logger.debug("SH p-value appended to %s" % shout.name)
 

src/jcvi/formats/gff.py

@@ -1,11 +1,12 @@
 #!/usr/bin/env python
 # -*- coding: UTF-8 -*-
 
-from collections import defaultdict
 import os
 import os.path as op
 import re
 import sys
+
+from collections import defaultdict
 from urllib.parse import quote, unquote
 
 from ..annotation.reformat import atg_name
@@ -67,7 +68,12 @@
 )
 multiple_gff_attributes = ("Parent", "Alias", "Dbxref", "Ontology_term")
 safechars = " /:?~#+!$'@()*[]|"
-VALID_HUMAN_CHROMOSMES = set([str(x) for x in range(1, 23)] + ["X", "Y"])
+VALID_HUMAN_CHROMOSOMES = {
+    **{str(i): f"chr{i}" for i in range(1, 23)},
+    "X": "chrX",
+    "Y": "chrY",
+    "MT": "chrM",
+}
 
 
 class GffLine(object):
@@ -209,6 +215,7 @@ def update_attributes(self, skipEmpty=True, gff3=True, gtf=None, urlquote=True):
         attributes = []
         if gtf:
             gff3 = None
+            urlquote = False
         elif gff3 is None:
             gff3 = self.gff3
 
@@ -1495,7 +1502,6 @@ def chain(args):
             prev_gid = curr_gid
 
     for gkey, v in sorted(gffdict.items()):
-        gseqid, key = gkey
         seqid = v["seqid"]
         source = v["source"]
         type = v["type"]
@@ -1672,7 +1678,12 @@ def format(args):
         action="store_true",
         help="Store entire GFF file in memory during first iteration",
     )
-
+    p.add_argument(
+        "--human_chr",
+        default=False,
+        action="store_true",
+        help="Only allow 1-22XY/MT, and add `chr` prefix to seqid",
+    )
     p.set_outfile()
     p.set_SO_opts()
 
@@ -1843,6 +1854,10 @@ def format(args):
         so, _ = GODag_from_SO()
         valid_soterm = {}
 
+    gtf = ".gtf" in gffile and not opts.gff3
+    if gtf:
+        logger.info("Output in GTF format")
+
     fw = must_open(outfile, "w")
     if not make_gff_store:
         gff = Gff(gffile, keep_attr_order=(not opts.no_keep_attr_order), strict=strict)
@@ -2021,6 +2036,11 @@ def format(args):
                     _parent = [parent, "{0}-Protein".format(parent)]
                     g.set_attr("Parent", _parent)
 
+        if opts.human_chr:
+            if g.seqid not in VALID_HUMAN_CHROMOSOMES:
+                continue
+            g.seqid = VALID_HUMAN_CHROMOSOMES[g.seqid]
+
         pp = g.get_attr("Parent", first=False)
         if (
             opts.multiparents == "split" and (pp and len(pp) > 1) and g.type != "CDS"
@@ -2034,9 +2054,8 @@ def format(args):
                     fix_gsac(g, notes)
                 print(g, file=fw)
         else:
-            if g.gff3 and not opts.gff3:
-                opts.gff3 = True
-            g.update_attributes(gff3=opts.gff3)
+            gff3 = g.gff3 or opts.gff3
+            g.update_attributes(gff3=gff3, gtf=gtf)
             if gsac:
                 fix_gsac(g, notes)
             if duptype == g.type and skip[(g.seqid, g.idx, id, g.start, g.end)] == 1:
@@ -2576,7 +2595,6 @@ def gtf(args):
     gff = Gff(gffile, strict=not opts.nostrict)
     transcript_info = AutoVivification()
     for g in gff:
-        transcript_id = None
         if g.type.endswith(("RNA", "transcript")):
             if "ID" in g.attributes and "Parent" in g.attributes:
                 transcript_id = g.get_attr("ID")
@@ -3024,7 +3042,7 @@ def bed(args):
         "--human_chr",
         default=False,
         action="store_true",
-        help="Only allow 1-22XY, and add `chr` prefix to seqid",
+        help="Only allow 1-22XY/MT, and add `chr` prefix to seqid",
     )
     p.add_argument(
         "--ensembl_cds",
@@ -3090,9 +3108,9 @@ def bed(args):
         if span:
             bl.score = bl.span
         if human_chr:
-            if bl.seqid not in VALID_HUMAN_CHROMOSMES:
+            if bl.seqid not in VALID_HUMAN_CHROMOSOMES:
                 continue
-            bl.seqid = "chr" + bl.seqid
+            bl.seqid = VALID_HUMAN_CHROMOSOMES[bl.seqid]
         if ensembl_cds:
             if g.get_attr("gene_biotype") != "protein_coding":
                 continue

src/jcvi/formats/sam.py

@@ -195,9 +195,9 @@ def get_minibam_bed(bamfile, bedfile, minibam=None):
 
     cmd = "cat {}".format(bedfile)
     cmd += " | perl -lane 'print \"$F[0]:$F[1]-$F[2]\"'"
-    cmd += " | xargs -n1 -t -I \{\}"
+    cmd += r" | xargs -n1 -t -I \{\}"
     cmd += " samtools view {}".format(bamfile)
-    cmd += " \{\} >> " + minisamfile
+    cmd += r" \{\} >> " + minisamfile
     sh(cmd)
 
     cmd = "samtools view {} -b".format(minisamfile)

src/jcvi/projects/allmaps.py

@@ -256,7 +256,7 @@ def estimategaps(args):
     root.text(
         sum(markers) / 2,
         ypos + pad,
-        "Distance: 1.29cM $\Leftrightarrow$ 211,824bp (6.1 cM/Mb)",
+        r"Distance: 1.29cM $\Leftrightarrow$ 211,824bp (6.1 cM/Mb)",
         **fontprop,
     )
 
@@ -310,7 +310,7 @@ def lms(args):
     xdata = [x + randint(-3, 3) for x in range(10, 110, 10)]
     ydata = [x + randint(-3, 3) for x in range(10, 110, 10)]
     ydata[3:7] = ydata[3:7][::-1]
-    xydata = zip(xdata, ydata)
+    xydata = list(zip(xdata, ydata))
     lis = xydata[:3] + [xydata[4]] + xydata[7:]
     lds = xydata[3:7]
     xlis, ylis = zip(*lis)
@@ -425,7 +425,7 @@ def subplot_twinx(
     legend=None,
     loc="upper left",
 ):
-    columned_data = zip(*data)
+    columned_data = list(zip(*data))
     x, yy = columned_data[0], columned_data[1:]
     assert len(ylabels) == 2
     assert len(yy) == 2
@@ -464,7 +464,7 @@ def subplot_twinx(
 def subplot(
     ax, data, xlabel, ylabel, xlim=None, ylim=1.1, xcast=float, ycast=float, legend=None
 ):
-    columned_data = zip(*data)
+    columned_data = list(zip(*data))
     x, yy = columned_data[0], columned_data[1:]
     lines = []
     for y, m in zip(yy, "o^x"):

src/jcvi/projects/misc.py

@@ -99,8 +99,8 @@ def grabseeds(args):
 
     fig = plt.figure(1, (iopts.w, iopts.h))
     ax = fig.add_subplot(1, 1, 1)
-    ax.set_xlabel(f"Principal Component 1 ({pc1_var * 100:.0f}\%)", fontsize=15)
-    ax.set_ylabel(f"Principal Component 2 ({pc2_var * 100:.0f}\%)", fontsize=15)
+    ax.set_xlabel(f"Principal Component 1 ({pc1_var * 100:.0f}%)", fontsize=15)
+    ax.set_ylabel(f"Principal Component 2 ({pc2_var * 100:.0f}%)", fontsize=15)
     ax.set_title("Sorghum kernels, PCA Plot", fontsize=20)
     ax.scatter("PC1", "PC2", s="ScatterSize", c="Color", data=final_df)
     set_helvetica_axis(ax)

src/jcvi/utils/db.py

@@ -35,7 +35,7 @@ def db_defaults(connector="Sybase"):
 def get_profile(
     sqshrc="~/.sqshrc", connector="Sybase", hostname=None, username=None, password=None
 ):
-    """
+    r"""
     get database, username, password from .sqshrc file e.g.
     \set username="user"
     """

tests/formats/test_obo.py

@@ -11,9 +11,9 @@ def test_oboreader():
     r4 = go["GO:0000010"]
     assert r1.item_id == "GO:0000001"
     assert r1.name == "mitochondrion inheritance"
-    assert r2.item_id == "GO:0000002"
-    assert r2.namespace == "biological_process"
-    assert r3.item_id == "GO:0000006"
+    assert r2.item_id == "GO:0000006"
+    assert r2.namespace == "molecular_function"
+    assert r3.item_id == "GO:0000007"
     assert tuple(sorted(r4.alt_ids)) == ("GO:0036422",)
 
     if os.path.exists(obo_file):