|
28
|
1 #!/usr/bin/env python
|
|
|
2
|
|
0
|
3 import argparse
|
|
|
4 import os
|
|
|
5 import pandas
|
|
|
6 import pypandoc
|
|
|
7 import re
|
|
|
8 import subprocess
|
|
|
9 import sys
|
|
|
10
|
|
|
11 from Bio import SeqIO
|
|
|
12 from datetime import date
|
|
|
13 from mdutils.mdutils import MdUtils
|
|
14
|
14 # FIXME: TableOfContents doesn't work.
|
|
|
15 # from mdutils.tools import TableOfContents
|
|
0
|
16
|
|
|
17 CDC_ADVISORY = 'The analysis and report presented here should be treated as preliminary. Please contact the CDC/BDRD with any results regarding _Bacillus anthracis_.'
|
|
|
18
|
|
|
19
|
|
|
20 class PimaReport:
|
|
|
21
|
|
21
|
22 def __init__(self, analysis_name=None, amr_deletions_file=None, amr_matrix_files=None, assembler_version=None,
|
|
|
23 assembly_fasta_file=None, assembly_name=None, bedtools_version=None, blastn_version=None,
|
|
|
24 circos_files=None, compute_sequence_length_file=None, contig_coverage_file=None, dbkey=None,
|
|
26
|
25 dnadiff_snps_file=None, dnadiff_version=None, errors_file=None, feature_bed_files=None,
|
|
21
|
26 feature_png_files=None, flye_assembly_info_file=None, genome_insertions_file=None, gzipped=None,
|
|
26
|
27 illumina_forward_read_file=None, illumina_reverse_read_file=None, kraken2_report_file=None,
|
|
28
|
28 kraken2_version=None, lrn_risk_amr_file=None, lrn_risk_blacklist_file=None, lrn_risk_vf_file=None,
|
|
|
29 minimap2_version=None, mutation_regions_bed_file=None, mutation_regions_tsv_files=None,
|
|
|
30 ont_file=None, pima_css=None, plasmids_file=None, quast_report_file=None, read_type=None,
|
|
|
31 reference_insertions_file=None, samtools_version=None, varscan_version=None):
|
|
0
|
32 self.ofh = open("process_log.txt", "w")
|
|
|
33
|
|
1
|
34 self.ofh.write("amr_deletions_file: %s\n" % str(amr_deletions_file))
|
|
|
35 self.ofh.write("amr_matrix_files: %s\n" % str(amr_matrix_files))
|
|
0
|
36 self.ofh.write("analysis_name: %s\n" % str(analysis_name))
|
|
21
|
37 self.ofh.write("assembler_version: %s\n" % str(assembler_version))
|
|
0
|
38 self.ofh.write("assembly_fasta_file: %s\n" % str(assembly_fasta_file))
|
|
|
39 self.ofh.write("assembly_name: %s\n" % str(assembly_name))
|
|
12
|
40 self.ofh.write("bedtools_version: %s\n" % str(bedtools_version))
|
|
2
|
41 self.ofh.write("blastn_version: %s\n" % str(blastn_version))
|
|
13
|
42 self.ofh.write("circos_files: %s\n" % str(circos_files))
|
|
1
|
43 self.ofh.write("compute_sequence_length_file: %s\n" % str(compute_sequence_length_file))
|
|
|
44 self.ofh.write("contig_coverage_file: %s\n" % str(contig_coverage_file))
|
|
|
45 self.ofh.write("dbkey: %s\n" % str(dbkey))
|
|
|
46 self.ofh.write("dnadiff_snps_file: %s\n" % str(dnadiff_snps_file))
|
|
2
|
47 self.ofh.write("dnadiff_version: %s\n" % str(dnadiff_version))
|
|
18
|
48 self.ofh.write("errors_file: %s\n" % str(errors_file))
|
|
0
|
49 self.ofh.write("feature_bed_files: %s\n" % str(feature_bed_files))
|
|
|
50 self.ofh.write("feature_png_files: %s\n" % str(feature_png_files))
|
|
1
|
51 self.ofh.write("flye_assembly_info_file: %s\n" % str(flye_assembly_info_file))
|
|
0
|
52 self.ofh.write("gzipped: %s\n" % str(gzipped))
|
|
1
|
53 self.ofh.write("genome_insertions_file: %s\n" % str(genome_insertions_file))
|
|
26
|
54 self.ofh.write("illumina_forward_read_file: %s\n" % str(illumina_forward_read_file))
|
|
|
55 self.ofh.write("illumina_reverse_read_file: %s\n" % str(illumina_reverse_read_file))
|
|
2
|
56 self.ofh.write("kraken2_report_file: %s\n" % str(kraken2_report_file))
|
|
|
57 self.ofh.write("kraken2_version: %s\n" % str(kraken2_version))
|
|
28
|
58 self.ofh.write("lrn_risk_amr_file: %s\n" % str(lrn_risk_amr_file))
|
|
|
59 self.ofh.write("lrn_risk_blacklist_file: %s\n" % str(lrn_risk_blacklist_file))
|
|
|
60 self.ofh.write("lrn_risk_vf_file: %s\n" % str(lrn_risk_vf_file))
|
|
12
|
61 self.ofh.write("minimap2_version: %s\n" % str(minimap2_version))
|
|
0
|
62 self.ofh.write("mutation_regions_bed_file: %s\n" % str(mutation_regions_bed_file))
|
|
|
63 self.ofh.write("mutation_regions_tsv_files: %s\n" % str(mutation_regions_tsv_files))
|
|
26
|
64 self.ofh.write("ont_file: %s\n" % str(ont_file))
|
|
0
|
65 self.ofh.write("pima_css: %s\n" % str(pima_css))
|
|
1
|
66 self.ofh.write("plasmids_file: %s\n" % str(plasmids_file))
|
|
13
|
67 self.ofh.write("quast_report_file: %s\n" % str(quast_report_file))
|
|
18
|
68 self.ofh.write("read_type: %s\n" % str(read_type))
|
|
1
|
69 self.ofh.write("reference_insertions_file: %s\n" % str(reference_insertions_file))
|
|
12
|
70 self.ofh.write("samtools_version: %s\n" % str(samtools_version))
|
|
|
71 self.ofh.write("varscan_version: %s\n" % str(varscan_version))
|
|
0
|
72
|
|
|
73 # General
|
|
|
74 self.doc = None
|
|
|
75 self.report_md = 'pima_report.md'
|
|
|
76
|
|
|
77 # Inputs
|
|
1
|
78 self.amr_deletions_file = amr_deletions_file
|
|
|
79 self.amr_matrix_files = amr_matrix_files
|
|
13
|
80 self.analysis_name = analysis_name.split('_')[0]
|
|
|
81 self.ofh.write("self.analysis_name: %s\n" % str(self.analysis_name))
|
|
21
|
82 if assembler_version is None:
|
|
|
83 self.assembler_version = 'assembler (version unknown)'
|
|
|
84 else:
|
|
|
85 if read_type == 'ont':
|
|
|
86 # Assembler is flye.
|
|
|
87 assembler_version = assembler_version.rstrip(' _assembly info_')
|
|
|
88 else:
|
|
|
89 # Assembler is spades.
|
|
|
90 assembler_version = assembler_version.rstrip(' _contigs')
|
|
|
91 self.assembler_version = re.sub('_', '.', assembler_version)
|
|
0
|
92 self.assembly_fasta_file = assembly_fasta_file
|
|
12
|
93 self.assembly_name = re.sub('_', '.', assembly_name.rstrip(' _consensus_'))
|
|
|
94 if bedtools_version is None:
|
|
|
95 self.bedtools_version = 'bedtools (version unknown)'
|
|
|
96 else:
|
|
|
97 self.bedtools_version = re.sub('_', '.', bedtools_version.rstrip(' _genome insertions'))
|
|
|
98 if blastn_version is None:
|
|
|
99 self.blastn_version = 'blastn (version unknown)'
|
|
|
100 else:
|
|
|
101 self.blastn_version = re.sub('_', '.', blastn_version.rstrip(' _features_'))
|
|
13
|
102 self.circos_files = circos_files
|
|
1
|
103 self.compute_sequence_length_file = compute_sequence_length_file
|
|
|
104 self.contig_coverage_file = contig_coverage_file
|
|
|
105 self.dbkey = dbkey
|
|
|
106 self.dnadiff_snps_file = dnadiff_snps_file
|
|
12
|
107 if dnadiff_version is None:
|
|
|
108 self.dnadiff_version = 'dnadiff (version unknown)'
|
|
|
109 else:
|
|
|
110 self.dnadiff_version = re.sub('_', '.', dnadiff_version.rstrip(' _snps_'))
|
|
18
|
111 self.errors_file = errors_file
|
|
0
|
112 self.feature_bed_files = feature_bed_files
|
|
|
113 self.feature_png_files = feature_png_files
|
|
1
|
114 self.flye_assembly_info_file = flye_assembly_info_file
|
|
0
|
115 self.gzipped = gzipped
|
|
1
|
116 self.genome_insertions_file = genome_insertions_file
|
|
26
|
117 self.illumina_forward_read_file = illumina_forward_read_file
|
|
|
118 self.illumina_reverse_read_file = illumina_reverse_read_file
|
|
2
|
119 self.kraken2_report_file = kraken2_report_file
|
|
12
|
120 if kraken2_version is None:
|
|
|
121 self.kraken2_version = 'kraken2 (version unknown)'
|
|
|
122 else:
|
|
|
123 self.kraken2_version = re.sub('_', '.', kraken2_version.rstrip(' _report_'))
|
|
28
|
124 self.lrn_risk_amr_file = lrn_risk_amr_file
|
|
|
125 self.lrn_risk_blacklist_file = lrn_risk_blacklist_file
|
|
|
126 self.lrn_risk_vf_file = lrn_risk_vf_file
|
|
12
|
127 if minimap2_version is None:
|
|
|
128 self.minimap2_version = 'minimap2 (version unknown)'
|
|
|
129 else:
|
|
|
130 self.minimap2_version = re.sub('_', '.', minimap2_version)
|
|
0
|
131 self.mutation_regions_bed_file = mutation_regions_bed_file
|
|
|
132 self.mutation_regions_tsv_files = mutation_regions_tsv_files
|
|
|
133 self.pima_css = pima_css
|
|
1
|
134 self.plasmids_file = plasmids_file
|
|
13
|
135 self.quast_report_file = quast_report_file
|
|
18
|
136 self.read_type = read_type.upper()
|
|
13
|
137 self.reference_insertions_file = reference_insertions_file
|
|
1
|
138 self.reference_insertions_file = reference_insertions_file
|
|
12
|
139 if samtools_version is None:
|
|
|
140 self.samtools_version = 'samtools (version unknown)'
|
|
|
141 else:
|
|
|
142 self.samtools_version = re.sub('_', '.', samtools_version)
|
|
|
143 if varscan_version is None:
|
|
|
144 self.varscan_version = 'varscan (version unknown)'
|
|
|
145 else:
|
|
|
146 self.varscan_version = re.sub('_', '.', varscan_version)
|
|
0
|
147
|
|
|
148 # Titles
|
|
|
149 self.alignment_title = 'Comparison with reference'
|
|
|
150 self.alignment_notes_title = 'Alignment notes'
|
|
|
151 self.amr_matrix_title = 'AMR matrix'
|
|
|
152 self.assembly_methods_title = 'Assembly'
|
|
|
153 self.assembly_notes_title = 'Assembly notes'
|
|
|
154 self.basecalling_title = 'Basecalling'
|
|
|
155 self.basecalling_methods_title = 'Basecalling'
|
|
|
156 self.contamination_methods_title = 'Contamination check'
|
|
|
157 self.contig_alignment_title = 'Alignment vs. reference contigs'
|
|
|
158 self.feature_title = 'Features found in the assembly'
|
|
|
159 self.feature_methods_title = 'Feature annotation'
|
|
|
160 self.feature_plot_title = 'Feature annotation plots'
|
|
|
161 self.large_indel_title = 'Large insertions & deletions'
|
|
28
|
162 self.lrn_risk_title = 'LRNRisk isolate classification'
|
|
0
|
163 self.methods_title = 'Methods'
|
|
18
|
164 self.mutation_errors_title = 'Errors finding mutations in the sample'
|
|
0
|
165 self.mutation_title = 'Mutations found in the sample'
|
|
|
166 self.mutation_methods_title = 'Mutation screening'
|
|
|
167 self.plasmid_methods_title = 'Plasmid annotation'
|
|
1
|
168 self.plasmid_title = 'Plasmid annotation'
|
|
21
|
169 self.reference_genome_title = 'Reference genome'
|
|
0
|
170 self.reference_methods_title = 'Reference comparison'
|
|
|
171 self.snp_indel_title = 'SNPs and small indels'
|
|
13
|
172 self.summary_title = 'Summary'
|
|
0
|
173
|
|
|
174 # Methods
|
|
|
175 self.methods = pandas.Series(dtype='float64')
|
|
|
176 self.methods[self.contamination_methods_title] = pandas.Series(dtype='float64')
|
|
|
177 self.methods[self.assembly_methods_title] = pandas.Series(dtype='float64')
|
|
21
|
178 self.methods[self.reference_genome_title] = pandas.Series(dtype='float64')
|
|
0
|
179 self.methods[self.reference_methods_title] = pandas.Series(dtype='float64')
|
|
|
180 self.methods[self.mutation_methods_title] = pandas.Series(dtype='float64')
|
|
|
181 self.methods[self.feature_methods_title] = pandas.Series(dtype='float64')
|
|
|
182 self.methods[self.plasmid_methods_title] = pandas.Series(dtype='float64')
|
|
|
183
|
|
|
184 # Notes
|
|
|
185 self.assembly_notes = pandas.Series(dtype=object)
|
|
|
186 self.alignment_notes = pandas.Series(dtype=object)
|
|
|
187 self.contig_alignment = pandas.Series(dtype=object)
|
|
|
188
|
|
|
189 # Values
|
|
|
190 self.assembly_size = 0
|
|
|
191 self.contig_info = None
|
|
|
192 self.feature_hits = pandas.Series(dtype='float64')
|
|
26
|
193 self.ont_fast5 = None
|
|
|
194 self.ont_file = ont_file
|
|
|
195 self.ont_n50 = None
|
|
|
196 self.ont_read_count = None
|
|
14
|
197 # TODO: should the following be passed as a parameter?
|
|
|
198 self.ont_coverage_min = 30
|
|
|
199 # TODO: should the following be passed as a parameter?
|
|
1
|
200 self.ont_n50_min = 2500
|
|
26
|
201
|
|
21
|
202 if self.read_type == 'ONT':
|
|
|
203 self.ont_raw_fastq = self.analysis_name
|
|
26
|
204 self.ont_bases = 0
|
|
|
205 self.illumina_bases = None
|
|
21
|
206 self.illumina_fastq = None
|
|
26
|
207 self.illumina_length_mean = None
|
|
|
208 self.illumina_read_count = None
|
|
21
|
209 else:
|
|
26
|
210 self.illumina_fastq = self.analysis_name
|
|
|
211 self.illumina_bases = 0
|
|
|
212 self.illumina_length_mean = 0
|
|
|
213 self.illumina_read_count = 0
|
|
|
214 self.ont_bases = None
|
|
21
|
215 self.ont_raw_fastq = None
|
|
0
|
216
|
|
|
217 # Actions
|
|
|
218 self.did_guppy_ont_fast5 = False
|
|
|
219 self.did_qcat_ont_fastq = False
|
|
21
|
220 self.ofh.write("self.read_type: %s\n" % str(self.read_type))
|
|
|
221 if self.read_type == 'ONT':
|
|
26
|
222 self.info_ont_fastq(self.ont_file)
|
|
21
|
223 else:
|
|
26
|
224 self.info_illumina_fastq([self.illumina_forward_read_file, self.illumina_reverse_read_file])
|
|
0
|
225 self.load_contig_info()
|
|
|
226
|
|
|
227 def run_command(self, command):
|
|
|
228 self.ofh.write("\nXXXXXX In run_command, command:\n%s\n\n" % str(command))
|
|
|
229 try:
|
|
|
230 return re.split('\\n', subprocess.check_output(command, shell=True).decode('utf-8'))
|
|
|
231 except Exception:
|
|
|
232 message = 'Command %s failed: exiting...' % command
|
|
|
233 sys.exit(message)
|
|
|
234
|
|
|
235 def format_kmg(self, number, decimals=0):
|
|
|
236 self.ofh.write("\nXXXXXX In format_kmg, number:\n%s\n" % str(number))
|
|
|
237 self.ofh.write("XXXXXX In format_kmg, decimals:\n%s\n\n" % str(decimals))
|
|
|
238 if number == 0:
|
|
|
239 return '0'
|
|
|
240 magnitude_powers = [10**9, 10**6, 10**3, 1]
|
|
|
241 magnitude_units = ['G', 'M', 'K', '']
|
|
|
242 for i in range(len(magnitude_units)):
|
|
|
243 if number >= magnitude_powers[i]:
|
|
|
244 magnitude_power = magnitude_powers[i]
|
|
|
245 magnitude_unit = magnitude_units[i]
|
|
|
246 return ('{:0.' + str(decimals) + 'f}').format(number / magnitude_power) + magnitude_unit
|
|
|
247
|
|
|
248 def load_contig_info(self):
|
|
|
249 self.contig_info = pandas.Series(dtype=object)
|
|
|
250 self.contig_info[self.read_type] = pandas.read_csv(self.contig_coverage_file, header=None, index_col=None, sep='\t').sort_values(1, axis=0, ascending=False)
|
|
|
251 self.contig_info[self.read_type].columns = ['contig', 'size', 'coverage']
|
|
26
|
252 mean_coverage = (self.contig_info[self.read_type].iloc[:, 1] * self.contig_info[self.read_type].iloc[:, 2]).sum() / self.contig_info[self.read_type].iloc[:, 1].sum()
|
|
|
253 if mean_coverage <= self.ont_coverage_min:
|
|
|
254 warning = '%s mean coverage ({:.0f}X) is less than the recommended minimum ({:.0f}X).'.format(mean_coverage, self.ont_coverage_min) % self.read_type
|
|
1
|
255 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
|
256 # Report if some contigs have low coverage.
|
|
|
257 low_coverage = self.contig_info[self.read_type].loc[self.contig_info[self.read_type]['coverage'] < self.ont_coverage_min, :]
|
|
|
258 if low_coverage.shape[0] >= 0:
|
|
|
259 for contig_i in range(low_coverage.shape[0]):
|
|
|
260 warning = '%s coverage of {:s} ({:.0f}X) is less than the recommended minimum ({:.0f}X).'.format(low_coverage.iloc[contig_i, 0], low_coverage.iloc[contig_i, 2], self.ont_coverage_min) % self.read_type
|
|
|
261 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
26
|
262 # See if some contigs have anonymously low coverage.
|
|
|
263 fold_coverage = self.contig_info[self.read_type]['coverage'] / mean_coverage
|
|
1
|
264 low_coverage = self.contig_info[self.read_type].loc[fold_coverage < 1 / 5, :]
|
|
8
|
265 if low_coverage.shape[0] >= 0:
|
|
1
|
266 for contig_i in range(low_coverage.shape[0]):
|
|
26
|
267 warning = '%s coverage of {:s} ({:.0f}X) is less than 1/5 the mean coverage ({:.0f}X).'.format(low_coverage.iloc[contig_i, 0], low_coverage.iloc[contig_i, 2], mean_coverage) % self.read_type
|
|
1
|
268 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
0
|
269
|
|
|
270 def load_fasta(self, fasta):
|
|
|
271 sequence = pandas.Series(dtype=object)
|
|
|
272 for contig in SeqIO.parse(fasta, 'fasta'):
|
|
|
273 sequence[contig.id] = contig
|
|
|
274 return sequence
|
|
|
275
|
|
|
276 def load_assembly(self):
|
|
|
277 self.assembly = self.load_fasta(self.assembly_fasta_file)
|
|
|
278 self.num_assembly_contigs = len(self.assembly)
|
|
26
|
279 self.assembly_size = self.format_kmg(sum([len(x) for x in self.assembly]), decimals=1)
|
|
0
|
280
|
|
26
|
281 def info_illumina_fastq(self, illumina_read_files):
|
|
0
|
282 self.ofh.write("\nXXXXXX In info_illumina_fastq\n\n")
|
|
|
283 if self.gzipped:
|
|
|
284 opener = 'gunzip -c'
|
|
|
285 else:
|
|
|
286 opener = 'cat'
|
|
26
|
287 for fastq_file in illumina_read_files:
|
|
|
288 command = ' '.join([opener,
|
|
|
289 fastq_file,
|
|
|
290 '| awk \'{getline;s += length($1);getline;getline;}END{print s/(NR/4)"\t"(NR/4)"\t"s}\''])
|
|
|
291 output = self.run_command(command)
|
|
|
292 self.ofh.write("output:\n%s\n" % str(output))
|
|
|
293 self.ofh.write("re.split('\\t', self.run_command(command)[0]:\n%s\n" % str(re.split('\\t', self.run_command(command)[0])))
|
|
|
294 values = []
|
|
|
295 for i in re.split('\\t', self.run_command(command)[0]):
|
|
|
296 if i == '':
|
|
|
297 values.append(float('nan'))
|
|
|
298 else:
|
|
|
299 values.append(float(i))
|
|
|
300 self.ofh.write("values:\n%s\n" % str(values))
|
|
|
301 self.ofh.write("values[0]:\n%s\n" % str(values[0]))
|
|
|
302 self.illumina_length_mean += values[0]
|
|
|
303 self.ofh.write("values[1]:\n%s\n" % str(values[1]))
|
|
|
304 self.illumina_read_count += int(values[1])
|
|
|
305 self.ofh.write("values[2]:\n%s\n" % str(values[2]))
|
|
|
306 self.illumina_bases += int(values[2])
|
|
|
307 self.illumina_length_mean /= 2
|
|
0
|
308 self.illumina_bases = self.format_kmg(self.illumina_bases, decimals=1)
|
|
|
309
|
|
|
310 def start_doc(self):
|
|
13
|
311 header_text = 'Analysis of ' + self.analysis_name
|
|
|
312 self.doc = MdUtils(file_name=self.report_md, title=header_text)
|
|
0
|
313
|
|
14
|
314 def add_table_of_contents(self):
|
|
|
315 self.doc.create_marker(text_marker="TableOfContents")
|
|
|
316 self.doc.new_line()
|
|
|
317 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
318 self.doc.new_line()
|
|
|
319
|
|
0
|
320 def add_run_information(self):
|
|
|
321 self.ofh.write("\nXXXXXX In add_run_information\n\n")
|
|
|
322 self.doc.new_line()
|
|
|
323 self.doc.new_header(1, 'Run information')
|
|
|
324 # Tables in md.utils are implemented as a wrapping function.
|
|
|
325 Table_list = [
|
|
|
326 "Category",
|
|
|
327 "Information",
|
|
|
328 "Date",
|
|
|
329 date.today(),
|
|
|
330 "ONT FAST5",
|
|
14
|
331 self.wordwrap_markdown(self.ont_fast5),
|
|
0
|
332 "ONT FASTQ",
|
|
14
|
333 self.wordwrap_markdown(self.ont_raw_fastq),
|
|
0
|
334 "Illumina FASTQ",
|
|
21
|
335 self.wordwrap_markdown(self.illumina_fastq),
|
|
0
|
336 "Assembly",
|
|
|
337 self.wordwrap_markdown(self.assembly_name),
|
|
|
338 "Reference",
|
|
|
339 self.wordwrap_markdown(self.dbkey),
|
|
|
340 ]
|
|
|
341 self.doc.new_table(columns=2, rows=7, text=Table_list, text_align='left')
|
|
|
342 self.doc.new_line()
|
|
14
|
343 # FIXME: the following doesn't work.
|
|
|
344 # self.add_table_of_contents()
|
|
0
|
345 self.doc.new_line()
|
|
|
346
|
|
|
347 def add_ont_library_information(self):
|
|
|
348 self.ofh.write("\nXXXXXX In add_ont_library_information\n\n")
|
|
|
349 if self.ont_n50 is None:
|
|
|
350 return
|
|
|
351 self.doc.new_line()
|
|
|
352 self.doc.new_header(2, 'ONT library statistics')
|
|
|
353 Table_List = [
|
|
|
354 "Category",
|
|
|
355 "Quantity",
|
|
|
356 "ONT N50",
|
|
|
357 '{:,}'.format(self.ont_n50),
|
|
|
358 "ONT reads",
|
|
|
359 '{:,}'.format(self.ont_read_count),
|
|
|
360 "ONT bases",
|
|
|
361 '{:s}'.format(self.ont_bases),
|
|
|
362 "Illumina FASTQ",
|
|
17
|
363 "N/A",
|
|
0
|
364 "Assembly",
|
|
|
365 self.wordwrap_markdown(self.assembly_name),
|
|
|
366 "Reference",
|
|
|
367 self.wordwrap_markdown(self.dbkey),
|
|
|
368 ]
|
|
|
369 self.doc.new_table(columns=2, rows=7, text=Table_List, text_align='left')
|
|
|
370 self.doc.new_line()
|
|
|
371
|
|
|
372 def add_illumina_library_information(self):
|
|
|
373 self.ofh.write("\nXXXXXX In add_illumina_library_information\n\n")
|
|
14
|
374 if self.illumina_length_mean is None:
|
|
0
|
375 return
|
|
|
376 self.doc.new_line()
|
|
|
377 self.doc.new_header(2, 'Illumina library statistics')
|
|
|
378 Table_List = [
|
|
|
379 "Illumina Info.",
|
|
|
380 "Quantity",
|
|
|
381 'Illumina mean length',
|
|
|
382 '{:.1f}'.format(self.illumina_length_mean),
|
|
|
383 'Illumina reads',
|
|
|
384 '{:,}'.format(self.illumina_read_count),
|
|
|
385 'Illumina bases',
|
|
|
386 '{:s}'.format(self.illumina_bases)
|
|
|
387 ]
|
|
|
388 self.doc.new_table(columns=2, rows=4, text=Table_List, text_align='left')
|
|
|
389
|
|
8
|
390 def evaluate_assembly(self):
|
|
1
|
391 assembly_info = pandas.read_csv(self.compute_sequence_length_file, sep='\t', header=None)
|
|
|
392 assembly_info.columns = ['contig', 'length']
|
|
|
393 self.contig_sizes = assembly_info
|
|
|
394 # Take a look at the number of contigs, their sizes,
|
|
|
395 # and circularity. Warn if things don't look good.
|
|
|
396 if assembly_info.shape[0] > 4:
|
|
|
397 warning = 'Assembly produced {:d} contigs, more than ususally expected; assembly may be fragmented'.format(assembly_info.shape[0])
|
|
|
398 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
|
399 small_contigs = assembly_info.loc[assembly_info['length'] <= 3000, :]
|
|
|
400 if small_contigs.shape[0] > 0:
|
|
|
401 warning = 'Assembly produced {:d} small contigs ({:s}); assembly may include spurious sequences.'.format(small_contigs.shape[0], ', '.join(small_contigs['contig']))
|
|
|
402 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
|
403
|
|
0
|
404 def add_assembly_information(self):
|
|
|
405 self.ofh.write("\nXXXXXX In add_assembly_information\n\n")
|
|
|
406 if self.assembly_fasta_file is None:
|
|
|
407 return
|
|
|
408 self.load_assembly()
|
|
|
409 self.doc.new_line()
|
|
|
410 self.doc.new_header(2, 'Assembly statistics')
|
|
|
411 Table_List = [
|
|
|
412 "Category",
|
|
|
413 "Information",
|
|
|
414 "Contigs",
|
|
|
415 str(self.num_assembly_contigs),
|
|
|
416 "Assembly size",
|
|
|
417 str(self.assembly_size),
|
|
|
418 ]
|
|
|
419 self.doc.new_table(columns=2, rows=3, text=Table_List, text_align='left')
|
|
|
420
|
|
|
421 def info_ont_fastq(self, fastq_file):
|
|
|
422 self.ofh.write("\nXXXXXX In info_ont_fastq, fastq_file:\n%s\n\n" % str(fastq_file))
|
|
|
423 opener = 'cat'
|
|
|
424 if self.gzipped:
|
|
|
425 opener = 'gunzip -c'
|
|
|
426 else:
|
|
|
427 opener = 'cat'
|
|
|
428 command = ' '.join([opener,
|
|
|
429 fastq_file,
|
|
|
430 '| awk \'{getline;print length($0);s += length($1);getline;getline;}END{print "+"s}\'',
|
|
|
431 '| sort -gr',
|
|
|
432 '| awk \'BEGIN{bp = 0;f = 0}',
|
|
|
433 '{if(NR == 1){sub(/+/, "", $1);s=$1}else{bp += $1;if(bp > s / 2 && f == 0){n50 = $1;f = 1}}}',
|
|
|
434 'END{printf "%d\\t%d\\t%d\\n", n50, (NR - 1), s;exit}\''])
|
|
|
435 result = list(re.split('\\t', self.run_command(command)[0]))
|
|
|
436 if result[1] == '0':
|
|
21
|
437 warning = 'No ONT reads found'
|
|
|
438 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
14
|
439 self.ont_n50, self.ont_read_count, ont_raw_bases = [int(i) for i in result]
|
|
0
|
440 command = ' '.join([opener,
|
|
|
441 fastq_file,
|
|
|
442 '| awk \'{getline;print length($0);getline;getline;}\''])
|
|
|
443 result = self.run_command(command)
|
|
|
444 result = list(filter(lambda x: x != '', result))
|
|
14
|
445 self.ont_bases = self.format_kmg(ont_raw_bases, decimals=1)
|
|
|
446 if self.ont_n50 <= self.ont_n50_min:
|
|
|
447 warning = 'ONT N50 (%s) is less than the recommended minimum (%s)' % (str(self.ont_n50), str(self.ont_n50_min))
|
|
1
|
448 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
0
|
449
|
|
|
450 def wordwrap_markdown(self, string):
|
|
|
451 if string:
|
|
|
452 if len(string) < 35:
|
|
|
453 return string
|
|
|
454 else:
|
|
|
455 if '/' in string:
|
|
|
456 adjust = string.split('/')
|
|
|
457 out = ''
|
|
|
458 max = 35
|
|
|
459 for i in adjust:
|
|
|
460 out = out + '/' + i
|
|
|
461 if len(out) > max:
|
|
|
462 out += '<br>'
|
|
|
463 max += 35
|
|
|
464 return out
|
|
|
465 else:
|
|
|
466 out = [string[i:i + 35] for i in range(0, len(string), 50)]
|
|
|
467 return '<br>'.join(out)
|
|
|
468 else:
|
|
|
469 return string
|
|
|
470
|
|
|
471 def add_contig_info(self):
|
|
|
472 self.ofh.write("\nXXXXXX In add_contig_info\n\n")
|
|
26
|
473 if self.contig_info is None or self.read_type not in self.contig_info.index:
|
|
0
|
474 return
|
|
26
|
475 self.doc.new_line()
|
|
|
476 self.doc.new_header(2, 'Assembly coverage by ' + self.read_type)
|
|
|
477 Table_List = ["Contig", "Length (bp)", "Coverage (X)"]
|
|
|
478 formatted = self.contig_info[self.read_type].copy()
|
|
|
479 formatted.iloc[:, 1] = formatted.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
|
|
|
480 for i in range(self.contig_info[self.read_type].shape[0]):
|
|
|
481 Table_List = Table_List + formatted.iloc[i, :].values.tolist()
|
|
|
482 row_count = int(len(Table_List) / 3)
|
|
|
483 self.doc.new_table(columns=3, rows=row_count, text=Table_List, text_align='left')
|
|
0
|
484
|
|
|
485 def add_assembly_notes(self):
|
|
|
486 self.ofh.write("\nXXXXXX In add_assembly_notes\n\n")
|
|
|
487 if len(self.assembly_notes) == 0:
|
|
|
488 return
|
|
|
489 self.doc.new_line()
|
|
|
490 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
491 self.doc.new_line()
|
|
|
492 self.doc.new_header(2, self.assembly_notes_title)
|
|
1
|
493 for note in self.assembly_notes:
|
|
|
494 self.doc.new_line(note)
|
|
0
|
495
|
|
|
496 def add_contamination(self):
|
|
|
497 self.ofh.write("\nXXXXXX In add_contamination\n\n")
|
|
2
|
498 if self.kraken2_report_file is None:
|
|
0
|
499 return
|
|
2
|
500 # Read in the Kraken fractions and pull out the useful parts
|
|
8
|
501 kraken_fracs = pandas.read_csv(self.kraken2_report_file, delimiter='\t', header=None)
|
|
|
502 kraken_fracs.index = kraken_fracs.iloc[:, 4].values
|
|
|
503 kraken_fracs = kraken_fracs.loc[kraken_fracs.iloc[:, 3].str.match('[UG]1?'), :]
|
|
|
504 kraken_fracs = kraken_fracs.loc[(kraken_fracs.iloc[:, 0] >= 1) | (kraken_fracs.iloc[:, 3] == 'U'), :]
|
|
|
505 kraken_fracs = kraken_fracs.iloc[:, [0, 1, 3, 5]]
|
|
|
506 kraken_fracs.columns = ['Fraction', 'Reads', 'Level', 'Taxa']
|
|
|
507 kraken_fracs['Fraction'] = (kraken_fracs['Fraction'] / 100).round(4)
|
|
|
508 kraken_fracs.sort_values(by='Fraction', inplace=True, ascending=False)
|
|
|
509 kraken_fracs['Taxa'] = kraken_fracs['Taxa'].str.lstrip()
|
|
0
|
510 self.doc.new_line()
|
|
|
511 self.doc.new_header(2, 'Contamination check')
|
|
10
|
512 self.doc.new_line(self.read_type + ' classifications')
|
|
|
513 self.doc.new_line()
|
|
|
514 Table_List = ["Percent of Reads", "Reads", "Level", "Label"]
|
|
|
515 for index, row in kraken_fracs.iterrows():
|
|
|
516 Table_List = Table_List + row.tolist()
|
|
|
517 row_count = int(len(Table_List) / 4)
|
|
|
518 self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left')
|
|
|
519 if self.contamination_methods_title not in self.methods:
|
|
|
520 self.methods[self.contamination_methods_title] = ''
|
|
11
|
521 method = '%s was used to assign the raw reads into taxa.' % self.kraken2_version.rstrip('report')
|
|
0
|
522 self.methods[self.contamination_methods_title] = self.methods[self.contamination_methods_title].append(pandas.Series(method))
|
|
|
523
|
|
|
524 def add_alignment(self):
|
|
|
525 self.ofh.write("\nXXXXXX In add_alignment\n\n")
|
|
13
|
526 if self.quast_report_file is not None:
|
|
|
527 # Process quast values.
|
|
|
528 quast_report = pandas.read_csv(self.quast_report_file, header=0, index_col=0, sep='\t')
|
|
|
529 quast_mismatches = int(float(quast_report.loc['# mismatches per 100 kbp', :][0]) * (float(quast_report.loc['Total length (>= 0 bp)', :][0]) / 100000.))
|
|
|
530 quast_indels = int(float(quast_report.loc['# indels per 100 kbp', :][0]) * (float(quast_report.loc['Total length (>= 0 bp)', :][0]) / 100000.))
|
|
|
531 self.doc.new_line()
|
|
|
532 self.doc.new_header(level=2, title=self.alignment_title)
|
|
|
533 self.doc.new_line()
|
|
|
534 self.doc.new_header(level=3, title=self.snp_indel_title)
|
|
|
535 Table_1 = [
|
|
|
536 "Category",
|
|
|
537 "Quantity",
|
|
|
538 'SNPs',
|
|
|
539 '{:,}'.format(quast_mismatches),
|
|
|
540 'Small indels',
|
|
|
541 '{:,}'.format(quast_indels)
|
|
|
542 ]
|
|
0
|
543 self.doc.new_table(columns=2, rows=3, text=Table_1, text_align='left')
|
|
|
544 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
545 self.doc.new_line()
|
|
13
|
546 # TODO: self.alignment_notes is not currently populated.
|
|
0
|
547 if len(self.alignment_notes) > 0:
|
|
|
548 self.doc.new_header(level=3, title=self.alignment_notes_title)
|
|
|
549 for note in self.alignment_notes:
|
|
|
550 self.doc.new_line(note)
|
|
13
|
551 if len(self.circos_files) > 0:
|
|
|
552 # Add circos PNG files.
|
|
|
553 for circos_file in self.circos_files:
|
|
|
554 contig = os.path.basename(circos_file)
|
|
|
555 contig_title = 'Alignment to %s' % contig
|
|
|
556 self.doc.new_line()
|
|
|
557 self.doc.new_header(level=3, title=contig_title)
|
|
25
|
558 self.doc.new_line('Blue color indicates query sequences aligned to the reference sequence, which is shown in yellow')
|
|
13
|
559 self.doc.new_line(self.doc.new_inline_image(text='contig_title', path=os.path.abspath(circos_file)))
|
|
|
560 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
561 self.doc.new_line()
|
|
21
|
562 if self.dbkey == 'ref_genome':
|
|
|
563 headers = ["* Chromosome - NC_007530.2 Bacillus anthracis str. 'Ames Ancestor', complete sequence",
|
|
|
564 "* pXO1 - NC_007322.2 Bacillus anthracis str. 'Ames Ancestor' plasmid pXO1, complete sequence",
|
|
|
565 "* pXO2 - NC_007323.3 Bacillus anthracis str. 'Ames Ancestor' plasmid pXO2, complete sequence"]
|
|
|
566 method = '\n'.join(headers)
|
|
|
567 self.methods[self.reference_genome_title] = self.methods[self.reference_genome_title].append(pandas.Series(method))
|
|
|
568 method = 'The genome assembly was aligned against the reference sequence using %s.' % self.dnadiff_version
|
|
0
|
569 self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method))
|
|
|
570
|
|
|
571 def add_features(self):
|
|
|
572 self.ofh.write("\nXXXXXX In add_features\n\n")
|
|
|
573 if len(self.feature_bed_files) == 0:
|
|
|
574 return
|
|
|
575 for bbf in self.feature_bed_files:
|
|
|
576 if os.path.getsize(bbf) > 0:
|
|
|
577 best = pandas.read_csv(filepath_or_buffer=bbf, sep='\t', header=None)
|
|
|
578 self.feature_hits[os.path.basename(bbf)] = best
|
|
|
579 if len(self.feature_hits) == 0:
|
|
|
580 return
|
|
|
581 self.ofh.write("self.feature_hits: %s\n" % str(self.feature_hits))
|
|
|
582 self.doc.new_line()
|
|
|
583 self.doc.new_header(level=2, title=self.feature_title)
|
|
|
584 for feature_name in self.feature_hits.index.tolist():
|
|
|
585 self.ofh.write("feature_name: %s\n" % str(feature_name))
|
|
|
586 features = self.feature_hits[feature_name].copy()
|
|
|
587 self.ofh.write("features: %s\n" % str(features))
|
|
|
588 if features.shape[0] == 0:
|
|
|
589 continue
|
|
|
590 features.iloc[:, 1] = features.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
|
|
|
591 features.iloc[:, 2] = features.iloc[:, 2].apply(lambda x: '{:,}'.format(x))
|
|
|
592 self.doc.new_line()
|
|
|
593 self.doc.new_header(level=3, title=feature_name)
|
|
|
594 if (features.shape[0] == 0):
|
|
|
595 continue
|
|
|
596 for contig in pandas.unique(features.iloc[:, 0]):
|
|
|
597 self.ofh.write("contig: %s\n" % str(contig))
|
|
|
598 self.doc.new_line(contig)
|
|
|
599 contig_features = features.loc[(features.iloc[:, 0] == contig), :]
|
|
|
600 self.ofh.write("contig_features: %s\n" % str(contig_features))
|
|
|
601 Table_List = ['Start', 'Stop', 'Feature', 'Identity (%)', 'Strand']
|
|
|
602 for i in range(contig_features.shape[0]):
|
|
|
603 self.ofh.write("i: %s\n" % str(i))
|
|
|
604 feature = contig_features.iloc[i, :].copy(deep=True)
|
|
|
605 self.ofh.write("feature: %s\n" % str(feature))
|
|
|
606 feature[4] = '{:.3f}'.format(feature[4])
|
|
2
|
607 self.ofh.write("feature[1:].values.tolist(): %s\n" % str(feature[1:].values.tolist()))
|
|
|
608 Table_List = Table_List + feature[1:].values.tolist()
|
|
0
|
609 self.ofh.write("Table_List: %s\n" % str(Table_List))
|
|
|
610 row_count = int(len(Table_List) / 5)
|
|
|
611 self.ofh.write("row_count: %s\n" % str(row_count))
|
|
|
612 self.doc.new_line()
|
|
1
|
613 self.ofh.write("Before new_table, len(Table_List):: %s\n" % str(len(Table_List)))
|
|
|
614 self.doc.new_table(columns=5, rows=row_count, text=Table_List, text_align='left')
|
|
12
|
615 blastn_version = 'The genome assembly was queried for features using %s.' % self.blastn_version
|
|
|
616 bedtools_version = 'Feature hits were clustered using %s and the highest scoring hit for each cluster was reported.' % self.bedtools_version
|
|
0
|
617 method = '%s %s' % (blastn_version, bedtools_version)
|
|
|
618 self.methods[self.feature_methods_title] = self.methods[self.feature_methods_title].append(pandas.Series(method))
|
|
|
619
|
|
|
620 def add_feature_plots(self):
|
|
|
621 self.ofh.write("\nXXXXXX In add_feature_plots\n\n")
|
|
|
622 if len(self.feature_png_files) == 0:
|
|
|
623 return
|
|
|
624 self.doc.new_line()
|
|
|
625 self.doc.new_header(level=2, title='Feature Plots')
|
|
|
626 self.doc.new_paragraph('Only contigs with features are shown')
|
|
|
627 for feature_png_file in self.feature_png_files:
|
|
|
628 self.doc.new_line(self.doc.new_inline_image(text='Analysis', path=os.path.abspath(feature_png_file)))
|
|
|
629
|
|
|
630 def add_mutations(self):
|
|
|
631 self.ofh.write("\nXXXXXX In add_mutations\n\n")
|
|
|
632 if len(self.mutation_regions_tsv_files) == 0:
|
|
|
633 return
|
|
8
|
634 try:
|
|
0
|
635 mutation_regions = pandas.read_csv(self.mutation_regions_bed_file, sep='\t', header=0, index_col=False)
|
|
|
636 except Exception:
|
|
|
637 # Likely an empty file.
|
|
|
638 return
|
|
|
639 amr_mutations = pandas.Series(dtype=object)
|
|
|
640 for region_i in range(mutation_regions.shape[0]):
|
|
|
641 region = mutation_regions.iloc[region_i, :]
|
|
|
642 region_name = str(region['name'])
|
|
|
643 self.ofh.write("Processing mutations for region %s\n" % region_name)
|
|
|
644 region_mutations_tsv_name = '%s_mutations.tsv' % region_name
|
|
|
645 if region_mutations_tsv_name not in self.mutation_regions_tsv_files:
|
|
|
646 continue
|
|
|
647 region_mutations_tsv = self.mutation_regions_tsv_files[region_mutations_tsv_name]
|
|
8
|
648 try:
|
|
0
|
649 region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False)
|
|
|
650 except Exception:
|
|
|
651 region_mutations = pandas.DataFrame()
|
|
|
652 if region_mutations.shape[0] == 0:
|
|
|
653 continue
|
|
1
|
654 # Figure out what kind of mutations are in this region.
|
|
0
|
655 region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index)
|
|
|
656 region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel'
|
|
|
657 region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index)
|
|
|
658 region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index)
|
|
|
659 region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1)
|
|
|
660 region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']]
|
|
|
661 amr_mutations[region['name']] = region_mutations
|
|
2
|
662 if (amr_mutations.shape[0] > 0):
|
|
|
663 # Report the mutations.
|
|
0
|
664 self.doc.new_line()
|
|
2
|
665 self.doc.new_header(level=2, title=self.mutation_title)
|
|
|
666 for region_name in amr_mutations.index.tolist():
|
|
|
667 region_mutations = amr_mutations[region_name].copy()
|
|
|
668 self.doc.new_line()
|
|
|
669 self.doc.new_header(level=3, title=region_name)
|
|
|
670 if (region_mutations.shape[0] == 0):
|
|
|
671 self.doc.append('None')
|
|
|
672 continue
|
|
|
673 region_mutations.iloc[:, 1] = region_mutations.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
|
|
|
674 Table_List = ['Reference contig', 'Position', 'Reference', 'Alternate', 'Drug', 'Note']
|
|
|
675 for i in range(region_mutations.shape[0]):
|
|
|
676 Table_List = Table_List + region_mutations.iloc[i, [0, 1, 3, 4, 5, 6]].values.tolist()
|
|
|
677 row_count = int(len(Table_List) / 6)
|
|
|
678 self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left')
|
|
18
|
679 if os.path.getsize(self.errors_file) > 0:
|
|
|
680 # Report the errors encountered when attempting
|
|
|
681 # to find mutations in the sample.
|
|
|
682 self.doc.new_line()
|
|
|
683 self.doc.new_header(level=2, title=self.mutation_errors_title)
|
|
|
684 with open(self.errors_file, 'r') as efh:
|
|
|
685 for i, line in enumerate(efh):
|
|
|
686 line = line.strip()
|
|
|
687 if line:
|
|
|
688 self.doc.new_line('* %s' % line)
|
|
12
|
689 method = '%s reads were mapped to the reference sequence using %s.' % (self.read_type, self.minimap2_version)
|
|
0
|
690 self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
|
|
13
|
691 method = 'Mutations were identified using %s and %s.' % (self.samtools_version, self.varscan_version)
|
|
0
|
692 self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
|
|
|
693
|
|
|
694 def add_amr_matrix(self):
|
|
|
695 self.ofh.write("\nXXXXXX In add_amr_matrix\n\n")
|
|
|
696 # Make sure that we have an AMR matrix to plot
|
|
1
|
697 if len(self.amr_matrix_files) == 0:
|
|
|
698 return
|
|
|
699 self.doc.new_line()
|
|
|
700 self.doc.new_header(level=2, title=self.amr_matrix_title)
|
|
|
701 self.doc.new_line('AMR genes and mutations with their corresponding drugs')
|
|
|
702 for amr_matrix_file in self.amr_matrix_files:
|
|
|
703 self.doc.new_line(self.doc.new_inline_image(text='AMR genes and mutations with their corresponding drugs',
|
|
|
704 path=os.path.abspath(amr_matrix_file)))
|
|
0
|
705
|
|
|
706 def add_large_indels(self):
|
|
|
707 self.ofh.write("\nXXXXXX In add_large_indels\n\n")
|
|
1
|
708 large_indels = pandas.Series(dtype='float64')
|
|
|
709 # Pull in insertions.
|
|
|
710 try:
|
|
|
711 reference_insertions = pandas.read_csv(filepath_or_buffer=self.reference_insertions_file, sep='\t', header=None)
|
|
|
712 except Exception:
|
|
|
713 reference_insertions = pandas.DataFrame()
|
|
|
714 try:
|
|
|
715 genome_insertions = pandas.read_csv(filepath_or_buffer=self.genome_insertions_file, sep='\t', header=None)
|
|
|
716 except Exception:
|
|
|
717 genome_insertions = pandas.DataFrame()
|
|
|
718 large_indels['Reference insertions'] = reference_insertions
|
|
|
719 large_indels['Query insertions'] = genome_insertions
|
|
|
720 # Pull in deletions.
|
|
|
721 try:
|
|
|
722 amr_deletions = pandas.read_csv(filepath_or_buffer=self.amr_deletion_file, sep='\t', header=None)
|
|
|
723 except Exception:
|
|
|
724 amr_deletions = pandas.DataFrame()
|
|
|
725 if amr_deletions.shape[0] > 0:
|
|
|
726 amr_deletions.columns = ['contig', 'start', 'stop', 'name', 'type', 'drug', 'note']
|
|
|
727 amr_deletions = amr_deletions.loc[amr_deletions['type'].isin(['large-deletion', 'any']), :]
|
|
|
728 self.doc.new_line()
|
|
|
729 self.doc.new_header(level=2, title=self.large_indel_title)
|
|
25
|
730 self.doc.new_line('This section is informative only when your isolates were identified as *Bacillus anthracis* strains')
|
|
1
|
731 for genome in ['Reference insertions', 'Query insertions']:
|
|
|
732 genome_indels = large_indels[genome].copy()
|
|
|
733 self.doc.new_line()
|
|
|
734 self.doc.new_header(level=3, title=genome)
|
|
|
735 if (genome_indels.shape[0] == 0):
|
|
|
736 continue
|
|
|
737 genome_indels.iloc[:, 1] = genome_indels.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
|
|
|
738 genome_indels.iloc[:, 2] = genome_indels.iloc[:, 2].apply(lambda x: '{:,}'.format(x))
|
|
|
739 genome_indels.iloc[:, 3] = genome_indels.iloc[:, 3].apply(lambda x: '{:,}'.format(x))
|
|
|
740 Table_List = [
|
|
|
741 'Reference contig', 'Start', 'Stop', 'Size (bp)'
|
|
|
742 ]
|
|
|
743 for i in range(genome_indels.shape[0]):
|
|
|
744 Table_List = Table_List + genome_indels.iloc[i, :].values.tolist()
|
|
|
745 row_count = int(len(Table_List) / 4)
|
|
|
746 self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left')
|
|
12
|
747 method = 'Large insertions or deletions were found as the complement of aligned regions using %s.' % self.bedtools_version
|
|
1
|
748 self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method))
|
|
|
749 self.doc.new_line()
|
|
|
750 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
751 self.doc.new_line()
|
|
0
|
752
|
|
28
|
753 def add_lrn_risk_info(self):
|
|
|
754 self.ofh.write("\nXXXXXX In add_lrn_risk_info\n\n")
|
|
|
755 if self.lrn_risk_amr_file is None and self.lrn_risk_blacklist_file is None and self.lrn_risk_vf_file is None:
|
|
|
756 return
|
|
30
|
757 if not os.path.isfile(self.lrn_risk_amr_file) and not os.path.isfile(self.lrn_risk_blacklist_file) and not os.path.isfile(self.lrn_risk_vf_file):
|
|
|
758 return
|
|
|
759 if os.path.getsize(self.lrn_risk_amr_file) == 0 and os.path.getsize(self.lrn_risk_blacklist_file) == 0 and os.path.getsize(self.lrn_risk_vf_file) == 0:
|
|
|
760 return
|
|
28
|
761 self.doc.new_line()
|
|
|
762 self.doc.new_header(level=2, title=self.lrn_risk_title)
|
|
|
763 # Process self.lrn_risk_amr_file.
|
|
|
764 try:
|
|
|
765 lrn_risk_amr = pandas.read_csv(filepath_or_buffer=self.lrn_risk_amr_file, sep='\t', header=0)
|
|
|
766 except Exception:
|
|
|
767 lrn_risk_amr = pandas.DataFrame()
|
|
|
768 if lrn_risk_amr.shape[0] > 0:
|
|
|
769 self.doc.new_line()
|
|
|
770 self.doc.new_header(level=2, title="AMR Determinant Distribution")
|
|
|
771 self.doc.new_line()
|
|
|
772 Table_List = ["Gene", "Contig", "% Identity", "% Coverage", "E-Value", "Annotation", "Comparison to Publicly Available Genomes"]
|
|
|
773 for index, row in lrn_risk_amr.iterrows():
|
|
|
774 Table_List = Table_List + row.tolist()
|
|
|
775 row_count = int(len(Table_List) / 7)
|
|
|
776 self.doc.new_table(columns=7, rows=row_count, text=Table_List, text_align='left')
|
|
|
777 # Process self.lrn_risk_blacklist_file.
|
|
|
778 try:
|
|
|
779 lrn_risk_blacklist = pandas.read_csv(filepath_or_buffer=self.lrn_risk_blacklist_file, sep='\t', header=0)
|
|
|
780 except Exception:
|
|
|
781 lrn_risk_blacklist = pandas.DataFrame()
|
|
|
782 if lrn_risk_blacklist.shape[0] > 0:
|
|
|
783 self.doc.new_line()
|
|
|
784 self.doc.new_header(level=2, title="Blacklisted High-risk Virulence Factors")
|
|
|
785 self.doc.new_line()
|
|
|
786 Table_List = ["Blacklisted Gene", "Reason", "Risk Category"]
|
|
|
787 for index, row in lrn_risk_blacklist.iterrows():
|
|
|
788 Table_List = Table_List + row.tolist()
|
|
|
789 row_count = int(len(Table_List) / 3)
|
|
|
790 self.doc.new_table(columns=3, rows=row_count, text=Table_List, text_align='left')
|
|
|
791 # Process self.lrn_risk_vf_file.
|
|
|
792 try:
|
|
|
793 lrn_risk_vf = pandas.read_csv(filepath_or_buffer=self.lrn_risk_vf_file, sep='\t', header=0)
|
|
|
794 except Exception:
|
|
|
795 lrn_risk_vf = pandas.DataFrame()
|
|
|
796 if lrn_risk_vf.shape[0] > 0:
|
|
|
797 self.doc.new_line()
|
|
|
798 self.doc.new_header(level=2, title="Virulence Factor Distribution")
|
|
|
799 self.doc.new_line()
|
|
|
800 Table_List = ["Gene", "Contig", "% Identity", "% Coverage", "E-Value", "Annotation", "Comparison to Publicly Available Genomes"]
|
|
|
801 for index, row in lrn_risk_vf.iterrows():
|
|
|
802 Table_List = Table_List + row.tolist()
|
|
|
803 row_count = int(len(Table_List) / 7)
|
|
|
804 self.doc.new_table(columns=7, rows=row_count, text=Table_List, text_align='left')
|
|
|
805 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
806 self.doc.new_line()
|
|
|
807
|
|
0
|
808 def add_plasmids(self):
|
|
8
|
809 try:
|
|
1
|
810 plasmids = pandas.read_csv(filepath_or_buffer=self.plasmids_file, sep='\t', header=0)
|
|
|
811 except Exception:
|
|
0
|
812 return
|
|
|
813 plasmids = plasmids.copy()
|
|
|
814 self.doc.new_line()
|
|
1
|
815 self.doc.new_header(level=2, title=self.plasmid_title)
|
|
0
|
816 if (plasmids.shape[0] == 0):
|
|
|
817 self.doc.new_line('None')
|
|
|
818 return
|
|
|
819 plasmids.iloc[:, 3] = plasmids.iloc[:, 3].apply(lambda x: '{:,}'.format(x))
|
|
|
820 plasmids.iloc[:, 4] = plasmids.iloc[:, 4].apply(lambda x: '{:,}'.format(x))
|
|
|
821 plasmids.iloc[:, 5] = plasmids.iloc[:, 5].apply(lambda x: '{:,}'.format(x))
|
|
1
|
822 Table_List = ['Genome contig', 'Plasmid hit', 'Plasmid acc.', 'Contig size', 'Aliged', 'Plasmid size']
|
|
0
|
823 for i in range(plasmids.shape[0]):
|
|
|
824 Table_List = Table_List + plasmids.iloc[i, 0:6].values.tolist()
|
|
|
825 row_count = int(len(Table_List) / 6)
|
|
|
826 self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left')
|
|
12
|
827 method = 'The plasmid reference database was queried against the genome assembly using %s.' % self.minimap2_version
|
|
0
|
828 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
|
|
2
|
829 method = 'The resulting BAM was converted to a PSL using a custom version of sam2psl.'
|
|
0
|
830 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
|
|
|
831 method = 'Plasmid-to-genome hits were resolved using the pChunks algorithm.'
|
|
|
832 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
|
|
|
833
|
|
|
834 def add_methods(self):
|
|
|
835 self.ofh.write("\nXXXXXX In add_methods\n\n")
|
|
|
836 if len(self.methods) == 0:
|
|
|
837 return
|
|
|
838 self.doc.new_line()
|
|
|
839 self.doc.new_header(level=2, title=self.methods_title)
|
|
|
840 for methods_section in self.methods.index.tolist():
|
|
|
841 if self.methods[methods_section] is None or len(self.methods[methods_section]) == 0:
|
|
|
842 continue
|
|
|
843 self.doc.new_line()
|
|
|
844 self.doc.new_header(level=3, title=methods_section)
|
|
|
845 self.doc.new_paragraph(' '.join(self.methods[methods_section]))
|
|
24
|
846 self.doc.new_line('<div style="page-break-after: always;"></div>')
|
|
|
847 self.doc.new_line()
|
|
0
|
848
|
|
|
849 def add_summary(self):
|
|
|
850 self.ofh.write("\nXXXXXX In add_summary\n\n")
|
|
|
851 # Add summary title
|
|
|
852 self.doc.new_header(level=1, title=self.summary_title)
|
|
|
853 # First section of Summary
|
|
|
854 self.doc.new_header(level=1, title='CDC Advisory')
|
|
|
855 self.doc.new_paragraph(CDC_ADVISORY)
|
|
|
856 self.doc.new_line()
|
|
|
857 self.add_run_information()
|
|
|
858 self.add_ont_library_information()
|
|
|
859 methods = []
|
|
|
860 if self.did_guppy_ont_fast5:
|
|
|
861 methods += ['ONT reads were basecalled using guppy']
|
|
|
862 if self.did_qcat_ont_fastq:
|
|
|
863 methods += ['ONT reads were demultiplexed and trimmed using qcat']
|
|
|
864 self.methods[self.basecalling_methods_title] = pandas.Series(methods)
|
|
26
|
865 self.add_illumina_library_information()
|
|
|
866 self.add_assembly_information()
|
|
1
|
867 self.add_contig_info()
|
|
|
868 self.evaluate_assembly()
|
|
21
|
869 if self.assembler_version is not None:
|
|
|
870 if self.read_type == 'ONT':
|
|
|
871 method = 'ONT reads were assembled using %s' % self.assembler_version
|
|
|
872 self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.Series(method))
|
|
|
873 # Pull in the assembly summary and look at the coverage.
|
|
|
874 assembly_info = pandas.read_csv(self.flye_assembly_info_file, header=0, index_col=0, sep='\t')
|
|
|
875 # Look for non-circular contigs.
|
|
|
876 open_contigs = assembly_info.loc[assembly_info['circ.'] == 'N', :]
|
|
|
877 if open_contigs.shape[0] > 0:
|
|
|
878 open_contig_ids = open_contigs.index.values
|
|
|
879 warning = 'Flye reported {:d} open contigs ({:s}); assembly may be incomplete.'.format(open_contigs.shape[0], ', '.join(open_contig_ids))
|
|
|
880 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
|
|
|
881 else:
|
|
|
882 method = 'Illumina reads were assembled using %s' % self.assembler_version
|
|
26
|
883 method = 'The genome assembly was polished using ONT reads and medaka.'
|
|
27
|
884 self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.Series(method))
|
|
1
|
885 self.add_assembly_notes()
|
|
0
|
886
|
|
|
887 def make_tex(self):
|
|
|
888 self.doc.new_table_of_contents(table_title='detailed run information', depth=2, marker="tableofcontents")
|
|
|
889 text = self.doc.file_data_text
|
|
|
890 text = text.replace("##--[", "")
|
|
|
891 text = text.replace("]--##", "")
|
|
|
892 self.doc.file_data_text = text
|
|
|
893 self.doc.create_md_file()
|
|
|
894
|
|
|
895 def make_report(self):
|
|
|
896 self.ofh.write("\nXXXXXX In make_report\n\n")
|
|
|
897 self.start_doc()
|
|
|
898 self.add_summary()
|
|
|
899 self.add_contamination()
|
|
|
900 self.add_alignment()
|
|
|
901 self.add_features()
|
|
|
902 self.add_feature_plots()
|
|
|
903 self.add_mutations()
|
|
|
904 self.add_large_indels()
|
|
|
905 self.add_plasmids()
|
|
|
906 self.add_amr_matrix()
|
|
28
|
907 self.add_lrn_risk_info()
|
|
0
|
908 # self.add_snps()
|
|
|
909 self.add_methods()
|
|
|
910 self.make_tex()
|
|
|
911 # It took me quite a long time to find out that the value of the -t
|
|
|
912 # (implied) argument in the following command must be 'html' instead of
|
|
|
913 # the more logical 'pdf'. see the answer from snsn in this thread:
|
|
|
914 # https://github.com/jessicategner/pypandoc/issues/186
|
|
|
915 self.ofh.write("\nXXXXX In make_report, calling pypandoc.convert_file...\n\n")
|
|
|
916 pypandoc.convert_file(self.report_md,
|
|
|
917 'html',
|
|
|
918 extra_args=['--pdf-engine=weasyprint', '-V', '-css=%s' % self.pima_css],
|
|
|
919 outputfile='pima_report.pdf')
|
|
|
920 self.ofh.close()
|
|
|
921
|
|
|
922
|
|
|
923 parser = argparse.ArgumentParser()
|
|
|
924
|
|
1
|
925 parser.add_argument('--amr_deletions_file', action='store', dest='amr_deletions_file', help='AMR deletions BED file')
|
|
|
926 parser.add_argument('--amr_matrix_png_dir', action='store', dest='amr_matrix_png_dir', help='Directory of AMR matrix PNG files')
|
|
0
|
927 parser.add_argument('--analysis_name', action='store', dest='analysis_name', help='Sample identifier')
|
|
21
|
928 parser.add_argument('--assembler_version', action='store', dest='assembler_version', default=None, help='Assembler version string')
|
|
0
|
929 parser.add_argument('--assembly_fasta_file', action='store', dest='assembly_fasta_file', help='Assembly fasta file')
|
|
|
930 parser.add_argument('--assembly_name', action='store', dest='assembly_name', help='Assembly identifier')
|
|
12
|
931 parser.add_argument('--bedtools_version', action='store', dest='bedtools_version', default=None, help='Bedtools version string')
|
|
2
|
932 parser.add_argument('--blastn_version', action='store', dest='blastn_version', default=None, help='Blastn version string')
|
|
13
|
933 parser.add_argument('--circos_png_dir', action='store', dest='circos_png_dir', help='Directory of circos PNG files')
|
|
1
|
934 parser.add_argument('--compute_sequence_length_file', action='store', dest='compute_sequence_length_file', help='Comnpute sequence length tabular file')
|
|
|
935 parser.add_argument('--contig_coverage_file', action='store', dest='contig_coverage_file', help='Contig coverage TSV file')
|
|
|
936 parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference genome identifier')
|
|
|
937 parser.add_argument('--dnadiff_snps_file', action='store', dest='dnadiff_snps_file', help='DNAdiff snps tabular file')
|
|
12
|
938 parser.add_argument('--dnadiff_version', action='store', dest='dnadiff_version', default=None, help='DNAdiff version string')
|
|
18
|
939 parser.add_argument('--errors_file', action='store', dest='errors_file', default=None, help='AMR mutations errors encountered txt file')
|
|
0
|
940 parser.add_argument('--feature_bed_dir', action='store', dest='feature_bed_dir', help='Directory of best feature hits bed files')
|
|
|
941 parser.add_argument('--feature_png_dir', action='store', dest='feature_png_dir', help='Directory of best feature hits png files')
|
|
1
|
942 parser.add_argument('--flye_assembly_info_file', action='store', dest='flye_assembly_info_file', default=None, help='Flye assembly info tabular file')
|
|
|
943 parser.add_argument('--genome_insertions_file', action='store', dest='genome_insertions_file', help='Genome insertions BED file')
|
|
26
|
944 parser.add_argument('--gzipped', action='store_true', dest='gzipped', default=False, help='Sample(s) is/are gzipped')
|
|
|
945 parser.add_argument('--illumina_forward_read_file', action='store', dest='illumina_forward_read_file', help='Illumina forward read file')
|
|
|
946 parser.add_argument('--illumina_reverse_read_file', action='store', dest='illumina_reverse_read_file', help='Illumina reverse read file')
|
|
2
|
947 parser.add_argument('--kraken2_report_file', action='store', dest='kraken2_report_file', default=None, help='kraken2 report file')
|
|
|
948 parser.add_argument('--kraken2_version', action='store', dest='kraken2_version', default=None, help='kraken2 version string')
|
|
28
|
949 parser.add_argument('--lrn_risk_amr_file', action='store', dest='lrn_risk_amr_file', default=None, help='LRN RISK AMR TSV file')
|
|
|
950 parser.add_argument('--lrn_risk_blacklist_file', action='store', dest='lrn_risk_blacklist_file', default=None, help='LRN RISK blacklist TSV file')
|
|
|
951 parser.add_argument('--lrn_risk_vf_file', action='store', dest='lrn_risk_vf_file', default=None, help='LRN RISK virulence factors TSV file')
|
|
12
|
952 parser.add_argument('--minimap2_version', action='store', dest='minimap2_version', default=None, help='minimap2 version string')
|
|
0
|
953 parser.add_argument('--mutation_regions_bed_file', action='store', dest='mutation_regions_bed_file', help='AMR mutation regions BRD file')
|
|
|
954 parser.add_argument('--mutation_regions_dir', action='store', dest='mutation_regions_dir', help='Directory of mutation regions TSV files')
|
|
26
|
955 parser.add_argument('--ont_file', action='store', dest='ont_file', help='ONT single read file')
|
|
0
|
956 parser.add_argument('--pima_css', action='store', dest='pima_css', help='PIMA css stypesheet')
|
|
23
|
957 parser.add_argument('--plasmids_file', action='store', dest='plasmids_file', default=None, help='pChunks plasmids TSV file')
|
|
13
|
958 parser.add_argument('--quast_report_file', action='store', dest='quast_report_file', help='Quast report tabular file')
|
|
18
|
959 parser.add_argument('--read_type', action='store', dest='read_type', help='Sample read type (ONT or Illumina)')
|
|
1
|
960 parser.add_argument('--reference_insertions_file', action='store', dest='reference_insertions_file', help='Reference insertions BED file')
|
|
12
|
961 parser.add_argument('--samtools_version', action='store', dest='samtools_version', default=None, help='Samtools version string')
|
|
|
962 parser.add_argument('--varscan_version', action='store', dest='varscan_version', default=None, help='Varscan version string')
|
|
0
|
963
|
|
|
964 args = parser.parse_args()
|
|
|
965
|
|
1
|
966 # Prepare the AMR matrix PNG files.
|
|
|
967 amr_matrix_files = []
|
|
|
968 for file_name in sorted(os.listdir(args.amr_matrix_png_dir)):
|
|
|
969 file_path = os.path.abspath(os.path.join(args.amr_matrix_png_dir, file_name))
|
|
|
970 amr_matrix_files.append(file_path)
|
|
13
|
971 # Prepare the circos PNG files.
|
|
|
972 circos_files = []
|
|
|
973 for file_name in sorted(os.listdir(args.circos_png_dir)):
|
|
|
974 file_path = os.path.abspath(os.path.join(args.circos_png_dir, file_name))
|
|
|
975 circos_files.append(file_path)
|
|
0
|
976 # Prepare the features BED files.
|
|
|
977 feature_bed_files = []
|
|
|
978 for file_name in sorted(os.listdir(args.feature_bed_dir)):
|
|
|
979 file_path = os.path.abspath(os.path.join(args.feature_bed_dir, file_name))
|
|
|
980 feature_bed_files.append(file_path)
|
|
|
981 # Prepare the features PNG files.
|
|
|
982 feature_png_files = []
|
|
|
983 for file_name in sorted(os.listdir(args.feature_png_dir)):
|
|
|
984 file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name))
|
|
|
985 feature_png_files.append(file_path)
|
|
|
986 # Prepare the mutation regions TSV files.
|
|
|
987 mutation_regions_files = []
|
|
|
988 for file_name in sorted(os.listdir(args.mutation_regions_dir)):
|
|
|
989 file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name))
|
|
|
990 mutation_regions_files.append(file_path)
|
|
|
991
|
|
|
992 markdown_report = PimaReport(args.analysis_name,
|
|
1
|
993 args.amr_deletions_file,
|
|
|
994 amr_matrix_files,
|
|
21
|
995 args.assembler_version,
|
|
0
|
996 args.assembly_fasta_file,
|
|
|
997 args.assembly_name,
|
|
12
|
998 args.bedtools_version,
|
|
2
|
999 args.blastn_version,
|
|
13
|
1000 circos_files,
|
|
1
|
1001 args.compute_sequence_length_file,
|
|
|
1002 args.contig_coverage_file,
|
|
|
1003 args.dbkey,
|
|
|
1004 args.dnadiff_snps_file,
|
|
2
|
1005 args.dnadiff_version,
|
|
18
|
1006 args.errors_file,
|
|
0
|
1007 feature_bed_files,
|
|
|
1008 feature_png_files,
|
|
1
|
1009 args.flye_assembly_info_file,
|
|
|
1010 args.genome_insertions_file,
|
|
0
|
1011 args.gzipped,
|
|
26
|
1012 args.illumina_forward_read_file,
|
|
|
1013 args.illumina_reverse_read_file,
|
|
2
|
1014 args.kraken2_report_file,
|
|
|
1015 args.kraken2_version,
|
|
28
|
1016 args.lrn_risk_amr_file,
|
|
|
1017 args.lrn_risk_blacklist_file,
|
|
|
1018 args.lrn_risk_vf_file,
|
|
12
|
1019 args.minimap2_version,
|
|
0
|
1020 args.mutation_regions_bed_file,
|
|
|
1021 mutation_regions_files,
|
|
26
|
1022 args.ont_file,
|
|
1
|
1023 args.pima_css,
|
|
|
1024 args.plasmids_file,
|
|
13
|
1025 args.quast_report_file,
|
|
18
|
1026 args.read_type,
|
|
12
|
1027 args.reference_insertions_file,
|
|
|
1028 args.samtools_version,
|
|
|
1029 args.varscan_version)
|
|
0
|
1030 markdown_report.make_report()
|
