Repository 'cpt_intron_detection'
hg clone https://toolshed.g2.bx.psu.edu/repos/cpt/cpt_intron_detection

Changeset 0:1a19092729be (2022-05-13)
Next changeset 1:6ed12d7c2458 (2022-05-20)
Commit message:
Uploaded
added:
cpt_intron_detect/cpt-macros.xml
cpt_intron_detect/gff3.py
cpt_intron_detect/intron_detection.py
cpt_intron_detect/intron_detection.xml
cpt_intron_detect/macros.xml
cpt_intron_detect/test-data/T7_CLEAN.gff3
cpt_intron_detect/test-data/T7_INTRON.gff3
cpt_intron_detect/test-data/T7_NR.blastxml
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/cpt-macros.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/cpt-macros.xml Fri May 13 05:08:54 2022 +0000
[
@@ -0,0 +1,115 @@
+<?xml version="1.0"?>
+<macros>
+ <xml name="gff_requirements">
+ <requirements>
+ <requirement type="package" version="2.7">python</requirement>
+ <requirement type="package" version="1.65">biopython</requirement>
+ <requirement type="package" version="2.12.1">requests</requirement>
+ <yield/>
+ </requirements>
+ <version_command>
+ <![CDATA[
+ cd $__tool_directory__ && git rev-parse HEAD
+ ]]>
+ </version_command>
+ </xml>
+ <xml name="citation/mijalisrasche">
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+ <citation type="bibtex">@unpublished{galaxyTools,
+ author = {E. Mijalis, H. Rasche},
+ title = {CPT Galaxy Tools},
+ year = {2013-2017},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation>
+ </xml>
+ <xml name="citations">
+ <citations>
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+ <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {E. Mijalis, H. Rasche},
+ title = {CPT Galaxy Tools},
+ year = {2013-2017},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation> 
+ <yield/>
+ </citations>
+ </xml>
+     <xml name="citations-crr">
+ <citations>
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+ <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {C. Ross},
+ title = {CPT Galaxy Tools},
+ year = {2020-},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation>
+ <yield/>
+ </citations>
+ </xml>
+        <xml name="citations-2020">
+ <citations>
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+ <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {E. Mijalis, H. Rasche},
+ title = {CPT Galaxy Tools},
+ year = {2013-2017},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation>
+                        <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {A. Criscione},
+ title = {CPT Galaxy Tools},
+ year = {2019-2021},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+                        </citation>
+                        <yield/>
+ </citations>
+ </xml>
+        <xml name="citations-2020-AJC-solo">
+ <citations>
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+                        <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {A. Criscione},
+ title = {CPT Galaxy Tools},
+ year = {2019-2021},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+                        </citation>
+                        <yield/>
+ </citations>
+ </xml>
+        <xml name="citations-clm">
+ <citations>
+ <citation type="doi">10.1371/journal.pcbi.1008214</citation>
+ <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {C. Maughmer},
+ title = {CPT Galaxy Tools},
+ year = {2017-2020},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation>
+                        <yield/>
+ </citations>
+ </xml>
+        <xml name="sl-citations-clm">
+ <citation type="bibtex">
+ @unpublished{galaxyTools,
+ author = {C. Maughmer},
+ title = {CPT Galaxy Tools},
+ year = {2017-2020},
+ note = {https://github.com/tamu-cpt/galaxy-tools/}
+ }
+ </citation>
+                        <yield/>
+ </xml>
+</macros>
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/gff3.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/gff3.py Fri May 13 05:08:54 2022 +0000
[
b'@@ -0,0 +1,346 @@\n+import copy\n+import logging\n+\n+log = logging.getLogger()\n+log.setLevel(logging.WARN)\n+\n+\n+def feature_lambda(\n+    feature_list,\n+    test,\n+    test_kwargs,\n+    subfeatures=True,\n+    parent=None,\n+    invert=False,\n+    recurse=True,\n+):\n+    """Recursively search through features, testing each with a test function, yielding matches.\n+\n+    GFF3 is a hierachical data structure, so we need to be able to recursively\n+    search through features. E.g. if you\'re looking for a feature with\n+    ID=\'bob.42\', you can\'t just do a simple list comprehension with a test\n+    case. You don\'t know how deeply burried bob.42 will be in the feature tree. This is where feature_lambda steps in.\n+\n+    :type feature_list: list\n+    :param feature_list: an iterable of features\n+\n+    :type test: function reference\n+    :param test: a closure with the method signature (feature, **kwargs) where\n+                 the kwargs are those passed in the next argument. This\n+                 function should return True or False, True if the feature is\n+                 to be yielded as part of the main feature_lambda function, or\n+                 False if it is to be ignored. This function CAN mutate the\n+                 features passed to it (think "apply").\n+\n+    :type test_kwargs: dictionary\n+    :param test_kwargs: kwargs to pass to your closure when it is called.\n+\n+    :type subfeatures: boolean\n+    :param subfeatures: when a feature is matched, should just that feature be\n+                        yielded to the caller, or should the entire sub_feature\n+                        tree for that feature be included? subfeatures=True is\n+                        useful in cases such as searching for a gene feature,\n+                        and wanting to know what RBS/Shine_Dalgarno_sequences\n+                        are in the sub_feature tree (which can be accomplished\n+                        with two feature_lambda calls). subfeatures=False is\n+                        useful in cases when you want to process (and possibly\n+                        return) the entire feature tree, such as applying a\n+                        qualifier to every single feature.\n+\n+    :type invert: boolean\n+    :param invert: Negate/invert the result of the filter.\n+\n+    :rtype: yielded list\n+    :return: Yields a list of matching features.\n+    """\n+    # Either the top level set of [features] or the subfeature attribute\n+    for feature in feature_list:\n+        feature._parent = parent\n+        if not parent:\n+            # Set to self so we cannot go above root.\n+            feature._parent = feature\n+        test_result = test(feature, **test_kwargs)\n+        # if (not invert and test_result) or (invert and not test_result):\n+        if invert ^ test_result:\n+            if not subfeatures:\n+                feature_copy = copy.deepcopy(feature)\n+                feature_copy.sub_features = list()\n+                yield feature_copy\n+            else:\n+                yield feature\n+\n+        if recurse and hasattr(feature, "sub_features"):\n+            for x in feature_lambda(\n+                feature.sub_features,\n+                test,\n+                test_kwargs,\n+                subfeatures=subfeatures,\n+                parent=feature,\n+                invert=invert,\n+                recurse=recurse,\n+            ):\n+                yield x\n+\n+\n+def fetchParent(feature):\n+    if not hasattr(feature, "_parent") or feature._parent is None:\n+        return feature\n+    else:\n+        return fetchParent(feature._parent)\n+\n+\n+def feature_test_true(feature, **kwargs):\n+    return True\n+\n+\n+def feature_test_type(feature, **kwargs):\n+    if "type" in kwargs:\n+        return str(feature.type).upper() == str(kwargs["type"]).upper()\n+    elif "types" in kwargs:\n+      for x in kwargs["types"]:\n+        if str(feature.type).upper() == str(x).upper():\n+          return True\n+      return False\n+    raise Exception("Incorrect feature_test_type call, ne'..b'feature.location.start,\n+        # feature.location.end,\n+        # feature.location.strand\n+        # )\n+    return result\n+\n+\n+def get_gff3_id(gene):\n+    return gene.qualifiers.get("Name", [gene.id])[0]\n+\n+\n+def ensure_location_in_bounds(start=0, end=0, parent_length=0):\n+    # This prevents frameshift errors\n+    while start < 0:\n+        start += 3\n+    while end < 0:\n+        end += 3\n+    while start > parent_length:\n+        start -= 3\n+    while end > parent_length:\n+        end -= 3\n+    return (start, end)\n+\n+\n+def coding_genes(feature_list):\n+    for x in genes(feature_list):\n+        if (\n+            len(\n+                list(\n+                    feature_lambda(\n+                        x.sub_features,\n+                        feature_test_type,\n+                        {"type": "CDS"},\n+                        subfeatures=False,\n+                    )\n+                )\n+            )\n+            > 0\n+        ):\n+            yield x\n+\n+\n+def genes(feature_list, feature_type="gene", sort=False):\n+    """\n+    Simple filter to extract gene features from the feature set.\n+    """\n+\n+    if not sort:\n+        for x in feature_lambda(\n+            feature_list, feature_test_type, {"type": feature_type}, subfeatures=True\n+        ):\n+            yield x\n+    else:\n+        data = list(genes(feature_list, feature_type=feature_type, sort=False))\n+        data = sorted(data, key=lambda feature: feature.location.start)\n+        for x in data:\n+            yield x\n+\n+\n+def wa_unified_product_name(feature):\n+    """\n+    Try and figure out a name. We gave conflicting instructions, so\n+    this isn\'t as trivial as it should be. Sometimes it will be in\n+    \'product\' or \'Product\', othertimes in \'Name\'\n+    """\n+    # Manually applied tags.\n+    protein_product = feature.qualifiers.get(\n+        "product", feature.qualifiers.get("Product", [None])\n+    )[0]\n+\n+    # If neither of those are available ...\n+    if protein_product is None:\n+        # And there\'s a name...\n+        if "Name" in feature.qualifiers:\n+            if not is_uuid(feature.qualifiers["Name"][0]):\n+                protein_product = feature.qualifiers["Name"][0]\n+\n+    return protein_product\n+\n+\n+def is_uuid(name):\n+    return name.count("-") == 4 and len(name) == 36\n+\n+\n+def get_rbs_from(gene):\n+    # Normal RBS annotation types\n+    rbs_rbs = list(\n+        feature_lambda(\n+            gene.sub_features, feature_test_type, {"type": "RBS"}, subfeatures=False\n+        )\n+    )\n+    rbs_sds = list(\n+        feature_lambda(\n+            gene.sub_features,\n+            feature_test_type,\n+            {"type": "Shine_Dalgarno_sequence"},\n+            subfeatures=False,\n+        )\n+    )\n+    # Fraking apollo\n+    apollo_exons = list(\n+        feature_lambda(\n+            gene.sub_features, feature_test_type, {"type": "exon"}, subfeatures=False\n+        )\n+    )\n+    apollo_exons = [x for x in apollo_exons if len(x) < 10]\n+    # These are more NCBI\'s style\n+    regulatory_elements = list(\n+        feature_lambda(\n+            gene.sub_features,\n+            feature_test_type,\n+            {"type": "regulatory"},\n+            subfeatures=False,\n+        )\n+    )\n+    rbs_regulatory = list(\n+        feature_lambda(\n+            regulatory_elements,\n+            feature_test_quals,\n+            {"regulatory_class": ["ribosome_binding_site"]},\n+            subfeatures=False,\n+        )\n+    )\n+    # Here\'s hoping you find just one ;)\n+    return rbs_rbs + rbs_sds + rbs_regulatory + apollo_exons\n+\n+\n+def nice_name(record):\n+    """\n+    get the real name rather than NCBI IDs and so on. If fails, will return record.id\n+    """\n+    name = record.id\n+    likely_parental_contig = list(genes(record.features, feature_type="contig"))\n+    if len(likely_parental_contig) == 1:\n+        name = likely_parental_contig[0].qualifiers.get("organism", [name])[0]\n+    return name\n+\n+\n+def fsort(it):\n+    for i in sorted(it, key=lambda x: int(x.location.start)):\n+        yield i\n'
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/intron_detection.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/intron_detection.py Fri May 13 05:08:54 2022 +0000
[
b'@@ -0,0 +1,498 @@\n+#!/usr/bin/env python\n+import sys\n+import re\n+import itertools\n+import argparse\n+import hashlib\n+import copy\n+from CPT_GFFParser import gffParse, gffWrite, gffSeqFeature\n+from Bio.Blast import NCBIXML\n+from Bio.SeqFeature import SeqFeature, FeatureLocation\n+from gff3 import feature_lambda\n+from collections import OrderedDict\n+import logging\n+\n+logging.basicConfig(level=logging.DEBUG)\n+log = logging.getLogger()\n+\n+\n+def parse_xml(blastxml, thresh):\n+    """ Parses xml file to get desired info (genes, hits, etc) """\n+    blast = []\n+    discarded_records = 0\n+    totLen = 0\n+    for iter_num, blast_record in enumerate(NCBIXML.parse(blastxml), 1):\n+        blast_gene = []\n+        align_num = 0\n+        for alignment in blast_record.alignments:\n+            align_num += 1\n+            # hit_gis = alignment.hit_id + alignment.hit_def\n+            # gi_nos = [str(gi) for gi in re.findall(\'(?<=gi\\|)\\d{9,10}\', hit_gis)]\n+            gi_nos = str(alignment.accession)\n+\n+            for hsp in alignment.hsps:\n+                x = float(hsp.identities - 1) / ((hsp.query_end) - hsp.query_start)\n+                if x < thresh:\n+                    discarded_records += 1\n+                    continue\n+                nice_name = blast_record.query\n+\n+                if " " in nice_name:\n+                    nice_name = nice_name[0 : nice_name.index(" ")]\n+\n+                blast_gene.append(\n+                    {\n+                        "gi_nos": gi_nos,\n+                        "sbjct_length": alignment.length,\n+                        "query_length": blast_record.query_length,\n+                        "sbjct_range": (hsp.sbjct_start, hsp.sbjct_end),\n+                        "query_range": (hsp.query_start, hsp.query_end),\n+                        "name": nice_name,\n+                        "evalue": hsp.expect,\n+                        "identity": hsp.identities,\n+                        "identity_percent": x,\n+                        "hit_num": align_num,\n+                        "iter_num": iter_num,\n+                        "match_id": alignment.title.partition(">")[0],\n+                    }\n+                )\n+\n+        blast.append(blast_gene)\n+        totLen += len(blast_gene)\n+    log.debug("parse_blastxml %s -> %s", totLen + discarded_records, totLen)\n+    return blast\n+\n+\n+def filter_lone_clusters(clusters):\n+    """ Removes all clusters with only one member and those with no hits """\n+    filtered_clusters = {}\n+    for key in clusters:\n+        if len(clusters[key]) > 1 and len(key) > 0:\n+            filtered_clusters[key] = clusters[key]\n+    log.debug("filter_lone_clusters %s -> %s", len(clusters), len(filtered_clusters))\n+    return filtered_clusters\n+\n+\n+def test_true(feature, **kwargs):\n+    return True\n+\n+\n+def parse_gff(gff3):\n+    """ Extracts strand and start location to be used in cluster filtering """\n+    log.debug("parse_gff3")\n+    gff_info = {}\n+    _rec = None\n+    for rec in gffParse(gff3):\n+        endBase = len(rec.seq)\n+\n+        _rec = rec\n+        _rec.annotations = {}\n+        for feat in feature_lambda(rec.features, test_true, {}, subfeatures=False):\n+            if feat.type == "CDS":\n+                if "Name" in feat.qualifiers.keys():\n+                    CDSname = feat.qualifiers["Name"]\n+                else:\n+                    CDSname = feat.qualifiers["ID"]\n+                gff_info[feat.id] = {\n+                    "strand": feat.strand,\n+                    "start": feat.location.start,\n+                    "end": feat.location.end,\n+                    "loc": feat.location,\n+                    "feat": feat,\n+                    "name": CDSname,\n+                }\n+\n+    gff_info = OrderedDict(sorted(gff_info.items(), key=lambda k: k[1]["start"]))\n+    # endBase = 0\n+    for i, feat_id in enumerate(gff_info):\n+        gff_info[feat_id].update({"index": i})\n+        if gff_info[feat_id]["loc"].end > endBase:\n+            endBase = gff_info[feat_id]["loc"].end\n+\n+    re'..b'                   / abs(\n+                        cluster_elem["query_range"][1] - cluster_elem["query_range"][0]\n+                    )\n+                )\n+\n+                tempLoc = FeatureLocation(\n+                    cds.location.start + (3 * (cluster_elem["query_range"][0] - 1)),\n+                    cds.location.start + (3 * (cluster_elem["query_range"][1])),\n+                    cds.location.strand,\n+                )\n+                cds.location = tempLoc\n+                # Set the qualifiers appropriately\n+                cds.qualifiers = {\n+                    "ID": ["gp_%s.CDS.%s" % (cluster_idx, idx)],\n+                    "score": score,\n+                    "Name": self.gff_info[gene_name]["name"],\n+                    "evalue": cluster_elem["evalue"],\n+                    "Identity": cluster_elem["identity_percent"] * 100,\n+                    #\'|\'.join(cluster_elem[\'gi_nos\']) + "| title goes here."\n+                }\n+                # cds.location.start = cds.location.start +\n+                cdss.append(cds)\n+\n+            # And we attach the things properly.\n+            mRNA.sub_features = cdss\n+            mRNA.location = FeatureLocation(mRNA.location.start, mRNA.location.end, cds.location.strand)\n+            gene.sub_features = [mRNA]\n+            gene.location = FeatureLocation(gene.location.start, gene.location.end, cds.location.strand)\n+            \n+            # And append to our record\n+            rec.features.append(gene)\n+        return rec\n+\n+    def output_xml(self, clusters):\n+        threeLevel = {}\n+        # print((clusters.viewkeys()))\n+        # print(type(enumerate(clusters)))\n+        # print(type(clusters))\n+        for cluster_idx, cluster_id in enumerate(clusters):\n+            # print(type(cluster_id))\n+            # print(type(cluster_idx))\n+            # print(type(clusters[cluster_id][0][\'hit_num\']))\n+            if not (clusters[cluster_id][0]["iter_num"] in threeLevel.keys):\n+                threeLevel[clusters[cluster_id][0]["iter_num"]] = {}\n+        # for cluster_idx, cluster_id in enumerate(clusters):\n+        #    print(type(clusters[cluster_id]))\n+        #    b = {clusters[cluster_id][i]: clusters[cluster_id][i+1] for i in range(0, len(clusters[cluster_id]), 2)}\n+        #    print(type(b))#[\'name\']))\n+        # for hspList in clusters:\n+        # for x, idx in (enumerate(clusters)):#for hsp in hspList:\n+        #    print("In X")\n+\n+\n+if __name__ == "__main__":\n+    parser = argparse.ArgumentParser(description="Intron detection")\n+    parser.add_argument("gff3", type=argparse.FileType("r"), help="GFF3 gene calls")\n+    parser.add_argument(\n+        "blastp", type=argparse.FileType("r"), help="blast XML protein results"\n+    )\n+    parser.add_argument(\n+        "--minimum",\n+        help="Gap minimum (Default -1, set to a negative number to allow overlap)",\n+        default=-1,\n+        type=int,\n+    )\n+    parser.add_argument(\n+        "--maximum",\n+        help="Gap maximum in genome (Default 10000)",\n+        default=10000,\n+        type=int,\n+    )\n+    parser.add_argument(\n+        "--idThresh", help="ID Percent Threshold", default=0.4, type=float\n+    )\n+\n+    args = parser.parse_args()\n+\n+    threshCap = args.idThresh\n+    if threshCap > 1.00:\n+        threshCap = 1.00\n+    if threshCap < 0:\n+        threshCap = 0\n+\n+    # create new IntronFinder object based on user input\n+    ifinder = IntronFinder(args.gff3, args.blastp, threshCap)\n+    ifinder.create_clusters()\n+    ifinder.clusters = ifinder.check_strand()\n+    ifinder.clusters = ifinder.check_gene_gap(maximum=args.maximum)\n+    ifinder.clusters = ifinder.check_seq_overlap(minimum=args.minimum)\n+    # ifinder.output_xml(ifinder.clusters)\n+    # for x, idx in (enumerate(ifinder.clusters)):\n+    # print(ifinder.blast)\n+\n+    condensed_report = ifinder.cluster_report()\n+    rec = ifinder.output_gff3(ifinder.clusters)\n+    gffWrite([rec], sys.stdout)\n+\n+    # import pprint; pprint.pprint(ifinder.clusters)\n'
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/intron_detection.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/intron_detection.xml Fri May 13 05:08:54 2022 +0000
[
@@ -0,0 +1,57 @@
+<?xml version="1.1"?>
+<tool id="edu.tamu.cpt2.phage.intron_detection" name="Interrupted gene detection tool" version="19.1.0.0">
+  <description>based on nearby protein blast results</description>
+  <macros>
+    <import>macros.xml</import>
+ <import>cpt-macros.xml</import>
+  </macros>
+  <requirements>
+    <requirement type="package" version="3.6">python</requirement>
+    <requirement type="package" version="1.77">biopython</requirement>
+    <requirement type="package" version="1.1.3">cpt_gffparser</requirement>  
+  </requirements>
+  <command detect_errors="aggressive"><![CDATA[
+python $__tool_directory__/intron_detection.py
+
+"$gff3_data"
+"$blast"
+--minimum $minimum
+--maximum $maximum
+--idThresh $idThresh
+> $output]]></command>
+  <inputs>
+    <expand macro="gff3_input" />
+    <param label="Blastp Results" name="blast" type="data" format="blastxml"/>
+    <param label="Separation Minimum" name="minimum" type="text" value="-1"
+        help="Minimum number of amino acids a CDS pair must be separated by to detect an intron. Set to a negative number to instead allow CDS to overlap by up to that many AAs."/>
+    <param label="Separation Maximum" name="maximum" type="integer" value="10000"
+        help="Maximum number of bases a CDS pair may be separated by to detect an intron (eg, distance from end of 1 cds to the start of another)."/>
+    <param label="Percent Identity Threshold" name="idThresh" type="float" value=".30"
+        help="Minimum Identity Percentage for HSP to be considered for intron comparison (0.00 to 1.00)"/>
+  </inputs>
+  <outputs>
+    <data format="gff3" name="output" label="Potential interrupted genes"/>
+  </outputs>
+  <tests>
+      <test>
+          <param name="gff3_data" value="T7_CLEAN.gff3"/>
+          <param name="blast" value="T7_NR.blastxml"/>
+          <param name="minimum" value="10" />
+          <param name="maximum" value="2000" />
+          <param name="idThresh" value="0.3" />
+          <output name="output" file="T7_INTRON.gff3"/>
+      </test>
+  </tests>
+  <help><![CDATA[
+**What it does**
+
+This tool parses BLAST XML results and returns potential intron regions by 
+identifying query proteins with BLAST hits to the same target protein within 
+the maximum number of bases set above. This distance will also wraparound the 
+ends of the genome. The list of hits is returned as a gff3.
+
+(Tool formerly named Intron Detection, updated Nov 2020)
+
+      ]]></help>
+ <expand macro="citations-2020" />
+</tool>
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/macros.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/macros.xml Fri May 13 05:08:54 2022 +0000
b
@@ -0,0 +1,85 @@
+<?xml version="1.0"?>
+<macros>
+ <xml name="requirements">
+ <requirements>
+ <requirement type="package" version="3.6">python</requirement>
+ <requirement type="package" version="1.77">biopython</requirement>
+ <requirement type="package" version="1.1.3">cpt_gffparser</requirement>  
+ <yield/>
+ </requirements>
+ </xml>
+ <token name="@BLAST_TSV@">
+ "$blast_tsv"
+ </token>
+ <xml name="blast_tsv">
+ <param label="Blast Results" help="TSV/tabular (25 Column)"
+ name="blast_tsv" type="data" format="tabular" />
+ </xml>
+
+ <token name="@BLAST_XML@">
+ "$blast_xml"
+ </token>
+ <xml name="blast_xml">
+ <param label="Blast Results" help="XML format"
+ name="blast_xml" type="data" format="blastxml" />
+ </xml>
+ <xml name="gff3_with_fasta">
+ <param label="Genome Sequences" name="fasta" type="data" format="fasta" />
+ <param label="Genome Annotations" name="gff3" type="data" format="gff3" />
+ </xml>
+ <xml name="genome_selector">
+ <conditional name="reference_genome">
+ <param name="reference_genome_source" type="select" label="Reference Genome">
+ <option value="history" selected="True">From History</option>
+ <option value="cached">Locally Cached</option>
+ </param>
+ <when value="cached">
+ <param name="fasta_indexes" type="select" label="Source FASTA Sequence">
+ <options from_data_table="all_fasta"/>
+ </param>
+ </when>
+ <when value="history">
+ <param name="genome_fasta" type="data" format="fasta" label="Source FASTA Sequence"/>
+ </when>
+ </conditional>
+ </xml>
+ <xml name="gff3_input">
+ <param label="GFF3 Annotations" name="gff3_data" type="data" format="gff3"/>
+ </xml>
+ <xml name="input/gff3+fasta">
+ <expand macro="gff3_input" />
+ <expand macro="genome_selector" />
+ </xml>
+ <token name="@INPUT_GFF@">
+ "$gff3_data"
+ </token>
+ <token name="@INPUT_FASTA@">
+#if str($reference_genome.reference_genome_source) == 'cached':
+ "${reference_genome.fasta_indexes.fields.path}"
+#else if str($reference_genome.reference_genome_source) == 'history':
+ genomeref.fa
+#end if
+ </token>
+ <token name="@GENOME_SELECTOR_PRE@">
+#if $reference_genome.reference_genome_source == 'history':
+ ln -s $reference_genome.genome_fasta genomeref.fa;
+#end if
+ </token>
+ <token name="@GENOME_SELECTOR@">
+#if str($reference_genome.reference_genome_source) == 'cached':
+ "${reference_genome.fasta_indexes.fields.path}"
+#else if str($reference_genome.reference_genome_source) == 'history':
+ genomeref.fa
+#end if
+ </token>
+        <xml name="input/fasta">
+ <param label="Fasta file" name="sequences" type="data" format="fasta"/>
+ </xml>
+
+ <token name="@SEQUENCE@">
+ "$sequences"
+ </token>
+ <xml name="input/fasta/protein">
+ <param label="Protein fasta file" name="sequences" type="data" format="fasta"/>
+ </xml>
+</macros>
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/test-data/T7_CLEAN.gff3
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/test-data/T7_CLEAN.gff3 Fri May 13 05:08:54 2022 +0000
b
b'@@ -0,0 +1,172 @@\n+##gff-version 3\n+##sequence-region NC_001604 1 39937\n+NC_001604\tGenBank\tcontig\t1\t39937\t.\t+\t1\tDbxref=BioProject:PRJNA485481,taxon:10760;ID=NC_001604;Name=NC_001604;Note=Enterobacteria phage T7%2C complete genome.,VALIDATED REFSEQ: This record has undergone validation or preliminary review. The reference sequence was derived from V01146. The sequence was submitted by the authors %5B1%5D on magnetic tape and revised according to %5B3%5D,%5B4%5D,and %5B5%5D. %5B3%5D made changes at 8 positions in gene 1 without affecting the size of the total sequence but changing gene 1 amino acids 443,474,and 388 to 424. %5B4%5D inserted a T at nucleotide 17511,increasing the total sequence to 39937 bp. This change,originally found in T3 DNA %5B8%5D,revealed gene 5.9 and shortened gene 6. %5B5%5D changed the nucleotides at 11061 and 11062 from GT to TG,changing amino acid 119 of T7 lysozyme %28gene 3.5%29 from glycine to valine. Features have been extracted from %5B1%5D unless otherwise noted. The sequence shown is that of the l strand,which corresponds to the sequence of all mRNAs of known functional significance. Early mRNAs are produced by three major promoters for E. coli RNA polymerase A1,A2,and A3,located near the left end of the DNA. A fourth major E. coli promoter,A0 %28also called D%29,that would direct transcription leftward,and several minor E. coli promoters %28see Table 6 in %5B1%5D%29 function in vitro but have no known in vivo function. Late mRNAs are produced by 15 promoters for T7 RNA polymerase distributed across the right-most 85%25 of the DNA,and named e.g. phi10,for the first gene downstream of the promoter. There are also two T7 promoters,phiOL and phiOR,associated with possible origins of replication at the left and right ends of T7 DNA. The 23 base-pair consensus sequence for T7 promoters stretches from -17 to %2B6,where the initiating nucleotide is at %2B1. T7 DNA also contains a 160 base-pair terminal repetition. The beginning and end of RNAs are determined by the promoters,by a terminator for E. coli RNA polymerase,TE,located at the end of the early region,a terminator for T7 RNA polymerase,Tphi,located just downstream of gene 10,and a series of RNase III cleavage sites. Early mRNAs made by E. coli RNA polymerase are listed in Features. The many RNAs predicted to be made by T7 RNA polymerase are not listed but can be deduced from the position of the transcription signals %28see Tables 8 and 9 in %5B1%5D%29. Promoters are listed in Features by the known or predicted first nucleotide of the RNA,terminators by the last nucleotide of the RNA,and RNase III sites by the nucleotide 5%27 of the position of cleavage. Genes are numbered 0.3 to 19.5 in order of their left-to-right position on the genome. Proteins are named by the gene number,e.g.,the gene 1 protein,or by a functional name,T7 RNA polymerase. There is now genetic or biochemical evidence that proteins are produced from at least 52 of the 56 T7 genes. Gene 4 produces two proteins,4A and 4B,by initiating translation at two different sites in the same reading frame. Gene 10 produces two proteins,10A and 10B,by frameshifting during translation. Genes 0.6 and 5.5 probably also make two proteins by translational frameshifting,the gene 5.5 frameshift producing a gene 5.5-5.7 fusion protein. COMPLETENESS: full length.;comment1=VALIDATED REFSEQ: This record has undergone validation or preliminary review. The reference sequence was derived from V01146. The sequence was submitted by the authors %5B1%5D on magnetic tape and revised according to %5B3%5D%2C %5B4%5D%2C and %5B5%5D. %5B3%5D made changes at 8 positions in gene 1 without affecting the size of the total sequence but changing gene 1 amino acids 443%2C 474%2C and 388 to 424. %5B4%5D inserted a T at nucleotide 17511%2C increasing the total sequence to 39937 bp. This change%2C originally found in T3 DNA %5B8%5D%2C revealed gene 5.9 and shortened gene 6. %5B5%5D changed the nucleotides at 11061 and 1106'..b'mal lysis time. However%2C in phage T7 the holin protein gp17.5 does not appear to be essential and gp17.5 mutants only show a minor delay in lysis. Other names: gp17.5%3B lysis protein;codon_start=1;product=type II holin;protein_id=NP_042006.1;transl_table=11;translation=length.67\n+NC_001604\tGenBank\tgene\t36344\t36547\t.\t+\t1\tAlias=T7p53;Dbxref=GeneID:1261022;ID=T7p53.gene;Name=T7p53;Note=gene 17.5\n+NC_001604\tGenBank\tCDS\t36553\t36822\t.\t+\t1\tDbxref=GOA:P03693,UniProtKB/Swiss-Prot:P03693,GeneID:1261042;ID=T7p54;Name=T7p54;Note=involved in the packaging of genome monomers into a procapsid using head-to-tail concatemers of genomes. other names: DNA packaging protein A%3B DNA maturation protein A%3B terminase%2C small subunit;codon_start=1;product=DNA packaging protein%2C small subunit;protein_id=NP_042007.1;transl_table=11;translation=length.89\n+NC_001604\tGenBank\tregulatory\t36836\t36836\t.\t+\t1\tID=GenBank:regulatory:NC_001604:36836:36836;Note=E. coli promoter E%5B6%5D;regulatory_class=promoter\n+NC_001604\tGenBank\tsequence_secondary_structure\t36856\t36856\t.\t+\t1\tID=GenBank:sequence_secondary_structure:NC_001604:36856:36856;Note=RNase III site R18.5\n+NC_001604\tGenBank\tgene\t36553\t36822\t.\t+\t1\tAlias=T7p54;Dbxref=GeneID:1261042;ID=T7p54.gene;Name=T7p54;Note=gene 18\n+NC_001604\tGenBank\tCDS\t36917\t37348\t.\t+\t1\tDbxref=GOA:P03803,UniProtKB/Swiss-Prot:P03803,GeneID:1261067;ID=T7p55;Name=T7p55;Note=analog of phage lambda protein Rz%2C a cell lysis protein. Rz and gp18.5 share distant sequence similarity%2C similar function%2C and a similar genome neighborhood. In T7%2C gp18.5 interacts with gp18.7%2C a lambda RZ1-like lysis protein. Other names: gp18.5;codon_start=1;product=phage lambda Rz-like lysis protein;protein_id=NP_042008.1;transl_table=11;translation=length.143\n+NC_001604\tGenBank\tgene\t36917\t37348\t.\t+\t1\tAlias=T7p55;Dbxref=GeneID:1261067;ID=T7p55.gene;Name=T7p55;Note=gene 18.5\n+NC_001604\tGenBank\tCDS\t37032\t37283\t.\t+\t1\tDbxref=UniProtKB/Swiss-Prot:P03788,GeneID:1261057;ID=T7p56;Name=T7p56;Note=in Enterobacteria phage T7%2C this protein interacts with gp18.5 and is expressed from the -1 frame of a gene completely overlapping gene 18.5. This suggests that it may be an analog of lambda lysis protein Rz1. Other names: gp18.7.;codon_start=1;product=phage lambda Rz1-like protein;protein_id=NP_042009.1;transl_table=11;translation=length.83\n+NC_001604\tGenBank\tgene\t37032\t37283\t.\t+\t1\tAlias=T7p56;Dbxref=GeneID:1261057;ID=T7p56.gene;Name=T7p56;Note=gene 18.7\n+NC_001604\tGenBank\tCDS\t37370\t39130\t.\t+\t1\tDbxref=GOA:P03694,UniProtKB/Swiss-Prot:P03694,GeneID:1261062;ID=T7p57;Name=T7p57;Note=gene 19;codon_start=1;product=DNA maturation protein;protein_id=NP_042010.1;transl_table=11;translation=length.586\n+NC_001604\tGenBank\tgene\t37370\t39130\t.\t+\t1\tAlias=T7p57;Dbxref=GeneID:1261062;ID=T7p57.gene;Name=T7p57;Note=gene 19\n+NC_001604\tGenBank\tCDS\t38016\t38273\t.\t+\t1\tDbxref=UniProtKB/Swiss-Prot:P03789,GeneID:1261064;ID=T7p58;Name=T7p58;Note=gene 19.2;codon_start=1;product=hypothetical protein;protein_id=NP_042011.1;transl_table=11;translation=length.85\n+NC_001604\tGenBank\tgene\t38016\t38273\t.\t+\t1\tAlias=T7p58;Dbxref=GeneID:1261064;ID=T7p58.gene;Name=T7p58;Note=gene 19.2\n+NC_001604\tGenBank\tCDS\t38553\t38726\t.\t+\t1\tDbxref=UniProtKB/Swiss-Prot:P03790,GeneID:1261066;ID=T7p59;Name=T7p59;Note=gene 19.3;codon_start=1;product=hypothetical protein;protein_id=NP_042012.1;transl_table=11;translation=length.57\n+NC_001604\tGenBank\tregulatory\t39229\t39229\t.\t+\t1\tID=GenBank:regulatory:NC_001604:39229:39229;Note=T7 promoter phiOR;regulatory_class=promoter\n+NC_001604\tGenBank\tgene\t38553\t38726\t.\t+\t1\tAlias=T7p59;Dbxref=GeneID:1261066;ID=T7p59.gene;Name=T7p59;Note=gene 19.3\n+NC_001604\tGenBank\tCDS\t39389\t39538\t.\t+\t1\tDbxref=UniProtKB/Swiss-Prot:P03804,GeneID:1261068;ID=T7p60;Name=T7p60;Note=gene 19.5;codon_start=1;product=hypothetical protein;protein_id=NP_042013.1;transl_table=11;translation=length.49\n+NC_001604\tGenBank\tgene\t39389\t39538\t.\t+\t1\tAlias=T7p60;Dbxref=GeneID:1261068;ID=T7p60.gene;Name=T7p60;Note=gene 19.5\n'
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/test-data/T7_INTRON.gff3
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/test-data/T7_INTRON.gff3 Fri May 13 05:08:54 2022 +0000
[
@@ -0,0 +1,5 @@
+##gff-version 3
+NC_001604 feature gene 9158 10706 . + . ID=gp_0;Percent_Identities=T7p17 had 55%25 identity to NCBI Protein ID QFG06980,T7p19 had 70%25 identity to NCBI Protein ID QFG06980;Note=gb|QFG06980.1| putative helix-destabilizing protein [Escherichia virus ECG4];
+NC_001604 feature mRNA 9158 10706 . + . ID=gp_0.mRNA;note=T7p17 had 55%25 identity to NCBI Protein ID QFG06980,T7p19 had 70%25 identity to NCBI Protein ID QFG06980;Parent=gp_0;
+NC_001604 GenBank CDS 9317 9697 . + 1 ID=gp_0.CDS.0;score=563;Name=T7p17;evalue=6.62467e-40;Identity=55.55555555555556;Parent=gp_0.mRNA;
+NC_001604 GenBank CDS 10290 10649 . + 1 ID=gp_0.CDS.1;score=714;Name=T7p19;evalue=5.74854e-55;Identity=70.58823529411765;Parent=gp_0.mRNA;
b
diff -r 000000000000 -r 1a19092729be cpt_intron_detect/test-data/T7_NR.blastxml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_intron_detect/test-data/T7_NR.blastxml Fri May 13 05:08:54 2022 +0000
[
b'@@ -0,0 +1,379010 @@\n+<?xml version="1.0"?>\n+<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">\n+<BlastOutput>\n+  <BlastOutput_program>blastp</BlastOutput_program>\n+  <BlastOutput_version>BLASTP 2.9.0+</BlastOutput_version>\n+  <BlastOutput_reference>Stephen F. Altschul, Thomas L. Madden, Alejandro A. Sch&amp;auml;ffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), &quot;Gapped BLAST and PSI-BLAST: a new generation of protein database search programs&quot;, Nucleic Acids Res. 25:3389-3402.</BlastOutput_reference>\n+  <BlastOutput_db>/opt/blast/nr_v5/nr</BlastOutput_db>\n+  <BlastOutput_query-ID>Query_1</BlastOutput_query-ID>\n+  <BlastOutput_query-def>T7p01.p01 [Location=[924:1278](+);Name=T7p01]</BlastOutput_query-def>\n+  <BlastOutput_query-len>117</BlastOutput_query-len>\n+  <BlastOutput_param>\n+    <Parameters>\n+      <Parameters_matrix>BLOSUM62</Parameters_matrix>\n+      <Parameters_expect>0.001</Parameters_expect>\n+      <Parameters_gap-open>11</Parameters_gap-open>\n+      <Parameters_gap-extend>1</Parameters_gap-extend>\n+      <Parameters_filter>F</Parameters_filter>\n+    </Parameters>\n+  </BlastOutput_param>\n+<BlastOutput_iterations>\n+<Iteration>\n+  <Iteration_iter-num>1</Iteration_iter-num>\n+  <Iteration_query-ID>Query_1</Iteration_query-ID>\n+  <Iteration_query-def>T7p01.p01 [Location=[924:1278](+);Name=T7p01]</Iteration_query-def>\n+  <Iteration_query-len>117</Iteration_query-len>\n+<Iteration_hits>\n+<Hit>\n+  <Hit_num>1</Hit_num>\n+  <Hit_id>ref|NP_041954.1|</Hit_id>\n+  <Hit_def>hypothetical protein T7p01 [Escherichia phage T7] &gt;sp|P03775.1| RecName: Full=Overcome classical restriction gp0.3; AltName: Full=Gene product 0.3; Short=Gp0.3; AltName: Full=OCR [Escherichia phage T7] &gt;gb|AAZ32821.1| protein 0.3 [Enterobacteria phage T7.1] &gt;gb|AXQ60539.1| hypothetical protein [Synthetic phage] &gt;gb|AAP33909.1| gene 0.3 [Escherichia phage T7] &gt;gb|AVH85518.1| hypothetical protein mT7_p01 [Escherichia phage T7] &gt;emb|CAA24327.1| unnamed protein product [Escherichia phage T7]</Hit_def>\n+  <Hit_accession>NP_041954</Hit_accession>\n+  <Hit_len>117</Hit_len>\n+  <Hit_hsps>\n+    <Hsp>\n+      <Hsp_num>1</Hsp_num>\n+      <Hsp_bit-score>233.802</Hsp_bit-score>\n+      <Hsp_score>595</Hsp_score>\n+      <Hsp_evalue>1.54426e-77</Hsp_evalue>\n+      <Hsp_query-from>1</Hsp_query-from>\n+      <Hsp_query-to>117</Hsp_query-to>\n+      <Hsp_hit-from>1</Hsp_hit-from>\n+      <Hsp_hit-to>117</Hsp_hit-to>\n+      <Hsp_query-frame>0</Hsp_query-frame>\n+      <Hsp_hit-frame>0</Hsp_hit-frame>\n+      <Hsp_identity>117</Hsp_identity>\n+      <Hsp_positive>117</Hsp_positive>\n+      <Hsp_gaps>0</Hsp_gaps>\n+      <Hsp_align-len>117</Hsp_align-len>\n+      <Hsp_qseq>MAMSNMTYNNVFDHAYEMLKENIRYDDIRDTDDLHDAIHMAADNAVPHYYADIFSVMASEGIDLEFEDSGLMPDTKDVIRILQARIYEQLTIDLWEDAEDLLNEYLEEVEEYEEDEE</Hsp_qseq>\n+      <Hsp_hseq>MAMSNMTYNNVFDHAYEMLKENIRYDDIRDTDDLHDAIHMAADNAVPHYYADIFSVMASEGIDLEFEDSGLMPDTKDVIRILQARIYEQLTIDLWEDAEDLLNEYLEEVEEYEEDEE</Hsp_hseq>\n+      <Hsp_midline>MAMSNMTYNNVFDHAYEMLKENIRYDDIRDTDDLHDAIHMAADNAVPHYYADIFSVMASEGIDLEFEDSGLMPDTKDVIRILQARIYEQLTIDLWEDAEDLLNEYLEEVEEYEEDEE</Hsp_midline>\n+    </Hsp>\n+  </Hit_hsps>\n+</Hit>\n+<Hit>\n+  <Hit_num>2</Hit_num>\n+  <Hit_id>gb|ACY75829.1|</Hit_id>\n+  <Hit_def>0.3 protein [Escherichia phage T7]</Hit_def>\n+  <Hit_accession>ACY75829</Hit_accession>\n+  <Hit_len>117</Hit_len>\n+  <Hit_hsps>\n+    <Hsp>\n+      <Hsp_num>1</Hsp_num>\n+      <Hsp_bit-score>233.032</Hsp_bit-score>\n+      <Hsp_score>593</Hsp_score>\n+      <Hsp_evalue>3.75744e-77</Hsp_evalue>\n+      <Hsp_query-from>1</Hsp_query-from>\n+      <Hsp_query-to>117</Hsp_query-to>\n+      <Hsp_hit-from>1</Hsp_hit-from>\n+      <Hsp_hit-to>117</Hsp_hit-to>\n+      <Hsp_query-frame>0</Hsp_query-frame>\n+      <Hsp_hit-frame>0</Hsp_hit-frame>\n+      <Hsp_identity>116</Hsp_identity>\n+      <Hsp_positive>117</Hsp_positive>\n+      <Hsp_gaps>0</Hsp_gaps>\n+      <Hsp_align-len>117</Hsp_align-len>\n+     '..b' <Hit_num>48</Hit_num>\n+  <Hit_id>gb|AWY10495.1|</Hit_id>\n+  <Hit_def>hypothetical protein PFP1_45 [Pseudomonas phage PFP1]</Hit_def>\n+  <Hit_accession>AWY10495</Hit_accession>\n+  <Hit_len>57</Hit_len>\n+  <Hit_hsps>\n+    <Hsp>\n+      <Hsp_num>1</Hsp_num>\n+      <Hsp_bit-score>43.5134</Hsp_bit-score>\n+      <Hsp_score>101</Hsp_score>\n+      <Hsp_evalue>0.000414021</Hsp_evalue>\n+      <Hsp_query-from>4</Hsp_query-from>\n+      <Hsp_query-to>47</Hsp_query-to>\n+      <Hsp_hit-from>12</Hsp_hit-from>\n+      <Hsp_hit-to>55</Hsp_hit-to>\n+      <Hsp_query-frame>0</Hsp_query-frame>\n+      <Hsp_hit-frame>0</Hsp_hit-frame>\n+      <Hsp_identity>20</Hsp_identity>\n+      <Hsp_positive>30</Hsp_positive>\n+      <Hsp_gaps>0</Hsp_gaps>\n+      <Hsp_align-len>44</Hsp_align-len>\n+      <Hsp_qseq>LLLNLLRHRVTYRFLVVLCAALGYASLTGDLSSLESVVCSILTC</Hsp_qseq>\n+      <Hsp_hseq>VLVSLLKHRATYRFLAVLLVALGVANGEAIMAGVETFVCAFAGC</Hsp_hseq>\n+      <Hsp_midline>+L++LL+HR TYRFL VL  ALG A+    ++ +E+ VC+   C</Hsp_midline>\n+    </Hsp>\n+  </Hit_hsps>\n+</Hit>\n+<Hit>\n+  <Hit_num>49</Hit_num>\n+  <Hit_id>ref|YP_004306360.1|</Hit_id>\n+  <Hit_def>hypothetical protein phiIBB-PF7Ap40 [Pseudomonas phage phiIBB-PF7A] &gt;gb|ADV35705.1| conserved hypothetical protein [Pseudomonas phage phiIBB-PF7A]</Hit_def>\n+  <Hit_accession>YP_004306360</Hit_accession>\n+  <Hit_len>57</Hit_len>\n+  <Hit_hsps>\n+    <Hsp>\n+      <Hsp_num>1</Hsp_num>\n+      <Hsp_bit-score>43.1282</Hsp_bit-score>\n+      <Hsp_score>100</Hsp_score>\n+      <Hsp_evalue>0.000581524</Hsp_evalue>\n+      <Hsp_query-from>4</Hsp_query-from>\n+      <Hsp_query-to>47</Hsp_query-to>\n+      <Hsp_hit-from>12</Hsp_hit-from>\n+      <Hsp_hit-to>55</Hsp_hit-to>\n+      <Hsp_query-frame>0</Hsp_query-frame>\n+      <Hsp_hit-frame>0</Hsp_hit-frame>\n+      <Hsp_identity>20</Hsp_identity>\n+      <Hsp_positive>30</Hsp_positive>\n+      <Hsp_gaps>0</Hsp_gaps>\n+      <Hsp_align-len>44</Hsp_align-len>\n+      <Hsp_qseq>LLLNLLRHRVTYRFLVVLCAALGYASLTGDLSSLESVVCSILTC</Hsp_qseq>\n+      <Hsp_hseq>VLVSLAKHRATYRFLAVLLVALGISNGEAIMSGIETVACAYLGC</Hsp_hseq>\n+      <Hsp_midline>+L++L +HR TYRFL VL  ALG ++    +S +E+V C+ L C</Hsp_midline>\n+    </Hsp>\n+  </Hit_hsps>\n+</Hit>\n+<Hit>\n+  <Hit_num>50</Hit_num>\n+  <Hit_id>gb|ANY29061.1|</Hit_id>\n+  <Hit_def>hypothetical protein UNOSLW4_0230 [Pseudomonas phage UNO-SLW4] &gt;gb|ANY29108.1| hypothetical protein UNOSLW3_0235 [Pseudomonas phage UNO-SLW3] &gt;gb|ANY29155.1| hypothetical protein UNOSLW2_0235 [Pseudomonas phage UNO-SLW2] &gt;gb|ANY29202.1| hypothetical protein UNOSLW1_0235 [Pseudomonas phage UNO-SLW1]</Hit_def>\n+  <Hit_accession>ANY29061</Hit_accession>\n+  <Hit_len>57</Hit_len>\n+  <Hit_hsps>\n+    <Hsp>\n+      <Hsp_num>1</Hsp_num>\n+      <Hsp_bit-score>42.743</Hsp_bit-score>\n+      <Hsp_score>99</Hsp_score>\n+      <Hsp_evalue>0.000853372</Hsp_evalue>\n+      <Hsp_query-from>4</Hsp_query-from>\n+      <Hsp_query-to>47</Hsp_query-to>\n+      <Hsp_hit-from>12</Hsp_hit-from>\n+      <Hsp_hit-to>55</Hsp_hit-to>\n+      <Hsp_query-frame>0</Hsp_query-frame>\n+      <Hsp_hit-frame>0</Hsp_hit-frame>\n+      <Hsp_identity>18</Hsp_identity>\n+      <Hsp_positive>30</Hsp_positive>\n+      <Hsp_gaps>0</Hsp_gaps>\n+      <Hsp_align-len>44</Hsp_align-len>\n+      <Hsp_qseq>LLLNLLRHRVTYRFLVVLCAALGYASLTGDLSSLESVVCSILTC</Hsp_qseq>\n+      <Hsp_hseq>VMVSLLKHRATYRFLAVLLVALGVTNGEAIMAGIETIVCAFTGC</Hsp_hseq>\n+      <Hsp_midline>++++LL+HR TYRFL VL  ALG  +    ++ +E++VC+   C</Hsp_midline>\n+    </Hsp>\n+  </Hit_hsps>\n+</Hit>\n+</Iteration_hits>\n+  <Iteration_stat>\n+    <Statistics>\n+      <Statistics_db-num>251406437</Statistics_db-num>\n+      <Statistics_db-len>90675456348</Statistics_db-len>\n+      <Statistics_hsp-len>22</Statistics_hsp-len>\n+      <Statistics_eff-space>2298901897818</Statistics_eff-space>\n+      <Statistics_kappa>0.041</Statistics_kappa>\n+      <Statistics_lambda>0.267</Statistics_lambda>\n+      <Statistics_entropy>0.14</Statistics_entropy>\n+    </Statistics>\n+  </Iteration_stat>\n+</Iteration>\n+</BlastOutput_iterations>\n+</BlastOutput>\n+\n'