Better handling of duplicate sequences

Also save original fasta label in metadata

7 files changed, 59 insertions, 16 deletions

diff --git a/workflows/pangenome-generate/merge-metadata.cwl b/workflows/pangenome-generate/merge-metadata.cwl
index 9164c09..fcefe32 100644
--- a/workflows/pangenome-generate/merge-metadata.cwl
+++ b/workflows/pangenome-generate/merge-metadata.cwl
@@ -7,6 +7,8 @@ inputs:
   metadata: File[]
   metadataSchema: File
   subjects: string[]
+  dups: File?
+  originalLabels: File
 outputs:
   merged: stdout
 stdout: mergedmetadata.ttl
diff --git a/workflows/pangenome-generate/merge-metadata.py b/workflows/pangenome-generate/merge-metadata.py
index 64275b1..bfec781 100644
--- a/workflows/pangenome-generate/merge-metadata.py
+++ b/workflows/pangenome-generate/merge-metadata.py
@@ -1,9 +1,13 @@
+import re
 import schema_salad.schema
 import schema_salad.jsonld_context
+import json
 
 metadataSchema = '$(inputs.metadataSchema.path)'
 metadata = $(inputs.metadata)
 subjects = $(inputs.subjects)
+dups = json.loads('''$(inputs.dups)''')
+originalLabels = $(inputs.originalLabels)
 
 (document_loader,
  avsc_names,
@@ -11,7 +15,22 @@ subjects = $(inputs.subjects)
  metaschema_loader) = schema_salad.schema.load_schema(metadataSchema)
 
 for i, m in enumerate(metadata):
-    doc, metadata = schema_salad.schema.load_and_validate(document_loader, avsc_names, m["path"], True)
+    doc, metadata = schema_salad.schema.load_and_validate(document_loader, avsc_names, m["path"], False, False)
     doc["id"] = subjects[i]
     g = schema_salad.jsonld_context.makerdf(subjects[i], doc, document_loader.ctx)
     print(g.serialize(format="ntriples").decode("utf-8"))
+
+import logging
+
+if dups:
+    sameseqs = open(dups["path"], "rt")
+    for d in sameseqs:
+        logging.warn(d)
+        g = re.match(r"\\d+\\t(.*)", d)
+        logging.warn("%s", g.group(1))
+        sp = g.group(1).split(",")
+        for n in sp[1:]:
+            print("<%s> <http://biohackathon.org/bh20-seq-schema/has_duplicate_sequence> <%s> ." % (n.strip(), sp[0].strip()))
+
+orig = open(originalLabels["path"], "rt")
+print(orig.read())
diff --git a/workflows/pangenome-generate/pangenome-generate.cwl b/workflows/pangenome-generate/pangenome-generate.cwl
index 896f936..0cb1368 100644
--- a/workflows/pangenome-generate/pangenome-generate.cwl
+++ b/workflows/pangenome-generate/pangenome-generate.cwl
@@ -26,15 +26,11 @@ steps:
     in:
       readsFA: inputReads
       subjects: subjects
-    out: [relabeledSeqs]
+    out: [relabeledSeqs, originalLabels]
     run: relabel-seqs.cwl
-  common:
-    in: {readsFA: relabel/relabeledSeqs}
-    out: [duplicatedReads]
-    run: seqkit-common.cwl
   dedup:
     in: {readsFA: relabel/relabeledSeqs}
-    out: [readsMergeDedup]
+    out: [readsMergeDedup, dups]
     run: seqkit-rmdup.cwl
   overlapReads:
     in: {readsFA: dedup/readsMergeDedup}
@@ -63,5 +59,7 @@ steps:
       metadata: metadata
       metadataSchema: metadataSchema
       subjects: subjects
+      dups: dedup/dups
+      originalLabels: relabel/originalLabels
     out: [merged]
     run: merge-metadata.cwl
diff --git a/workflows/pangenome-generate/relabel-seqs.cwl b/workflows/pangenome-generate/relabel-seqs.cwl
index b5b7231..2b780d4 100644
--- a/workflows/pangenome-generate/relabel-seqs.cwl
+++ b/workflows/pangenome-generate/relabel-seqs.cwl
@@ -5,7 +5,13 @@ inputs:
   subjects: string[]
 outputs:
   relabeledSeqs:
-    type: stdout
+    type: File
+    outputBinding:
+      glob: relabeledSeqs.fasta
+  originalLabels:
+    type: File
+    outputBinding:
+      glob: originalLabels.ttl
 requirements:
   InlineJavascriptRequirement: {}
   InitialWorkDirRequirement:
@@ -15,5 +21,5 @@ requirements:
 hints:
   DockerRequirement:
     dockerPull: commonworkflowlanguage/cwltool_module
-stdout: relabeledSeqs.fasta
+stdout:
 baseCommand: [python, relabel-seqs.py]
diff --git a/workflows/pangenome-generate/relabel-seqs.py b/workflows/pangenome-generate/relabel-seqs.py
index 32f2386..b558fe2 100644
--- a/workflows/pangenome-generate/relabel-seqs.py
+++ b/workflows/pangenome-generate/relabel-seqs.py
@@ -1,13 +1,15 @@
-import sys
-
 reads = $(inputs.readsFA)
 subjects = $(inputs.subjects)
 
+relabeled_fasta = open("relabeledSeqs.fasta", "wt")
+original_labels = open("originalLabels.ttl", "wt")
+
 for i, r in enumerate(reads):
     with open(r["path"], "rt") as fa:
-        fa.readline()
-        print(">"+subjects[i])
+        label = fa.readline()
+        original_labels.write("<%s> <http://biohackathon.org/bh20-seq-schema/original_fasta_label> \\"%s\\" .\\n" % (subjects[i], label[1:].strip().replace('"', '\\\\"')))
+        relabeled_fasta.write(">"+subjects[i]+"\\n")
         data = fa.read(8096)
         while data:
-            sys.stdout.write(data)
+            relabeled_fasta.write(data)
             data = fa.read(8096)
diff --git a/workflows/pangenome-generate/seqkit-rmdup.cwl b/workflows/pangenome-generate/seqkit-rmdup.cwl
index 07184c3..071fa66 100644
--- a/workflows/pangenome-generate/seqkit-rmdup.cwl
+++ b/workflows/pangenome-generate/seqkit-rmdup.cwl
@@ -1,14 +1,14 @@
 cwlVersion: v1.1
 class: CommandLineTool
 inputs:
-  readsFA: File[]
+  readsFA: File
 outputs:
   readsMergeDedup:
     type: File
     outputBinding:
       glob: readsMergeDedup.fasta
   dups:
-    type: File
+    type: File?
     outputBinding:
       glob: dups.txt
 requirements:
diff --git a/workflows/pangenome-generate/testjob.yml b/workflows/pangenome-generate/testjob.yml
new file mode 100644
index 0000000..a48aff8
--- /dev/null
+++ b/workflows/pangenome-generate/testjob.yml
@@ -0,0 +1,16 @@
+inputReads:
+  - class: File
+    location: ../../example/sequence.fasta
+  - class: File
+    location: ../../example/sequence.fasta
+metadata:
+  - class: File
+    location: ../../example/metadata.yaml
+  - class: File
+    location: ../../example/metadata.yaml
+metadataSchema:
+  class: File
+  location: ../../bh20sequploader/bh20seq-schema.yml
+subjects:
+  - http://arvados.org/keep/seq1
+  - http://arvados.org/keep/seq2