From 9ddcfeacb3191638f42b08af999889d867f0f81c Mon Sep 17 00:00:00 2001 From: Peter Amstutz Date: Mon, 20 Apr 2020 14:57:25 -0400 Subject: Better handling of duplicate sequences Also save original fasta label in metadata --- workflows/pangenome-generate/merge-metadata.cwl | 2 ++ workflows/pangenome-generate/merge-metadata.py | 21 ++++++++++++++++++++- workflows/pangenome-generate/pangenome-generate.cwl | 10 ++++------ workflows/pangenome-generate/relabel-seqs.cwl | 10 ++++++++-- workflows/pangenome-generate/relabel-seqs.py | 12 +++++++----- workflows/pangenome-generate/seqkit-rmdup.cwl | 4 ++-- workflows/pangenome-generate/testjob.yml | 16 ++++++++++++++++ 7 files changed, 59 insertions(+), 16 deletions(-) create mode 100644 workflows/pangenome-generate/testjob.yml (limited to 'workflows/pangenome-generate') 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> <%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> \\"%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 -- cgit v1.2.3