Skip to content

Commit 4cceb7e

Browse files
tives82tives82
authored
Merge pull request #7 from UPHL-BioNGS/Tom_dev_240126
Update to Nextclade v3 and Enhancement of Coverage and Depth Analysis
2 parents 032d370 + aaf2ffe commit 4cceb7e

21 files changed

Lines changed: 522 additions & 50 deletions

‎README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,8 @@ The pipeline is built using [Nextflow](https://www.nextflow.io), a workflow tool
4949
> **Clean read data undergo assembly and influenza typing and subtyping. Based on the subtype information, Nextclade variables are gathered.**
5050
5151
* Assembly of influenza gene segments with (`IRMA`) using the built-in FLU module. Also, influenza typing and H/N subtype classifications are made.
52+
* Calculate the reference length, sequence length, and percent_coverage for segments assembled by IRMA with (`IRMA_SEGMENT_COVERAGE`)
53+
* Calculate the number of mapped reads and mean depth for segments assembled by IRMA with (`SAMTOOLS_MAPPED_READS`)
5254
* QC of consensus assembly (`IRMA_Consensus_QC`).
5355
* Generate IRMA consensus QC report (`IRMA_Consensus_QC_Reportsheet`)
5456
* Annotation of IRMA consensus sequences with (`VADR`)
@@ -69,6 +71,7 @@ The pipeline is built using [Nextflow](https://www.nextflow.io), a workflow tool
6971
> **Compiles report sheets from modules and outputs a pipeline summary report tsv file.**
7072
7173
* The (`Summary_Report`) consolidates and merges multiple report sheets into a single comprehensive summary report.
74+
* The (`merged_bam_coverage_results`) merges the gene segment report sheets detailing mapped reads, mean depth, reference length, sequence length, and percent_coverage.
7275

7376
## Quick Start
7477

‎bin/calc_percent_cov.py‎

Lines changed: 149 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,149 @@
1+
#!/usr/bin/env python
2+
3+
# version = '1.0.0'
4+
5+
# Modified from script https://github.com/CDPHE-bioinformatics/CDPHE-influenza/blob/main/scripts/calc_percent_cov.py
6+
7+
# import python modules
8+
import pandas as pd
9+
from datetime import date
10+
from Bio import SeqIO
11+
from Bio.SeqRecord import SeqRecord
12+
13+
import sys
14+
import argparse
15+
import subprocess
16+
17+
### Segment length dictionary
18+
ref_len_dict = {
19+
"A_MP": 982,
20+
"A_NP": 1497,
21+
"A_NS": 863,
22+
"A_PA": 2151,
23+
"A_PB1": 2274,
24+
"A_PB2": 2280,
25+
"A_HA_H1": 1704,
26+
"A_HA_H10": 1686,
27+
"A_HA_H11": 1698,
28+
"A_HA_H12": 1695,
29+
"A_HA_H13": 1701,
30+
"A_HA_H14": 1707,
31+
"A_HA_H15": 1713,
32+
"A_HA_H16": 1698,
33+
"A_HA_H2": 1689,
34+
"A_HA_H3": 1704,
35+
"A_HA_H4": 1695,
36+
"A_HA_H5": 1707,
37+
"A_HA_H6": 1704,
38+
"A_HA_H7": 1713,
39+
"A_HA_H8": 1701,
40+
"A_HA_H9": 1683,
41+
"A_NA_N1": 1413,
42+
"A_NA_N2": 1410,
43+
"A_NA_N3": 1410,
44+
"A_NA_N4": 1413,
45+
"A_NA_N5": 1422,
46+
"A_NA_N6": 1413,
47+
"A_NA_N7": 1416,
48+
"A_NA_N8": 1413,
49+
"A_NA_N9": 1413,
50+
"B_HA": 1758,
51+
"B_MP": 1139,
52+
"B_NA": 1408,
53+
"B_NP": 1683,
54+
"B_NS": 1034,
55+
"B_PA": 2181,
56+
"B_PB1": 2263,
57+
"B_PB2": 2313,
58+
}
59+
60+
61+
#### FUNCTIONS #####
62+
def getOptions():
63+
parser = argparse.ArgumentParser(description="Parses command.")
64+
parser.add_argument("fasta_files", help="Path to the fasta file")
65+
parser.add_argument("meta_id", help="Meta ID")
66+
options = parser.parse_args()
67+
return options
68+
69+
70+
def get_fasta_file_basename(fasta_file_path):
71+
basename = fasta_file_path.split("/")[-1] # strip directories
72+
return basename
73+
74+
75+
def get_segment_name(fasta_file_path):
76+
basename = fasta_file_path.split("/")[-1] # strip directories
77+
segment_name = basename.split(".")[0] # remove file extension
78+
return segment_name
79+
80+
81+
def get_gene_name(fasta_file_path):
82+
basename = fasta_file_path.split("/")[-1] # strip directories
83+
segment_name = basename.split(".")[0] # remove file extension
84+
gene_name = segment_name.split("_")[1] # extract gene name
85+
return gene_name
86+
87+
88+
def get_seq_length(fasta_file_path):
89+
# read in fasta file
90+
record = SeqIO.read(fasta_file_path, "fasta")
91+
92+
# get length of non ambigous bases
93+
seq = record.seq
94+
seq_length = seq.count("A") + seq.count("C") + seq.count("G") + seq.count("T")
95+
96+
return seq_length
97+
98+
99+
def calc_percent_cov(seq_length, ref_len_dict, segment_name):
100+
# calcuate per cov based on expected ref length
101+
expected_length = ref_len_dict[segment_name]
102+
percent_coverage = round(((seq_length / expected_length) * 100), 2)
103+
104+
return percent_coverage
105+
106+
107+
def create_output(meta_id, segment_name, seq_length, percent_coverage, reference_length):
108+
df = pd.DataFrame()
109+
df["Sample"] = [meta_id]
110+
df["segment_name"] = [segment_name]
111+
df["reference_length"] = [reference_length]
112+
df["seq_length"] = [seq_length]
113+
df["percent_coverage"] = [percent_coverage]
114+
115+
# Construct the output filename header
116+
output_header = "Sample\tsegment_name\treference_length\tseq_length\tpercent_coverage"
117+
118+
# Construct the output filename
119+
output_filename = f"{meta_id}.{segment_name}.perc_cov_results.tsv"
120+
121+
# Write the dataframe to the output file
122+
with open(output_filename, "w") as f:
123+
f.write(output_header + "\n")
124+
df.to_csv(f, sep="\t", index=False, header=False)
125+
126+
127+
#### MAIN ####
128+
if __name__ == "__main__":
129+
options = getOptions()
130+
fasta_file_path = options.fasta_files
131+
meta_id = options.meta_id
132+
133+
basename = get_fasta_file_basename(fasta_file_path=fasta_file_path)
134+
135+
segment_name = get_segment_name(fasta_file_path=fasta_file_path)
136+
gene_name = get_gene_name(fasta_file_path=fasta_file_path)
137+
138+
seq_length = get_seq_length(fasta_file_path=fasta_file_path)
139+
percent_coverage = calc_percent_cov(seq_length=seq_length, ref_len_dict=ref_len_dict, segment_name=segment_name)
140+
141+
reference_length = ref_len_dict[segment_name]
142+
143+
create_output(
144+
meta_id=meta_id,
145+
segment_name=segment_name,
146+
seq_length=seq_length,
147+
percent_coverage=percent_coverage,
148+
reference_length=reference_length,
149+
)

‎bin/flu_nextclade_variables.py‎

Lines changed: 13 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -4,31 +4,13 @@
44
import os
55
from os.path import exists
66
import argparse
7-
import pandas as pd
87

9-
# Dictionary containing data for various flu subtypes.
10-
# Each subtype has associated Nextclade dataset, reference, and tag variables.
8+
# Dictionary containing data for various flu subtypes with their associated datasets.
119
flu_subtypes = {
12-
"H1N1": {
13-
"dataset": "flu_h1n1pdm_ha",
14-
"reference": "CY121680",
15-
"tag": "2023-08-10T12:00:00Z",
16-
},
17-
"H3N2": {
18-
"dataset": "flu_h3n2_ha",
19-
"reference": "CY163680",
20-
"tag": "2023-08-10T12:00:00Z",
21-
},
22-
"Victoria": {
23-
"dataset": "flu_vic_ha",
24-
"reference": "KX058884",
25-
"tag": "2023-08-10T12:00:00Z",
26-
},
27-
"Yamagata": {
28-
"dataset": "flu_yam_ha",
29-
"reference": "JN993010",
30-
"tag": "2022-07-27T12:00:00Z",
31-
},
10+
"H1N1": {"dataset": "flu_h1n1pdm_ha"},
11+
"H3N2": {"dataset": "flu_h3n2_ha"},
12+
"Victoria": {"dataset": "flu_vic_ha"},
13+
"Yamagata": {"dataset": "flu_yam_ha"},
3214
}
3315

3416

@@ -57,12 +39,14 @@ def main():
5739
print(f"Error: Invalid flu subtype '{flu_subtype}' for sample '{args.sample}'")
5840
return
5941

60-
# For each variables (dataset, reference, tag) of the identified subtype, write it to a separate file and print variables.
61-
for item in ["dataset", "reference", "tag"]:
62-
file_path = flu_subtypes[flu_subtype][item]
63-
with open(file_path, "w") as f:
64-
f.write(f"{flu_subtypes[flu_subtype][item]}\n")
65-
print(f" {item}: {flu_subtypes[flu_subtype][item]} (output to {file_path})")
42+
# Prepare the dataset and file_path
43+
dataset = flu_subtypes[flu_subtype]["dataset"]
44+
file_path = dataset
45+
46+
# Write to the file and print information
47+
with open(file_path, "w") as f:
48+
f.write(f"{dataset}\n")
49+
print(f" {dataset}: {dataset} (output to {file_path})")
6650

6751

6852
if __name__ == "__main__":

‎bin/irma_consensus_qc.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -86,13 +86,13 @@ def main(consensus_fasta, meta_id):
8686
writer.writerow(
8787
[
8888
"Sample",
89-
"IRMA consensus ACTG Count",
90-
"IRMA consensus Degenerate Count",
91-
"IRMA consensus N Count",
92-
"IRMA consensus Total Count",
93-
"IRMA consensus Segment Count",
94-
"IRMA consensus N50",
95-
"IRMA consensus GC Content",
89+
"IRMA_consensus_ACTG_count",
90+
"IRMA_consensus_degenerate_count",
91+
"IRMA_consensus_N_count",
92+
"IRMA_consensus_total_count",
93+
"IRMA_consensus_segment_count",
94+
"IRMA_consensus_N50",
95+
"IRMA_consensus_GC_content",
9696
]
9797
)
9898

‎bin/merge_bam_coverage.py‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
import pandas as pd
2+
import sys
3+
4+
5+
def main():
6+
if len(sys.argv) != 4:
7+
print("Usage: python script.py <bam_results> <coverage_results> <output_file>")
8+
sys.exit(1)
9+
10+
bam_file_path = sys.argv[1]
11+
coverage_file_path = sys.argv[2]
12+
output_file_path = sys.argv[3]
13+
14+
# Read the TSV files into DataFrames
15+
bam_df = pd.read_csv(bam_file_path, sep="\t")
16+
coverage_df = pd.read_csv(coverage_file_path, sep="\t")
17+
18+
# Make sure 'Sample' is the first column in both DataFrames
19+
bam_df = bam_df[["Sample"] + [col for col in bam_df.columns if col != "Sample"]]
20+
coverage_df = coverage_df[["Sample"] + [col for col in coverage_df.columns if col != "Sample"]]
21+
22+
# Merge the DataFrames on the 'Sample' column
23+
merged_df = pd.merge(bam_df, coverage_df, on="Sample", how="outer")
24+
25+
# Round all numeric columns to 2 decimal places
26+
merged_df = merged_df.round(2)
27+
28+
# Make sure 'Sample' is the first column in the merged DataFrame
29+
merged_df = merged_df[["Sample"] + [col for col in merged_df.columns if col != "Sample"]]
30+
31+
# Write the merged DataFrame to a TSV file
32+
merged_df.to_csv(output_file_path, sep="\t", index=False)
33+
34+
35+
if __name__ == "__main__":
36+
main()

‎bin/merge_reports.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,9 @@ def merge_tsvs(files):
1111
temp_df = pd.read_csv(file, sep="\t")
1212
df = pd.merge(df, temp_df, on="Sample", how="outer")
1313

14+
# Round all numeric columns to 2 decimal places
15+
df = df.round(2)
16+
1417
return df
1518

1619

‎conf/modules.config‎

Lines changed: 52 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -290,13 +290,63 @@ process {
290290
pattern: "*"
291291
]
292292
}
293+
withName: 'IRMA_SEGMENT_COVERAGE' {
294+
ext.args = { "" }
295+
ext.when = { }
296+
publishDir = [
297+
enabled: true,
298+
mode: "${params.publish_dir_mode}",
299+
path: { "${params.outdir}/irma_segment_coverage/${meta.id}" },
300+
pattern: "*"
301+
]
302+
}
303+
withName: 'MERGE_COVERAGE_RESULTS' {
304+
ext.args = { "" }
305+
ext.when = { }
306+
publishDir = [
307+
enabled: true,
308+
mode: "${params.publish_dir_mode}",
309+
path: { "${params.outdir}/irma_segment_coverage/" },
310+
pattern: "*"
311+
]
312+
}
313+
withName: 'SAMTOOLS_MAPPED_READS' {
314+
ext.args = { "" }
315+
ext.when = { }
316+
publishDir = [
317+
enabled: true,
318+
mode: "${params.publish_dir_mode}",
319+
path: { "${params.outdir}/samtools_mapped_reads/${meta.id}" },
320+
pattern: "*"
321+
]
322+
}
323+
withName: 'MERGE_BAM_RESULTS' {
324+
ext.args = { "" }
325+
ext.when = { }
326+
publishDir = [
327+
enabled: true,
328+
mode: "${params.publish_dir_mode}",
329+
path: { "${params.outdir}/samtools_mapped_reads/" },
330+
pattern: "*"
331+
]
332+
}
333+
withName: 'MERGE_BAM_COVERAGE_RESULTS' {
334+
ext.args = { "" }
335+
ext.when = { }
336+
publishDir = [
337+
enabled: true,
338+
mode: "${params.publish_dir_mode}",
339+
path: { "${params.outdir}/SUMMARY_REPORTS" },
340+
pattern: "*"
341+
]
342+
}
293343
withName: COMBINED_SUMMARY_REPORT {
294344
ext.args = { "" }
295345
ext.when = { }
296346
publishDir = [
297347
enabled: true,
298348
mode: "${params.publish_dir_mode}",
299-
path: { "${params.outdir}/SUMMARY_REPORT" },
349+
path: { "${params.outdir}/SUMMARY_REPORTS" },
300350
pattern: "*"
301351
]
302352
}
@@ -306,7 +356,7 @@ process {
306356
publishDir = [
307357
enabled: true,
308358
mode: "${params.publish_dir_mode}",
309-
path: { "${params.outdir}/SUMMARY_REPORT" },
359+
path: { "${params.outdir}/SUMMARY_REPORTS" },
310360
pattern: "*"
311361
]
312362
}

‎modules/local/irma.nf‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,16 @@ process IRMA {
22
tag "$meta.id"
33
label 'process_high'
44

5-
container 'quay.io/staphb/irma:1.1.3'
5+
container 'quay.io/staphb/irma:1.1.4'
66

77
input:
88
tuple val(meta), path(reads)
99
val(irma_module)
1010

1111
output:
1212
tuple val(meta), path("${meta.id}/") , emit: irma
13+
tuple val(meta), path("${meta.id}/*.bam") , emit: irma_bam
14+
tuple val(meta), path("${meta.id}/*.fasta") , emit: irma_fasta
1315
tuple val(meta), path("*.irma.consensus.fasta") , optional:true, emit: assembly
1416
tuple val(meta), path("*_LOW_ABUNDANCE.txt") , optional:true, emit: failed_assembly
1517
tuple val(meta), path("*_HA.fasta") , optional:true, emit: HA
@@ -110,3 +112,4 @@ process IRMA {
110112
"""
111113
}
112114

115+

0 commit comments

Comments
 (0)