Merge branch 'master' of https://github.com/arvados/bh20-seq-resource

3 files changed, 60 insertions, 63 deletions

diff --git a/bh20sequploader/main.py b/bh20sequploader/main.py
index c442af0..fd0278d 100644
--- a/bh20sequploader/main.py
+++ b/bh20sequploader/main.py
@@ -22,30 +22,32 @@ ARVADOS_API_HOST='lugli.arvadosapi.com'
 ARVADOS_API_TOKEN='2fbebpmbo3rw3x05ueu2i6nx70zhrsb1p22ycu3ry34m4x4462'
 UPLOAD_PROJECT='lugli-j7d0g-n5clictpuvwk8aa'
 
-def qa_stuff(metadata, sequence_p1, sequence_p2):
+def qc_stuff(metadata, sequence_p1, sequence_p2, do_qc=True):
+    failed = False
     try:
-        log.debug("Checking metadata")
-        if not qc_metadata(metadata.name):
+        log.debug("Checking metadata" if do_qc else "Skipping metadata check")
+        if do_qc and not qc_metadata(metadata.name):
             log.warning("Failed metadata qc")
-            exit(1)
-    except ValueError as e:
+            failed = True
+    except Exception as e:
         log.debug(e)
-        log.debug("Failed metadata qc")
         print(e)
-        exit(1)
+        failed = True
 
     target = []
     try:
-        log.debug("Checking FASTA/FASTQ QC")
-        target.append(qc_fasta(sequence_p1))
+        log.debug("FASTA/FASTQ QC" if do_qc else "Limited FASTA/FASTQ QC")
+        target.append(qc_fasta(sequence_p1, check_with_clustalw=do_qc))
         if sequence_p2:
             target.append(qc_fasta(sequence_p2))
             target[0] = ("reads_1."+target[0][0][6:], target[0][1])
             target[1] = ("reads_2."+target[1][0][6:], target[0][1])
-    except ValueError as e:
+    except Exception as e:
         log.debug(e)
-        log.debug("Failed FASTA qc")
         print(e)
+        failed = True
+
+    if failed:
         exit(1)
 
     return target
@@ -62,13 +64,14 @@ def main():
     parser = argparse.ArgumentParser(description='Upload SARS-CoV-19 sequences for analysis')
     parser.add_argument('metadata', type=argparse.FileType('r'), help='sequence metadata json')
     parser.add_argument('sequence_p1', type=argparse.FileType('rb'), help='sequence FASTA/FASTQ')
-    parser.add_argument('sequence_p2', type=argparse.FileType('rb'), default=None, help='sequence FASTQ pair')
+    parser.add_argument('sequence_p2', type=argparse.FileType('rb'), default=None, nargs='?', help='sequence FASTQ pair')
     parser.add_argument("--validate", action="store_true", help="Dry run, validate only")
+    parser.add_argument("--skip-qc", action="store_true", help="Skip local qc check")
     args = parser.parse_args()
 
     api = arvados.api(host=ARVADOS_API_HOST, token=ARVADOS_API_TOKEN, insecure=True)
 
-    target = qa_stuff(args.metadata, args.sequence_p1, args.sequence_p2)
+    target = qc_stuff(args.metadata, args.sequence_p1, args.sequence_p2, not args.skip_qc)
     seqlabel = target[0][1]
 
     if args.validate:
diff --git a/bh20sequploader/qc_fasta.py b/bh20sequploader/qc_fasta.py
index e198430..37eb4e8 100644
--- a/bh20sequploader/qc_fasta.py
+++ b/bh20sequploader/qc_fasta.py
@@ -25,7 +25,7 @@ def read_fasta(sequence):
             raise ValueError("FASTA file contains multiple entries")
     return label, bases
 
-def qc_fasta(arg_sequence):
+def qc_fasta(arg_sequence, check_with_clustalw=True):
     log.debug("Starting qc_fasta")
     schema_resource = pkg_resources.resource_stream(__name__, "validation/formats")
     with tempfile.NamedTemporaryFile() as tmp:
@@ -51,47 +51,48 @@ def qc_fasta(arg_sequence):
     if seq_type == "text/fasta":
         # ensure that contains only one entry
         submitlabel, submitseq = read_fasta(sequence)
+        sequence.seek(0)
+        sequence.detach()
 
-        with tempfile.NamedTemporaryFile() as tmp1:
-            refstring = pkg_resources.resource_string(__name__, "SARS-CoV-2-reference.fasta")
-            tmp1.write(refstring)
-            tmp1.write(submitlabel.encode("utf8"))
-            tmp1.write(("".join(submitseq)).encode("utf8"))
-            tmp1.flush()
-            subbp = 0
-            refbp = 0
-            similarity = 0
-            try:
-                cmd = ["clustalw", "-infile="+tmp1.name,
-                       "-quicktree", "-iteration=none", "-type=DNA"]
-                print("QC checking similarity to reference")
-                print(" ".join(cmd))
-                result = subprocess.run(cmd, stdout=subprocess.PIPE)
-                res = result.stdout.decode("utf-8")
-                g1 = re.search(r"^Sequence 1: [^ ]+ +(\d+) bp$", res, flags=re.MULTILINE)
-                refbp = float(g1.group(1))
-                g2 = re.search(r"^Sequence 2: [^ ]+ +(\d+) bp$", res, flags=re.MULTILINE)
-                subbp = float(g2.group(1))
-                g3 = re.search(r"^Sequences \(1:2\) Aligned\. Score: (\d+(\.\d+)?)$", res, flags=re.MULTILINE)
-                similarity = float(g3.group(1))
+        if not check_with_clustalw:
+            return ("sequence.fasta"+gz, seqlabel)
 
-                print(g1.group(0))
-                print(g2.group(0))
-                print(g3.group(0))
-            except Exception as e:
-                logging.warn("Error trying to QC against reference sequence using 'clustalw': %s", e)
+        with tempfile.NamedTemporaryFile() as tmp1:
+            with tempfile.NamedTemporaryFile() as tmp2:
+                refstring = pkg_resources.resource_string(__name__, "SARS-CoV-2-reference.fasta")
+                tmp1.write(refstring)
+                tmp1.flush()
+                tmp2.write(submitlabel.encode("utf8"))
+                tmp2.write(("".join(submitseq)).encode("utf8"))
+                tmp2.flush()
+                subbp = 0
+                refbp = 0
+                similarity = 0
+                try:
+                    cmd = ["minimap2", "-c", tmp1.name, tmp2.name]
+                    logging.info("QC checking similarity to reference")
+                    logging.info(" ".join(cmd))
+                    result = subprocess.run(cmd, stdout=subprocess.PIPE)
+                    result.check_returncode()
+                    res = result.stdout.decode("utf-8")
+                    mm = res.split("\t")
+                    if len(mm) >= 10:
+                        # divide Number of matching bases in the mapping / Target sequence length
+                        similarity = (float(mm[9]) / float(mm[6])) * 100.0
+                    else:
+                        similarity = 0
+                except Exception as e:
+                    logging.warn("QC against reference sequence using 'minimap2': %s", e, exc_info=e)
 
-            if refbp and (subbp/refbp) < .7:
-                raise ValueError("QC fail: submit sequence length is shorter than 70% reference")
-            if refbp and (subbp/refbp) > 1.3:
-                raise ValueError("QC fail: submit sequence length is greater than 130% reference")
-            if similarity and similarity < 70.0:
-                raise ValueError("QC fail: submit similarity is less than 70%")
-            if refbp == 0 or similarity == 0:
-                raise ValueError("QC fail")
+                if similarity and similarity < 70.0:
+                    raise ValueError("QC fail: alignment to reference was less than 70%% (was %2.2f%%)" % (similarity))
+                if similarity == 0:
+                    raise ValueError("QC fail")
 
         return ("sequence.fasta"+gz, seqlabel)
     elif seq_type == "text/fastq":
+        sequence.seek(0)
+        sequence.detach()
         return ("reads.fastq"+gz, seqlabel)
     else:
         raise ValueError("Sequence file does not look like a DNA FASTA or FASTQ")
diff --git a/bh20sequploader/qc_metadata.py b/bh20sequploader/qc_metadata.py
index 9122ace..2b57991 100644
--- a/bh20sequploader/qc_metadata.py
+++ b/bh20sequploader/qc_metadata.py
@@ -21,20 +21,13 @@ def qc_metadata(metadatafile):
     shex = pkg_resources.resource_stream(__name__, "bh20seq-shex.rdf").read().decode("utf-8")
 
     if not isinstance(avsc_names, schema_salad.avro.schema.Names):
-        print(avsc_names)
-        return False
+        raise Exception(avsc_names)
 
-    try:
-        doc, metadata = schema_salad.schema.load_and_validate(document_loader, avsc_names, metadatafile, True)
-        g = schema_salad.jsonld_context.makerdf("workflow", doc, document_loader.ctx)
-        rslt, reason = evaluate(g, shex, doc["id"], "https://raw.githubusercontent.com/arvados/bh20-seq-resource/master/bh20sequploader/bh20seq-shex.rdf#submissionShape")
+    doc, metadata = schema_salad.schema.load_and_validate(document_loader, avsc_names, metadatafile, True)
+    g = schema_salad.jsonld_context.makerdf("workflow", doc, document_loader.ctx)
+    rslt, reason = evaluate(g, shex, doc["id"], "https://raw.githubusercontent.com/arvados/bh20-seq-resource/master/bh20sequploader/bh20seq-shex.rdf#submissionShape")
 
-        if not rslt:
-            log.debug(reason)
-            print(reason)
+    if not rslt:
+        raise Exception(reason)
 
-        return rslt
-    except Exception as e:
-        traceback.print_exc()
-        log.warn(e)
-    return False
+    return True