label source id duplicates, don't remove, keep UK-ENG format for public metadata

Author: Rachel Colquhoun

bin/add_to_uk_metadata.py

@@ -48,6 +48,11 @@ def add_sample_date(row, date_dict):
         except:
             row["sample_date"] = ""
 
+def add_source_id(row):
+    row["source_id"] = row["biosample_source_id"]
+    if row["root_biosample_source_id"] not in [None,""]:
+        row["source_id"] = row["root_biosample_source_id"]
+
 def add_pillar_2(row):
     if row['collection_pillar'] == 2 or row['central_sample_id'][0:4] in ["ALDP", "CAMC", "MILK", "QEUH"]:
         row["pillar_2"] = True
@@ -151,7 +156,7 @@ def main():
 
     date_dict = load_updated_dates(args.updated_date_file)
     accession_dict = load_accession(args.accession_file, log_handle)
-    new_columns = ["sample_date", "pillar_2", "sequence_name", "covv_accession_id", "edin_epi_week", "edin_epi_day", "why_excluded"]
+    new_columns = ["sample_date", "source_id", "pillar_2", "sequence_name", "covv_accession_id", "edin_epi_week", "edin_epi_day", "why_excluded"]
 
     with open(args.in_metadata, 'r', newline = '') as csv_in, \
          open(args.out_metadata, 'w', newline = '') as csv_out:
@@ -163,6 +168,7 @@ def main():
         for row in reader:
             try:
                 add_sample_date(row, date_dict)
+                add_source_id(row)
                 add_pillar_2(row)
                 add_sequence_name(row)
                 add_covv_accession_id(row, accession_dict)

modules/deduplicate_cog_uk.nf

@@ -148,43 +148,37 @@ process uk_unify_headers {
 }
 
 
-process uk_remove_duplicates_biosamplesourceid_by_date {
+process uk_label_sourceid_duplicates_to_omit {
     /**
-    * Where duplicate biosample_source_id, keeps the earliest
+    * Where duplicate source_id, labels all but the earliest as duplicates
     * @input uk_fasta, uk_metadata
     * @output uk_fasta_updated, uk_metadata_updated
     */
 
     publishDir "${publish_dev}/", pattern: "*.log", mode: 'copy'
 
     input:
-    path uk_fasta
     path uk_metadata
 
     output:
-    path "${uk_fasta.baseName}.deduplicated_by_biosamplesourceid.fa", emit: uk_fasta_updated
-    path "${uk_metadata.baseName}.deduplicated_by_biosamplesourceid.csv", emit: uk_metadata_updated
-    path "deduplicated_by_biosamplesourceid.log", emit: deduplicate_log
+    path "${uk_metadata.baseName}.deduplicated_by_sourceid.csv", emit: uk_metadata_updated
+    path "deduplicated_by_sourceid.log", emit: deduplicate_log
 
     script:
     """
     #!/usr/bin/env python3
     from Bio import SeqIO
     import csv
 
-    alignment = SeqIO.index("${uk_fasta}", "fasta")
-
     dup_dict = {}
     tokeep = set()
 
     with open("${uk_metadata}", 'r', newline = '') as csv_in:
         reader = csv.DictReader(csv_in, delimiter=",", quotechar='\"', dialect = "unix")
 
         for row in reader:
-            if row["why_excluded"]:
-                continue
             fasta_header = row["sequence_name"]
-            id = row["biosample_source_id"]
+            id = row["source_id"]
             epi_day = int(row["edin_epi_day"])
             completeness = float(row["unmapped_genome_completeness"])
 
@@ -200,7 +194,7 @@ process uk_remove_duplicates_biosamplesourceid_by_date {
             else:
                 dup_dict[id] = [{"fasta_header": fasta_header, "epi_day": epi_day, "completeness":completeness}]
 
-    with open("deduplicated_by_biosamplesourceid.log", "w") as log:
+    with open("deduplicated_by_sourceid.log", "w") as log:
         for k,v in dup_dict.items():
             tokeep.add(v[0]["fasta_header"])
             if len(v) > 1:
@@ -211,145 +205,39 @@ process uk_remove_duplicates_biosamplesourceid_by_date {
 
 
     with open("${uk_metadata}", 'r', newline = '') as csv_in, \
-         open("${uk_metadata.baseName}.deduplicated_by_biosamplesourceid.csv", 'w', newline = '') as csv_out, \
-         open("${uk_fasta.baseName}.deduplicated_by_biosamplesourceid.fa", 'w') as fasta_out:
+         open("${uk_metadata.baseName}.deduplicated_by_sourceid.csv", 'w', newline = '') as csv_out:
 
         reader = csv.DictReader(csv_in, delimiter=",", quotechar='\"', dialect = "unix")
-        writer = csv.DictWriter(csv_out, fieldnames = reader.fieldnames, delimiter=",", quotechar='\"', quoting=csv.QUOTE_MINIMAL, dialect = "unix")
+        writer = csv.DictWriter(csv_out, fieldnames = reader.fieldnames + ["duplicate"], delimiter=",", quotechar='\"', quoting=csv.QUOTE_MINIMAL, dialect = "unix")
         writer.writeheader()
 
         for row in reader:
-            if row["why_excluded"]:
-                writer.writerow(row)
-                continue
+            row["duplicate"] = None
             fasta_header = row["sequence_name"]
-            if fasta_header in tokeep:
-                writer.writerow(row)
-                seqrec = alignment[fasta_header]
-                fasta_out.write(">" + seqrec.id + "\\n")
-                fasta_out.write(str(seqrec.seq) + "\\n")
-            else:
-                row["why_excluded"] = "duplicate biosample_source_id"
-                writer.writerow(row)
+            if fasta_header not in tokeep:
+                row["duplicate"] = "True"
+            writer.writerow(row)
     """
 }
 
-process uk_remove_duplicates_rootbiosample_by_date {
-    /**
-    * Where duplicate root_biosample, keeps the oldest
-    * @input uk_fasta, uk_metadata
-    * @output uk_fasta_updated, uk_metadata_updated
-    * @params date
-    */
-
-    publishDir "${publish_dev}/", pattern: "*.log", mode: 'copy'
-
-    input:
-    path uk_fasta
-    path uk_metadata
-
-    output:
-    path "${uk_fasta.baseName}.deduplicated_by_rootbiosample.fa", emit: uk_fasta_updated
-    path "${uk_metadata.baseName}.deduplicated_by_rootbiosample.csv", emit: uk_metadata_updated
-    path "deduplicated_by_rootbiosample.log", emit: deduplicate_log
-
-    script:
-    """
-    #!/usr/bin/env python3
-    from Bio import SeqIO
-    import csv
 
-    alignment = SeqIO.index("${uk_fasta}", "fasta")
-
-    dup_dict = {}
-    tokeep = set()
-
-    with open("${uk_metadata}", 'r', newline = '') as csv_in:
-        reader = csv.DictReader(csv_in, delimiter=",", quotechar='\"', dialect = "unix")
-
-        for row in reader:
-            if row["why_excluded"]:
-                continue
-            fasta_header = row["sequence_name"]
-            id = row["root_biosample_source_id"]
-            epi_day = int(row["edin_epi_day"])
-            completeness = float(row["unmapped_genome_completeness"])
-
-            if id in ["None", "", None]:
-                tokeep.add(fasta_header)
-                continue
-
-            if id in dup_dict:
-                if epi_day < dup_dict[id][0]["epi_day"]:
-                    dup_dict[id].insert(0, {"fasta_header": fasta_header, "epi_day": epi_day, "completeness":completeness})
-                else:
-                    dup_dict[id].append({"fasta_header": fasta_header, "epi_day": epi_day, "completeness":completeness})
-            else:
-                dup_dict[id] = [{"fasta_header": fasta_header, "epi_day": epi_day, "completeness":completeness}]
-
-        with open("deduplicated_by_rootbiosample.log", "w") as log:
-            for k,v in dup_dict.items():
-                tokeep.add(v[0]["fasta_header"])
-                if len(v) > 1:
-                    for dup in v[1:]:
-                        log.write("For id %s, %s epi_day:%s completeness:%s kept, %s epi_day:%s completeness:%s removed as duplicate\\n" \
-                        %(k, v[0]["fasta_header"], v[0]["epi_day"], v[0]["completeness"], dup["fasta_header"], \
-                                    dup["epi_day"], dup["completeness"]))
-
-    with open("${uk_metadata}", 'r', newline = '') as csv_in, \
-         open("${uk_metadata.baseName}.deduplicated_by_rootbiosample.csv", 'w', newline = '') as csv_out, \
-         open("${uk_fasta.baseName}.deduplicated_by_rootbiosample.fa", 'w') as fasta_out:
-
-        reader = csv.DictReader(csv_in, delimiter=",", quotechar='\"', dialect = "unix")
-        writer = csv.DictWriter(csv_out, fieldnames = reader.fieldnames, delimiter=",", quotechar='\"', quoting=csv.QUOTE_MINIMAL, dialect = "unix")
-        writer.writeheader()
-
-        for row in reader:
-            if row["why_excluded"]:
-                writer.writerow(row)
-                continue
-            fasta_header = row["sequence_name"]
-            if fasta_header in tokeep:
-                writer.writerow(row)
-                seqrec = alignment[fasta_header]
-                fasta_out.write(">" + seqrec.id + "\\n")
-                fasta_out.write(str(seqrec.seq) + "\\n")
-            else:
-                row["why_excluded"] = "duplicate root_biosample_source_id"
-                writer.writerow(row)
-    """
-}
-
-
-workflow deduplicate_by_cogid_cog_uk {
+workflow deduplicate_cog_uk {
     take:
         uk_fasta
         uk_metadata
     main:
         uk_annotate_with_unmapped_genome_completeness(uk_fasta, uk_metadata)
         uk_remove_duplicates_COGID_by_proportionN(uk_fasta, uk_annotate_with_unmapped_genome_completeness.out)
         uk_unify_headers(uk_remove_duplicates_COGID_by_proportionN.out.uk_fasta_updated, uk_remove_duplicates_COGID_by_proportionN.out.uk_metadata_updated)
+        uk_label_sourceid_duplicates_to_omit(uk_remove_duplicates_COGID_by_proportionN.out.uk_metadata_updated)
     emit:
         fasta = uk_unify_headers.out
-        metadata = uk_remove_duplicates_COGID_by_proportionN.out.uk_metadata_updated
-}
-
-workflow deduplicate_by_biosample_cog_uk {
-    take:
-        uk_fasta
-        uk_metadata
-    main:
-        uk_remove_duplicates_biosamplesourceid_by_date(uk_fasta, uk_metadata)
-        uk_remove_duplicates_rootbiosample_by_date(uk_remove_duplicates_biosamplesourceid_by_date.out.uk_fasta_updated, uk_remove_duplicates_biosamplesourceid_by_date.out.uk_metadata_updated)
-    emit:
-        fasta = uk_remove_duplicates_rootbiosample_by_date.out.uk_fasta_updated
-        metadata = uk_remove_duplicates_rootbiosample_by_date.out.uk_metadata_updated
+        metadata = uk_label_sourceid_duplicates_to_omit.out.uk_metadata_updated
 }
 
 
 workflow {
     uk_fasta = file(params.uk_fasta)
     uk_metadata = file(params.uk_metadata)
-    deduplicate_by_cogid_cog_uk(uk_fasta, uk_metadata)
-    deduplicate_by_biosample_cog_uk(deduplicate_by_cogid_cog_uk.out.fasta, deduplicate_by_cogid_cog_uk.out.metadata)
+    deduplicate_cog_uk(uk_fasta, uk_metadata)
 }

modules/publish_all.nf

@@ -33,7 +33,7 @@ process combine_cog_gisaid {
           --filter-column fasta_header covv_accession_id central_sample_id biosample_source_id secondary_identifier root_sample_id \
                           pillar_2 \
                           sequence_name sample_date epi_week \
-                          country adm1 adm2 outer_postcode adm2_raw adm2_source nuts1 region latitude longitude location \
+                          country adm1 adm2 outer_postcode adm2_raw adm2_source NUTS1 region latitude longitude location \
                           submission_org_code is_surveillance is_community is_hcw \
                           is_travel_history travel_history \
                           lineage lineage_support lineages_version \
@@ -160,7 +160,8 @@ process uk_geography {
       --in-metadata ${uk_metadata} \
       --index-column sequence_name \
       --filter-column central_sample_id sequence_name sample_date epi_week \
-                      adm0 adm1 adm2 adm2_private \
+                      adm0 adm1 adm2 adm2_private adm1_UK \
+      --where-column adm1_UK=adm1 \
       --out-fasta geography_tmp/fetch.fa \
       --out-metadata geography_tmp/fetch.csv \
       --restrict

resources/publish_recipes.json

@@ -24,7 +24,7 @@
       "suffix": "public",
       "data": "cog_global",
       "metadata_fields": ["sequence_name","cog_id","gisaid_id","sample_date","epi_week","country","adm1","pillar_2","is_surveillance","is_travel_history","travel_history","lineage","lineage_support","lineages_version"],
-      "where": "gisaid_id=covv_accession_id cog_id=central_sample_id"
+      "where": "gisaid_id=covv_accession_id cog_id=central_sample_id adm1=adm1_UK"
     },
     {
       "suffix": "consortium",

workflows/process_cog_uk.nf

@@ -19,12 +19,11 @@ workflow process_cog_uk {
     main:
       preprocess_cog_uk(uk_fasta, uk_metadata, uk_accessions)
       pangolin_cog_uk(preprocess_cog_uk.out.fasta, preprocess_cog_uk.out.metadata)
-      deduplicate_by_cogid_cog_uk(preprocess_cog_uk.out.fasta, pangolin_cog_uk.out.metadata)
-      align_and_variant_call_cog_uk(deduplicate_by_cogid_cog_uk.out.fasta)
-      deduplicate_by_biosample_cog_uk(align_and_variant_call_cog_uk.out.fasta,deduplicate_by_cogid_cog_uk.out.metadata)
-      filter_and_trim_cog_uk(deduplicate_by_biosample_cog_uk.out.fasta, deduplicate_by_biosample_cog_uk.out.metadata)
+      deduplicate_cog_uk(preprocess_cog_uk.out.fasta, pangolin_cog_uk.out.metadata)
+      align_and_variant_call_cog_uk(deduplicate_cog_uk.out.fasta)
+      filter_and_trim_cog_uk(deduplicate_cog_uk.out.fasta, deduplicate_cog_uk.out.metadata)
     emit:
-      unaligned_fasta = deduplicate_by_cogid_cog_uk.out.fasta
+      unaligned_fasta = deduplicate_cog_uk.out.fasta
       aligned_fasta = align_and_variant_call_cog_uk.out.fasta
       trimmed_fasta = filter_and_trim_cog_uk.out.fasta
       metadata = filter_and_trim_cog_uk.out.metadata