0
+ − 1 import argparse
+ − 2 import os
+ − 3 import pandas
+ − 4 import pypandoc
+ − 5 import re
+ − 6 import subprocess
+ − 7 import sys
+ − 8
+ − 9 from Bio import SeqIO
+ − 10 from datetime import date
+ − 11 from mdutils.mdutils import MdUtils
+ − 12
+ − 13 CDC_ADVISORY = 'The analysis and report presented here should be treated as preliminary. Please contact the CDC/BDRD with any results regarding _Bacillus anthracis_.'
+ − 14
+ − 15
+ − 16 class PimaReport:
+ − 17
1
+ − 18 def __init__(self, analysis_name=None, amr_deletions_file=None, amr_matrix_files=None, assembly_fasta_file=None,
+ − 19 assembly_name=None, compute_sequence_length_file=None, contig_coverage_file=None, dbkey=None,
+ − 20 dnadiff_snps_file=None, feature_bed_files=None, feature_png_files=None, flye_assembly_info_file=None,
+ − 21 flye_version=None, genome_insertions_file=None, gzipped=None, illumina_fastq_file=None,
+ − 22 mutation_regions_bed_file=None, mutation_regions_tsv_files=None, pima_css=None, plasmids_file=None,
+ − 23 reference_insertions_file=None):
0
+ − 24 self.ofh = open("process_log.txt", "w")
+ − 25
1
+ − 26 self.ofh.write("amr_deletions_file: %s\n" % str(amr_deletions_file))
+ − 27 self.ofh.write("amr_matrix_files: %s\n" % str(amr_matrix_files))
0
+ − 28 self.ofh.write("analysis_name: %s\n" % str(analysis_name))
+ − 29 self.ofh.write("assembly_fasta_file: %s\n" % str(assembly_fasta_file))
+ − 30 self.ofh.write("assembly_name: %s\n" % str(assembly_name))
1
+ − 31 self.ofh.write("compute_sequence_length_file: %s\n" % str(compute_sequence_length_file))
+ − 32 self.ofh.write("contig_coverage_file: %s\n" % str(contig_coverage_file))
+ − 33 self.ofh.write("dbkey: %s\n" % str(dbkey))
+ − 34 self.ofh.write("dnadiff_snps_file: %s\n" % str(dnadiff_snps_file))
0
+ − 35 self.ofh.write("feature_bed_files: %s\n" % str(feature_bed_files))
+ − 36 self.ofh.write("feature_png_files: %s\n" % str(feature_png_files))
1
+ − 37 self.ofh.write("flye_assembly_info_file: %s\n" % str(flye_assembly_info_file))
+ − 38 self.ofh.write("flye_version: %s\n" % str(flye_version))
0
+ − 39 self.ofh.write("gzipped: %s\n" % str(gzipped))
1
+ − 40 self.ofh.write("genome_insertions_file: %s\n" % str(genome_insertions_file))
0
+ − 41 self.ofh.write("illumina_fastq_file: %s\n" % str(illumina_fastq_file))
+ − 42 self.ofh.write("mutation_regions_bed_file: %s\n" % str(mutation_regions_bed_file))
+ − 43 self.ofh.write("mutation_regions_tsv_files: %s\n" % str(mutation_regions_tsv_files))
+ − 44 self.ofh.write("pima_css: %s\n" % str(pima_css))
1
+ − 45 self.ofh.write("plasmids_file: %s\n" % str(plasmids_file))
+ − 46 # self.ofh.write("reference_genome: %s\n" % str(reference_genome))
+ − 47 self.ofh.write("reference_insertions_file: %s\n" % str(reference_insertions_file))
0
+ − 48
+ − 49 # General
+ − 50 self.doc = None
+ − 51 self.report_md = 'pima_report.md'
+ − 52
+ − 53 # Inputs
1
+ − 54 self.amr_deletions_file = amr_deletions_file
+ − 55 self.amr_matrix_files = amr_matrix_files
0
+ − 56 self.analysis_name = analysis_name
+ − 57 self.assembly_fasta_file = assembly_fasta_file
+ − 58 self.assembly_name = assembly_name
1
+ − 59 self.compute_sequence_length_file = compute_sequence_length_file
+ − 60 self.contig_coverage_file = contig_coverage_file
+ − 61 self.dbkey = dbkey
+ − 62 self.dnadiff_snps_file = dnadiff_snps_file
0
+ − 63 self.feature_bed_files = feature_bed_files
+ − 64 self.feature_png_files = feature_png_files
1
+ − 65 self.flye_assembly_info_file = flye_assembly_info_file
+ − 66 self.flye_version = flye_version
0
+ − 67 self.gzipped = gzipped
1
+ − 68 self.genome_insertions_file = genome_insertions_file
0
+ − 69 self.illumina_fastq_file = illumina_fastq_file
+ − 70 self.mutation_regions_bed_file = mutation_regions_bed_file
+ − 71 self.mutation_regions_tsv_files = mutation_regions_tsv_files
+ − 72 self.read_type = 'Illumina'
+ − 73 self.ont_bases = None
+ − 74 self.ont_n50 = None
+ − 75 self.ont_read_count = None
+ − 76 self.pima_css = pima_css
1
+ − 77 self.plasmids_file = plasmids_file
+ − 78 # self.reference_genome = reference_genome
+ − 79 self.reference_insertions_file = reference_insertions_file
0
+ − 80
+ − 81 # Titles
+ − 82 self.alignment_title = 'Comparison with reference'
+ − 83 self.alignment_notes_title = 'Alignment notes'
+ − 84 self.amr_matrix_title = 'AMR matrix'
+ − 85 self.assembly_methods_title = 'Assembly'
+ − 86 self.assembly_notes_title = 'Assembly notes'
+ − 87 self.basecalling_title = 'Basecalling'
+ − 88 self.basecalling_methods_title = 'Basecalling'
+ − 89 self.contamination_methods_title = 'Contamination check'
+ − 90 self.contig_alignment_title = 'Alignment vs. reference contigs'
+ − 91 self.feature_title = 'Features found in the assembly'
+ − 92 self.feature_methods_title = 'Feature annotation'
+ − 93 self.feature_plot_title = 'Feature annotation plots'
+ − 94 self.large_indel_title = 'Large insertions & deletions'
+ − 95 self.methods_title = 'Methods'
+ − 96 self.mutation_title = 'Mutations found in the sample'
+ − 97 self.mutation_methods_title = 'Mutation screening'
+ − 98 self.plasmid_methods_title = 'Plasmid annotation'
1
+ − 99 self.plasmid_title = 'Plasmid annotation'
0
+ − 100 self.reference_methods_title = 'Reference comparison'
+ − 101 self.snp_indel_title = 'SNPs and small indels'
+ − 102 self.summary_title = 'Analysis of %s' % analysis_name
+ − 103
+ − 104 # Methods
+ − 105 self.methods = pandas.Series(dtype='float64')
+ − 106 self.methods[self.contamination_methods_title] = pandas.Series(dtype='float64')
+ − 107 self.methods[self.assembly_methods_title] = pandas.Series(dtype='float64')
+ − 108 self.methods[self.reference_methods_title] = pandas.Series(dtype='float64')
+ − 109 self.methods[self.mutation_methods_title] = pandas.Series(dtype='float64')
+ − 110 self.methods[self.feature_methods_title] = pandas.Series(dtype='float64')
+ − 111 self.methods[self.plasmid_methods_title] = pandas.Series(dtype='float64')
+ − 112
+ − 113 # Contamination
+ − 114 self.kraken_fracs = pandas.Series(dtype=object)
+ − 115
+ − 116 # Notes
+ − 117 self.assembly_notes = pandas.Series(dtype=object)
+ − 118 self.alignment_notes = pandas.Series(dtype=object)
+ − 119 self.contig_alignment = pandas.Series(dtype=object)
+ − 120
+ − 121 # Values
+ − 122 self.assembly_size = 0
+ − 123 self.contig_info = None
+ − 124 self.did_medaka_ont_assembly = False
+ − 125 self.feature_hits = pandas.Series(dtype='float64')
+ − 126 self.illumina_length_mean = 0
+ − 127 self.illumina_read_count = 0
+ − 128 self.illumina_bases = 0
+ − 129 self.mean_coverage = 0
+ − 130 self.num_assembly_contigs = 0
1
+ − 131 # TODO: should the following 2 values be passed as parameters?
+ − 132 self.ont_n50_min = 2500
+ − 133 self.ont_coverage_min = 30
+ − 134 self.quast_indels = 0
+ − 135 self.quast_mismatches = 0
0
+ − 136
+ − 137 # Actions
+ − 138 self.did_guppy_ont_fast5 = False
+ − 139 self.did_qcat_ont_fastq = False
+ − 140 self.info_illumina_fastq()
+ − 141 self.load_contig_info()
+ − 142
+ − 143 def run_command(self, command):
+ − 144 self.ofh.write("\nXXXXXX In run_command, command:\n%s\n\n" % str(command))
+ − 145 try:
+ − 146 return re.split('\\n', subprocess.check_output(command, shell=True).decode('utf-8'))
+ − 147 except Exception:
+ − 148 message = 'Command %s failed: exiting...' % command
+ − 149 sys.exit(message)
+ − 150
+ − 151 def format_kmg(self, number, decimals=0):
+ − 152 self.ofh.write("\nXXXXXX In format_kmg, number:\n%s\n" % str(number))
+ − 153 self.ofh.write("XXXXXX In format_kmg, decimals:\n%s\n\n" % str(decimals))
+ − 154 if number == 0:
+ − 155 return '0'
+ − 156 magnitude_powers = [10**9, 10**6, 10**3, 1]
+ − 157 magnitude_units = ['G', 'M', 'K', '']
+ − 158 for i in range(len(magnitude_units)):
+ − 159 if number >= magnitude_powers[i]:
+ − 160 magnitude_power = magnitude_powers[i]
+ − 161 magnitude_unit = magnitude_units[i]
+ − 162 return ('{:0.' + str(decimals) + 'f}').format(number / magnitude_power) + magnitude_unit
+ − 163
+ − 164 def load_contig_info(self):
+ − 165 self.contig_info = pandas.Series(dtype=object)
+ − 166 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)
+ − 167 self.contig_info[self.read_type].columns = ['contig', 'size', 'coverage']
+ − 168 self.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()
1
+ − 169 if self.mean_coverage <= self.ont_coverage_min:
+ − 170 warning = '%s mean coverage ({:.0f}X) is less than the recommended minimum ({:.0f}X).'.format(self.mean_coverage, self.ont_coverage_min) % self.read_type
+ − 171 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+ − 172 # Report if some contigs have low coverage.
+ − 173 low_coverage = self.contig_info[self.read_type].loc[self.contig_info[self.read_type]['coverage'] < self.ont_coverage_min, :]
+ − 174 if low_coverage.shape[0] >= 0:
+ − 175 for contig_i in range(low_coverage.shape[0]):
+ − 176 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
+ − 177 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+ − 178 # See if some contigs have anolously low coverage.
+ − 179 fold_coverage = self.contig_info[self.read_type]['coverage'] / self.mean_coverage
+ − 180 low_coverage = self.contig_info[self.read_type].loc[fold_coverage < 1 / 5, :]
+ − 181 if low_coverage.shape[0] >= 0 :
+ − 182 for contig_i in range(low_coverage.shape[0]):
+ − 183 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], self.mean_coverage) % self.read_type
+ − 184 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
0
+ − 185
+ − 186 def load_fasta(self, fasta):
+ − 187 sequence = pandas.Series(dtype=object)
+ − 188 for contig in SeqIO.parse(fasta, 'fasta'):
+ − 189 sequence[contig.id] = contig
+ − 190 return sequence
+ − 191
+ − 192 def load_assembly(self):
+ − 193 self.assembly = self.load_fasta(self.assembly_fasta_file)
+ − 194 self.num_assembly_contigs = len(self.assembly)
+ − 195 for i in self.assembly:
+ − 196 self.assembly_size += len(i.seq)
+ − 197 self.assembly_size = self.format_kmg(self.assembly_size, decimals=1)
+ − 198
+ − 199 def info_illumina_fastq(self):
+ − 200 self.ofh.write("\nXXXXXX In info_illumina_fastq\n\n")
+ − 201 if self.gzipped:
+ − 202 opener = 'gunzip -c'
+ − 203 else:
+ − 204 opener = 'cat'
+ − 205 command = ' '.join([opener,
+ − 206 self.illumina_fastq_file,
+ − 207 '| awk \'{getline;s += length($1);getline;getline;}END{print s/(NR/4)"\t"(NR/4)"\t"s}\''])
+ − 208 output = self.run_command(command)
+ − 209 self.ofh.write("output:\n%s\n" % str(output))
+ − 210 self.ofh.write("re.split('\\t', self.run_command(command)[0]:\n%s\n" % str(re.split('\\t', self.run_command(command)[0])))
+ − 211 values = []
+ − 212 for i in re.split('\\t', self.run_command(command)[0]):
+ − 213 if i == '':
+ − 214 values.append(float('nan'))
+ − 215 else:
+ − 216 values.append(float(i))
+ − 217 self.ofh.write("values:\n%s\n" % str(values))
+ − 218 self.ofh.write("values[0]:\n%s\n" % str(values[0]))
+ − 219 self.illumina_length_mean += values[0]
+ − 220 self.ofh.write("values[1]:\n%s\n" % str(values[1]))
+ − 221 self.illumina_read_count += int(values[1])
+ − 222 self.ofh.write("values[2]:\n%s\n" % str(values[2]))
+ − 223 self.illumina_bases += int(values[2])
+ − 224 # The original PIMA code inserts self.illumina_fastq into
+ − 225 # a list for no apparent reason. We don't do that here.
+ − 226 # self.illumina_length_mean /= len(self.illumina_fastq)
+ − 227 self.illumina_length_mean /= 1
+ − 228 self.illumina_bases = self.format_kmg(self.illumina_bases, decimals=1)
+ − 229
+ − 230 def start_doc(self):
+ − 231 self.doc = MdUtils(file_name=self.report_md, title='')
+ − 232
+ − 233 def add_run_information(self):
+ − 234 self.ofh.write("\nXXXXXX In add_run_information\n\n")
+ − 235 self.doc.new_line()
+ − 236 self.doc.new_header(1, 'Run information')
+ − 237 # Tables in md.utils are implemented as a wrapping function.
+ − 238 Table_list = [
+ − 239 "Category",
+ − 240 "Information",
+ − 241 "Date",
+ − 242 date.today(),
+ − 243 "ONT FAST5",
+ − 244 "N/A",
+ − 245 "ONT FASTQ",
+ − 246 "N/A",
+ − 247 "Illumina FASTQ",
+ − 248 self.wordwrap_markdown(self.analysis_name),
+ − 249 "Assembly",
+ − 250 self.wordwrap_markdown(self.assembly_name),
+ − 251 "Reference",
+ − 252 self.wordwrap_markdown(self.dbkey),
+ − 253 ]
+ − 254 self.doc.new_table(columns=2, rows=7, text=Table_list, text_align='left')
+ − 255 self.doc.new_line()
+ − 256 self.doc.new_line()
+ − 257
+ − 258 def add_ont_library_information(self):
+ − 259 self.ofh.write("\nXXXXXX In add_ont_library_information\n\n")
+ − 260 if self.ont_n50 is None:
+ − 261 return
+ − 262 self.doc.new_line()
+ − 263 self.doc.new_header(2, 'ONT library statistics')
+ − 264 Table_List = [
+ − 265 "Category",
+ − 266 "Quantity",
+ − 267 "ONT N50",
+ − 268 '{:,}'.format(self.ont_n50),
+ − 269 "ONT reads",
+ − 270 '{:,}'.format(self.ont_read_count),
+ − 271 "ONT bases",
+ − 272 '{:s}'.format(self.ont_bases),
+ − 273 "Illumina FASTQ",
+ − 274 self.wordwrap_markdown(self.illumina_fastq_file),
+ − 275 "Assembly",
+ − 276 self.wordwrap_markdown(self.assembly_name),
+ − 277 "Reference",
+ − 278 self.wordwrap_markdown(self.dbkey),
+ − 279 ]
+ − 280 self.doc.new_table(columns=2, rows=7, text=Table_List, text_align='left')
+ − 281 self.doc.new_line()
+ − 282
+ − 283 def add_illumina_library_information(self):
+ − 284 self.ofh.write("\nXXXXXX In add_illumina_library_information\n\n")
+ − 285 if self.illumina_length_mean is None:
+ − 286 return
+ − 287 self.doc.new_line()
+ − 288 self.doc.new_header(2, 'Illumina library statistics')
+ − 289 Table_List = [
+ − 290 "Illumina Info.",
+ − 291 "Quantity",
+ − 292 'Illumina mean length',
+ − 293 '{:.1f}'.format(self.illumina_length_mean),
+ − 294 'Illumina reads',
+ − 295 '{:,}'.format(self.illumina_read_count),
+ − 296 'Illumina bases',
+ − 297 '{:s}'.format(self.illumina_bases)
+ − 298 ]
+ − 299 self.doc.new_table(columns=2, rows=4, text=Table_List, text_align='left')
+ − 300
1
+ − 301 def evaluate_assembly(self) :
+ − 302 assembly_info = pandas.read_csv(self.compute_sequence_length_file, sep='\t', header=None)
+ − 303 assembly_info.columns = ['contig', 'length']
+ − 304 self.contig_sizes = assembly_info
+ − 305 # Take a look at the number of contigs, their sizes,
+ − 306 # and circularity. Warn if things don't look good.
+ − 307 if assembly_info.shape[0] > 4:
+ − 308 warning = 'Assembly produced {:d} contigs, more than ususally expected; assembly may be fragmented'.format(assembly_info.shape[0])
+ − 309 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+ − 310 small_contigs = assembly_info.loc[assembly_info['length'] <= 3000, :]
+ − 311 if small_contigs.shape[0] > 0:
+ − 312 warning = 'Assembly produced {:d} small contigs ({:s}); assembly may include spurious sequences.'.format(small_contigs.shape[0], ', '.join(small_contigs['contig']))
+ − 313 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+ − 314
0
+ − 315 def add_assembly_information(self):
+ − 316 self.ofh.write("\nXXXXXX In add_assembly_information\n\n")
+ − 317 if self.assembly_fasta_file is None:
+ − 318 return
+ − 319 self.load_assembly()
+ − 320 self.doc.new_line()
+ − 321 self.doc.new_header(2, 'Assembly statistics')
+ − 322 Table_List = [
+ − 323 "Category",
+ − 324 "Information",
+ − 325 "Contigs",
+ − 326 str(self.num_assembly_contigs),
+ − 327 "Assembly size",
+ − 328 str(self.assembly_size),
+ − 329 ]
+ − 330 self.doc.new_table(columns=2, rows=3, text=Table_List, text_align='left')
+ − 331
+ − 332 def info_ont_fastq(self, fastq_file):
+ − 333 self.ofh.write("\nXXXXXX In info_ont_fastq, fastq_file:\n%s\n\n" % str(fastq_file))
+ − 334 opener = 'cat'
+ − 335 if self.gzipped:
+ − 336 opener = 'gunzip -c'
+ − 337 else:
+ − 338 opener = 'cat'
+ − 339 command = ' '.join([opener,
+ − 340 fastq_file,
+ − 341 '| awk \'{getline;print length($0);s += length($1);getline;getline;}END{print "+"s}\'',
+ − 342 '| sort -gr',
+ − 343 '| awk \'BEGIN{bp = 0;f = 0}',
+ − 344 '{if(NR == 1){sub(/+/, "", $1);s=$1}else{bp += $1;if(bp > s / 2 && f == 0){n50 = $1;f = 1}}}',
+ − 345 'END{printf "%d\\t%d\\t%d\\n", n50, (NR - 1), s;exit}\''])
+ − 346 result = list(re.split('\\t', self.run_command(command)[0]))
+ − 347 if result[1] == '0':
+ − 348 self.error_out('No ONT reads found')
+ − 349 ont_n50, ont_read_count, ont_raw_bases = [int(i) for i in result]
+ − 350 command = ' '.join([opener,
+ − 351 fastq_file,
+ − 352 '| awk \'{getline;print length($0);getline;getline;}\''])
+ − 353 result = self.run_command(command)
+ − 354 result = list(filter(lambda x: x != '', result))
1
+ − 355 # TODO: the following are not yet used...
+ − 356 # ont_read_lengths = [int(i) for i in result]
+ − 357 # ont_bases = self.format_kmg(ont_raw_bases, decimals=1)
+ − 358 if ont_n50 <= self.ont_n50_min:
+ − 359 warning = 'ONT N50 (%s) is less than the recommended minimum (%s)' % (str(ont_n50), str(self.ont_n50_min))
+ − 360 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
0
+ − 361
+ − 362 def wordwrap_markdown(self, string):
+ − 363 if string:
+ − 364 if len(string) < 35:
+ − 365 return string
+ − 366 else:
+ − 367 if '/' in string:
+ − 368 adjust = string.split('/')
+ − 369 out = ''
+ − 370 max = 35
+ − 371 for i in adjust:
+ − 372 out = out + '/' + i
+ − 373 if len(out) > max:
+ − 374 out += '<br>'
+ − 375 max += 35
+ − 376 return out
+ − 377 else:
+ − 378 out = [string[i:i + 35] for i in range(0, len(string), 50)]
+ − 379 return '<br>'.join(out)
+ − 380 else:
+ − 381 return string
+ − 382
+ − 383 def add_contig_info(self):
+ − 384 self.ofh.write("\nXXXXXX In add_contig_info\n\n")
+ − 385 if self.contig_info is None:
+ − 386 return
+ − 387 for method in ['ONT', 'Illumina']:
+ − 388 if method not in self.contig_info.index:
+ − 389 continue
+ − 390 self.doc.new_line()
+ − 391 self.doc.new_header(2, 'Assembly coverage by ' + method)
+ − 392 Table_List = ["Contig", "Length (bp)", "Coverage (X)"]
+ − 393 formatted = self.contig_info[method].copy()
+ − 394 formatted.iloc[:, 1] = formatted.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
+ − 395 for i in range(self.contig_info[method].shape[0]):
+ − 396 Table_List = Table_List + formatted.iloc[i, :].values.tolist()
+ − 397 row_count = int(len(Table_List) / 3)
+ − 398 self.doc.new_table(columns=3, rows=row_count, text=Table_List, text_align='left')
+ − 399
+ − 400 def add_assembly_notes(self):
+ − 401 self.ofh.write("\nXXXXXX In add_assembly_notes\n\n")
+ − 402 if len(self.assembly_notes) == 0:
+ − 403 return
+ − 404 self.doc.new_line()
+ − 405 self.doc.new_line('<div style="page-break-after: always;"></div>')
+ − 406 self.doc.new_line()
+ − 407 self.doc.new_header(2, self.assembly_notes_title)
1
+ − 408 for note in self.assembly_notes:
+ − 409 self.doc.new_line(note)
0
+ − 410
+ − 411 def add_contamination(self):
+ − 412 self.ofh.write("\nXXXXXX In add_contamination\n\n")
+ − 413 if len(self.kraken_fracs) == 0:
+ − 414 return
+ − 415 self.doc.new_line()
+ − 416 self.doc.new_header(2, 'Contamination check')
+ − 417 for read_type, kraken_fracs in self.kraken_fracs.iteritems():
+ − 418 self.doc.new_line(read_type + ' classifications')
+ − 419 self.doc.new_line()
+ − 420 Table_List = ["Percent of Reads", "Reads", "Level", "Label"]
+ − 421 for index, row in kraken_fracs.iterrows():
+ − 422 Table_List = Table_List + row.tolist()
+ − 423 row_count = int(len(Table_List) / 4)
+ − 424 self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left')
+ − 425 if self.contamination_methods_title not in self.methods:
+ − 426 self.methods[self.contamination_methods_title] = ''
+ − 427 method = 'Kraken2 was used to assign the raw reads into taxa.'
+ − 428 self.methods[self.contamination_methods_title] = self.methods[self.contamination_methods_title].append(pandas.Series(method))
+ − 429
+ − 430 def add_alignment(self):
+ − 431 self.ofh.write("\nXXXXXX In add_alignment\n\n")
+ − 432 # TODO: implement the draw_circos function for this.
+ − 433 if len(self.contig_alignment) > 0:
+ − 434 alignments = self.contig_alignment
+ − 435 else:
+ − 436 return
+ − 437 self.doc.new_line()
+ − 438 self.doc.new_header(level=2, title=self.alignment_title)
+ − 439 self.doc.new_line()
+ − 440 self.doc.new_header(level=3, title=self.snp_indel_title)
+ − 441 Table_1 = [
+ − 442 "Category",
+ − 443 "Quantity",
+ − 444 'SNPs',
1
+ − 445 '{:,}'.format(self.quast_mismatches),
0
+ − 446 'Small indels',
1
+ − 447 '{:,}'.format(self.quast_indels)
0
+ − 448 ]
+ − 449 self.doc.new_table(columns=2, rows=3, text=Table_1, text_align='left')
+ − 450 self.doc.new_line('<div style="page-break-after: always;"></div>')
+ − 451 self.doc.new_line()
+ − 452 if len(self.alignment_notes) > 0:
+ − 453 self.doc.new_header(level=3, title=self.alignment_notes_title)
+ − 454 for note in self.alignment_notes:
+ − 455 self.doc.new_line(note)
+ − 456 for contig in alignments.index.tolist():
+ − 457 contig_title = 'Alignment to %s' % contig
+ − 458 image_png = alignments[contig]
+ − 459 self.doc.new_line()
+ − 460 self.doc.new_header(level=3, title=contig_title)
+ − 461 self.doc.new_line(self.doc.new_inline_image(text='contig_title', path=os.path.abspath(image_png)))
+ − 462 self.doc.new_line('<div style="page-break-after: always;"></div>')
+ − 463 self.doc.new_line()
+ − 464 method = 'The genome assembly was aligned against the reference sequencing using dnadiff.'
+ − 465 self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method))
+ − 466
+ − 467 def add_features(self):
+ − 468 self.ofh.write("\nXXXXXX In add_features\n\n")
+ − 469 if len(self.feature_bed_files) == 0:
+ − 470 return
+ − 471 for bbf in self.feature_bed_files:
+ − 472 if os.path.getsize(bbf) > 0:
+ − 473 best = pandas.read_csv(filepath_or_buffer=bbf, sep='\t', header=None)
+ − 474 self.feature_hits[os.path.basename(bbf)] = best
+ − 475 if len(self.feature_hits) == 0:
+ − 476 return
+ − 477 self.ofh.write("self.feature_hits: %s\n" % str(self.feature_hits))
+ − 478 self.doc.new_line()
+ − 479 self.doc.new_header(level=2, title=self.feature_title)
+ − 480 for feature_name in self.feature_hits.index.tolist():
+ − 481 self.ofh.write("feature_name: %s\n" % str(feature_name))
+ − 482 features = self.feature_hits[feature_name].copy()
+ − 483 self.ofh.write("features: %s\n" % str(features))
+ − 484 if features.shape[0] == 0:
+ − 485 continue
+ − 486 features.iloc[:, 1] = features.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
+ − 487 features.iloc[:, 2] = features.iloc[:, 2].apply(lambda x: '{:,}'.format(x))
+ − 488 self.doc.new_line()
+ − 489 self.doc.new_header(level=3, title=feature_name)
+ − 490 if (features.shape[0] == 0):
+ − 491 continue
+ − 492 for contig in pandas.unique(features.iloc[:, 0]):
+ − 493 self.ofh.write("contig: %s\n" % str(contig))
+ − 494 self.doc.new_line(contig)
+ − 495 contig_features = features.loc[(features.iloc[:, 0] == contig), :]
+ − 496 self.ofh.write("contig_features: %s\n" % str(contig_features))
+ − 497 Table_List = ['Start', 'Stop', 'Feature', 'Identity (%)', 'Strand']
+ − 498 for i in range(contig_features.shape[0]):
+ − 499 self.ofh.write("i: %s\n" % str(i))
+ − 500 feature = contig_features.iloc[i, :].copy(deep=True)
+ − 501 self.ofh.write("feature: %s\n" % str(feature))
+ − 502 feature[4] = '{:.3f}'.format(feature[4])
1
+ − 503 # FIXME: Uncommenting the following line (which is
+ − 504 # https://github.com/appliedbinf/pima_md/blob/main/MarkdownReport.py#L317)
+ − 505 # will cause the job to fail with this exception:
+ − 506 # ValueError: columns * rows is not equal to text length
+ − 507 # Table_List = Table_List + feature[1:].values.tolist()
0
+ − 508 self.ofh.write("Table_List: %s\n" % str(Table_List))
+ − 509 row_count = int(len(Table_List) / 5)
+ − 510 self.ofh.write("row_count: %s\n" % str(row_count))
+ − 511 self.doc.new_line()
1
+ − 512 self.ofh.write("Before new_table, len(Table_List):: %s\n" % str(len(Table_List)))
+ − 513 self.doc.new_table(columns=5, rows=row_count, text=Table_List, text_align='left')
0
+ − 514 blastn_version = 'The genome assembly was queried for features using blastn.'
+ − 515 bedtools_version = 'Feature hits were clustered using bedtools and the highest scoring hit for each cluster was reported.'
+ − 516 method = '%s %s' % (blastn_version, bedtools_version)
+ − 517 self.methods[self.feature_methods_title] = self.methods[self.feature_methods_title].append(pandas.Series(method))
+ − 518
+ − 519 def add_feature_plots(self):
+ − 520 self.ofh.write("\nXXXXXX In add_feature_plots\n\n")
+ − 521 if len(self.feature_png_files) == 0:
+ − 522 return
+ − 523 self.doc.new_line()
+ − 524 self.doc.new_header(level=2, title='Feature Plots')
+ − 525 self.doc.new_paragraph('Only contigs with features are shown')
+ − 526 for feature_png_file in self.feature_png_files:
+ − 527 self.doc.new_line(self.doc.new_inline_image(text='Analysis', path=os.path.abspath(feature_png_file)))
+ − 528
+ − 529 def add_mutations(self):
+ − 530 self.ofh.write("\nXXXXXX In add_mutations\n\n")
+ − 531 if len(self.mutation_regions_tsv_files) == 0:
+ − 532 return
1
+ − 533 try :
0
+ − 534 mutation_regions = pandas.read_csv(self.mutation_regions_bed_file, sep='\t', header=0, index_col=False)
+ − 535 except Exception:
+ − 536 # Likely an empty file.
+ − 537 return
1
+ − 538 # TODO: this is the only place where reference_genome is used,
+ − 539 # so I'm commenting it out for now. We need to confirm if these
+ − 540 # errors that require the reference genmoe being passed are necessary.
+ − 541 # If so, we'll need to implement data tables in this tool.
+ − 542 # Make sure that the positions in the BED file fall within
+ − 543 # the chromosomes provided in the reference sequence.
+ − 544 """
+ − 545 for mutation_region in range(mutation_regions.shape[0]):
+ − 546 mutation_region = mutation_regions.iloc[mutation_region, :]
+ − 547 if not (mutation_region[0] in self.reference_genome):
+ − 548 self.ofh.write("\nMutation region: %s not found in reference genome.\n" % ' '.join(mutation_region.astype(str)))
+ − 549 continue
+ − 550 if not isinstance(mutation_region[1], int):
+ − 551 self.ofh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(mutation_region[1]))
+ − 552 break
+ − 553 elif not isinstance(mutation_region[2], int):
+ − 554 self.ofh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(mutation_region[2]))
+ − 555 break
+ − 556 if mutation_region[1] <= 0 or mutation_region[2] <= 0:
+ − 557 self.ofh.write("\nMutation region %s starts before the reference sequence.\n" % ' '.join(mutation_region.astype(str)))
+ − 558 if mutation_region[1] > len(self.reference_genome[mutation_region[0]].seq) or mutation_region[2] > len(self.reference_genome[mutation_region[0]].seq):
+ − 559 self.ofh.write("\nMutation region %s ends after the reference sequence.\n" % ' '.join(mutation_region.astype(str)))
+ − 560 """
0
+ − 561 amr_mutations = pandas.Series(dtype=object)
+ − 562 for region_i in range(mutation_regions.shape[0]):
+ − 563 region = mutation_regions.iloc[region_i, :]
+ − 564 region_name = str(region['name'])
+ − 565 self.ofh.write("Processing mutations for region %s\n" % region_name)
+ − 566 region_mutations_tsv_name = '%s_mutations.tsv' % region_name
+ − 567 if region_mutations_tsv_name not in self.mutation_regions_tsv_files:
+ − 568 continue
+ − 569 region_mutations_tsv = self.mutation_regions_tsv_files[region_mutations_tsv_name]
1
+ − 570 try :
0
+ − 571 region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False)
+ − 572 except Exception:
+ − 573 region_mutations = pandas.DataFrame()
+ − 574 if region_mutations.shape[0] == 0:
+ − 575 continue
1
+ − 576 # Figure out what kind of mutations are in this region.
0
+ − 577 region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index)
+ − 578 region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel'
+ − 579 region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index)
+ − 580 region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index)
+ − 581 region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1)
+ − 582 region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']]
+ − 583 amr_mutations[region['name']] = region_mutations
+ − 584 # Report the mutations.
+ − 585 self.doc.new_line()
+ − 586 self.doc.new_header(level=2, title=self.mutation_title)
+ − 587 for region_name in amr_mutations.index.tolist():
+ − 588 region_mutations = amr_mutations[region_name].copy()
+ − 589 self.doc.new_line()
+ − 590 self.doc.new_header(level=3, title=region_name)
+ − 591 if (region_mutations.shape[0] == 0):
+ − 592 self.doc.append('None')
+ − 593 continue
+ − 594 region_mutations.iloc[:, 1] = region_mutations.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
+ − 595 Table_List = ['Reference contig', 'Position', 'Reference', 'Alternate', 'Drug', 'Note']
+ − 596 for i in range(region_mutations.shape[0]):
+ − 597 Table_List = Table_List + region_mutations.iloc[i, [0, 1, 3, 4, 5, 6]].values.tolist()
+ − 598 row_count = int(len(Table_List) / 6)
+ − 599 self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left')
+ − 600 method = '%s reads were mapped to the reference sequence using minimap2.' % self.read_type
+ − 601 self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
+ − 602 method = 'Mutations were identified using samtools mpileup and varscan.'
+ − 603 self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
+ − 604
+ − 605 def add_amr_matrix(self):
+ − 606 self.ofh.write("\nXXXXXX In add_amr_matrix\n\n")
+ − 607 # Make sure that we have an AMR matrix to plot
1
+ − 608 if len(self.amr_matrix_files) == 0:
+ − 609 return
+ − 610 self.doc.new_line()
+ − 611 self.doc.new_header(level=2, title=self.amr_matrix_title)
+ − 612 self.doc.new_line('AMR genes and mutations with their corresponding drugs')
+ − 613 for amr_matrix_file in self.amr_matrix_files:
+ − 614 self.doc.new_line(self.doc.new_inline_image(text='AMR genes and mutations with their corresponding drugs',
+ − 615 path=os.path.abspath(amr_matrix_file)))
0
+ − 616
+ − 617 def add_large_indels(self):
+ − 618 self.ofh.write("\nXXXXXX In add_large_indels\n\n")
1
+ − 619 large_indels = pandas.Series(dtype='float64')
+ − 620 # Pull in insertions.
+ − 621 try:
+ − 622 reference_insertions = pandas.read_csv(filepath_or_buffer=self.reference_insertions_file, sep='\t', header=None)
+ − 623 except Exception:
+ − 624 reference_insertions = pandas.DataFrame()
+ − 625 try:
+ − 626 genome_insertions = pandas.read_csv(filepath_or_buffer=self.genome_insertions_file, sep='\t', header=None)
+ − 627 except Exception:
+ − 628 genome_insertions = pandas.DataFrame()
+ − 629 large_indels['Reference insertions'] = reference_insertions
+ − 630 large_indels['Query insertions'] = genome_insertions
+ − 631 # TODO: we don't seem to be reporting snps and deletions for some reason...
+ − 632 # Pull in the number of SNPs and small indels.
+ − 633 try:
+ − 634 snps = pandas.read_csv(filepath_or_buffer=self.dnadiff_snps_file, sep='\t', header=None)
+ − 635 # TODO: the following is not used...
+ − 636 # small_indels = snps.loc[(snps.iloc[:, 1] == '.') | (snps.iloc[:, 2] == '.'), :]
+ − 637 snps = snps.loc[(snps.iloc[:, 1] != '.') & (snps.iloc[:, 2] != '.'), :]
+ − 638 except Exception:
+ − 639 snps = pandas.DataFrame()
+ − 640 # Pull in deletions.
+ − 641 try:
+ − 642 amr_deletions = pandas.read_csv(filepath_or_buffer=self.amr_deletion_file, sep='\t', header=None)
+ − 643 except Exception:
+ − 644 amr_deletions = pandas.DataFrame()
+ − 645 if amr_deletions.shape[0] > 0:
+ − 646 amr_deletions.columns = ['contig', 'start', 'stop', 'name', 'type', 'drug', 'note']
+ − 647 amr_deletions = amr_deletions.loc[amr_deletions['type'].isin(['large-deletion', 'any']), :]
+ − 648 self.doc.new_line()
+ − 649 self.doc.new_header(level=2, title=self.large_indel_title)
+ − 650 for genome in ['Reference insertions', 'Query insertions']:
+ − 651 genome_indels = large_indels[genome].copy()
+ − 652 self.doc.new_line()
+ − 653 self.doc.new_header(level=3, title=genome)
+ − 654 if (genome_indels.shape[0] == 0):
+ − 655 continue
+ − 656 genome_indels.iloc[:, 1] = genome_indels.iloc[:, 1].apply(lambda x: '{:,}'.format(x))
+ − 657 genome_indels.iloc[:, 2] = genome_indels.iloc[:, 2].apply(lambda x: '{:,}'.format(x))
+ − 658 genome_indels.iloc[:, 3] = genome_indels.iloc[:, 3].apply(lambda x: '{:,}'.format(x))
+ − 659 Table_List = [
+ − 660 'Reference contig', 'Start', 'Stop', 'Size (bp)'
+ − 661 ]
+ − 662 for i in range(genome_indels.shape[0]):
+ − 663 Table_List = Table_List + genome_indels.iloc[i, :].values.tolist()
+ − 664 row_count = int(len(Table_List) / 4)
+ − 665 self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left')
+ − 666 method = 'Large insertions or deletions were found as the complement of aligned regions using bedtools.'
+ − 667 self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method))
+ − 668 self.doc.new_line()
+ − 669 self.doc.new_line('<div style="page-break-after: always;"></div>')
+ − 670 self.doc.new_line()
0
+ − 671
+ − 672 def add_plasmids(self):
1
+ − 673 try :
+ − 674 plasmids = pandas.read_csv(filepath_or_buffer=self.plasmids_file, sep='\t', header=0)
+ − 675 except Exception:
0
+ − 676 return
+ − 677 plasmids = plasmids.copy()
+ − 678 self.doc.new_line()
1
+ − 679 self.doc.new_header(level=2, title=self.plasmid_title)
0
+ − 680 if (plasmids.shape[0] == 0):
+ − 681 self.doc.new_line('None')
+ − 682 return
+ − 683 plasmids.iloc[:, 3] = plasmids.iloc[:, 3].apply(lambda x: '{:,}'.format(x))
+ − 684 plasmids.iloc[:, 4] = plasmids.iloc[:, 4].apply(lambda x: '{:,}'.format(x))
+ − 685 plasmids.iloc[:, 5] = plasmids.iloc[:, 5].apply(lambda x: '{:,}'.format(x))
1
+ − 686 Table_List = ['Genome contig', 'Plasmid hit', 'Plasmid acc.', 'Contig size', 'Aliged', 'Plasmid size']
0
+ − 687 for i in range(plasmids.shape[0]):
+ − 688 Table_List = Table_List + plasmids.iloc[i, 0:6].values.tolist()
+ − 689 row_count = int(len(Table_List) / 6)
+ − 690 self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left')
+ − 691 method = 'The plasmid reference database was queried against the genome assembly using minimap2.'
+ − 692 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
+ − 693 method = 'The resulting SAM was converted to a PSL using a custom version of sam2psl.'
+ − 694 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
+ − 695 method = 'Plasmid-to-genome hits were resolved using the pChunks algorithm.'
+ − 696 self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method))
+ − 697
+ − 698 def add_methods(self):
+ − 699 self.ofh.write("\nXXXXXX In add_methods\n\n")
+ − 700 self.doc.new_line('<div style="page-break-after: always;"></div>')
+ − 701 self.doc.new_line()
+ − 702 if len(self.methods) == 0:
+ − 703 return
+ − 704 self.doc.new_line()
+ − 705 self.doc.new_header(level=2, title=self.methods_title)
+ − 706 for methods_section in self.methods.index.tolist():
+ − 707 if self.methods[methods_section] is None or len(self.methods[methods_section]) == 0:
+ − 708 continue
+ − 709 self.doc.new_line()
+ − 710 self.doc.new_header(level=3, title=methods_section)
+ − 711 self.doc.new_paragraph(' '.join(self.methods[methods_section]))
+ − 712
+ − 713 def add_summary(self):
+ − 714 self.ofh.write("\nXXXXXX In add_summary\n\n")
+ − 715 # Add summary title
+ − 716 self.doc.new_header(level=1, title=self.summary_title)
+ − 717 # First section of Summary
+ − 718 self.doc.new_header(level=1, title='CDC Advisory')
+ − 719 self.doc.new_paragraph(CDC_ADVISORY)
+ − 720 self.doc.new_line()
+ − 721 self.add_run_information()
+ − 722 self.add_ont_library_information()
+ − 723 methods = []
+ − 724 if self.did_guppy_ont_fast5:
+ − 725 methods += ['ONT reads were basecalled using guppy']
+ − 726 if self.did_qcat_ont_fastq:
+ − 727 methods += ['ONT reads were demultiplexed and trimmed using qcat']
+ − 728 self.methods[self.basecalling_methods_title] = pandas.Series(methods)
+ − 729 self.add_illumina_library_information()
1
+ − 730 self.add_contig_info()
+ − 731 self.evaluate_assembly()
0
+ − 732 self.add_assembly_information()
1
+ − 733 if self.flye_assembly_info_file is not None:
+ − 734 method = 'ONT reads were assembled using %s' % self.flye_version
0
+ − 735 self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.Series(method))
1
+ − 736 # Pull in the assembly summary and look at the coverage.
+ − 737 assembly_info = pandas.read_csv(self.flye_assembly_info_file, header=0, index_col=0, sep='\t')
+ − 738 # Look for non-circular contigs.
+ − 739 open_contigs = assembly_info.loc[assembly_info['circ.'] == 'N', :]
+ − 740 if open_contigs.shape[0] > 0:
+ − 741 open_contig_ids = open_contigs.index.values
+ − 742 warning = 'Flye reported {:d} open contigs ({:s}); assembly may be incomplete.'.format(open_contigs.shape[0], ', '.join(open_contig_ids))
+ − 743 self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
0
+ − 744 if self.did_medaka_ont_assembly:
+ − 745 method = 'the genome assembly was polished using ont reads and medaka.'
+ − 746 self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.series(method))
1
+ − 747 self.info_ont_fastq(self.illumina_fastq_file)
+ − 748 self.add_assembly_notes()
0
+ − 749
+ − 750 def make_tex(self):
+ − 751 self.doc.new_table_of_contents(table_title='detailed run information', depth=2, marker="tableofcontents")
+ − 752 text = self.doc.file_data_text
+ − 753 text = text.replace("##--[", "")
+ − 754 text = text.replace("]--##", "")
+ − 755 self.doc.file_data_text = text
+ − 756 self.doc.create_md_file()
+ − 757
+ − 758 def make_report(self):
+ − 759 self.ofh.write("\nXXXXXX In make_report\n\n")
+ − 760 self.start_doc()
+ − 761 self.add_summary()
+ − 762 self.add_contamination()
+ − 763 self.add_alignment()
+ − 764 self.add_features()
+ − 765 self.add_feature_plots()
+ − 766 self.add_mutations()
+ − 767 self.add_large_indels()
+ − 768 self.add_plasmids()
1
+ − 769 # TODO stuff working to here...
0
+ − 770 self.add_amr_matrix()
+ − 771 # self.add_snps()
+ − 772 self.add_methods()
+ − 773 self.make_tex()
+ − 774 # It took me quite a long time to find out that the value of the -t
+ − 775 # (implied) argument in the following command must be 'html' instead of
+ − 776 # the more logical 'pdf'. see the answer from snsn in this thread:
+ − 777 # https://github.com/jessicategner/pypandoc/issues/186
+ − 778 self.ofh.write("\nXXXXX In make_report, calling pypandoc.convert_file...\n\n")
+ − 779 pypandoc.convert_file(self.report_md,
+ − 780 'html',
+ − 781 extra_args=['--pdf-engine=weasyprint', '-V', '-css=%s' % self.pima_css],
+ − 782 outputfile='pima_report.pdf')
+ − 783 self.ofh.close()
+ − 784
+ − 785
+ − 786 parser = argparse.ArgumentParser()
+ − 787
1
+ − 788 parser.add_argument('--amr_deletions_file', action='store', dest='amr_deletions_file', help='AMR deletions BED file')
+ − 789 parser.add_argument('--amr_matrix_png_dir', action='store', dest='amr_matrix_png_dir', help='Directory of AMR matrix PNG files')
0
+ − 790 parser.add_argument('--analysis_name', action='store', dest='analysis_name', help='Sample identifier')
+ − 791 parser.add_argument('--assembly_fasta_file', action='store', dest='assembly_fasta_file', help='Assembly fasta file')
+ − 792 parser.add_argument('--assembly_name', action='store', dest='assembly_name', help='Assembly identifier')
1
+ − 793 parser.add_argument('--compute_sequence_length_file', action='store', dest='compute_sequence_length_file', help='Comnpute sequence length tabular file')
+ − 794 parser.add_argument('--contig_coverage_file', action='store', dest='contig_coverage_file', help='Contig coverage TSV file')
+ − 795 parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference genome identifier')
+ − 796 parser.add_argument('--dnadiff_snps_file', action='store', dest='dnadiff_snps_file', help='DNAdiff snps tabular file')
0
+ − 797 parser.add_argument('--feature_bed_dir', action='store', dest='feature_bed_dir', help='Directory of best feature hits bed files')
+ − 798 parser.add_argument('--feature_png_dir', action='store', dest='feature_png_dir', help='Directory of best feature hits png files')
1
+ − 799 parser.add_argument('--flye_assembly_info_file', action='store', dest='flye_assembly_info_file', default=None, help='Flye assembly info tabular file')
+ − 800 parser.add_argument('--flye_version', action='store', dest='flye_version', default=None, help='Flye version string')
+ − 801 parser.add_argument('--genome_insertions_file', action='store', dest='genome_insertions_file', help='Genome insertions BED file')
0
+ − 802 parser.add_argument('--gzipped', action='store_true', dest='gzipped', default=False, help='Input sample is gzipped')
+ − 803 parser.add_argument('--illumina_fastq_file', action='store', dest='illumina_fastq_file', help='Input sample')
+ − 804 parser.add_argument('--mutation_regions_bed_file', action='store', dest='mutation_regions_bed_file', help='AMR mutation regions BRD file')
+ − 805 parser.add_argument('--mutation_regions_dir', action='store', dest='mutation_regions_dir', help='Directory of mutation regions TSV files')
+ − 806 parser.add_argument('--pima_css', action='store', dest='pima_css', help='PIMA css stypesheet')
1
+ − 807 parser.add_argument('--plasmids_file', action='store', dest='plasmids_file', help='pChunks plasmids TSV file')
+ − 808 parser.add_argument('--reference_insertions_file', action='store', dest='reference_insertions_file', help='Reference insertions BED file')
+ − 809 # parser.add_argument('--reference_genome', action='store', dest='reference_genome', help='Reference genome fasta file')
0
+ − 810
+ − 811 args = parser.parse_args()
+ − 812
1
+ − 813 # Prepare the AMR matrix PNG files.
+ − 814 amr_matrix_files = []
+ − 815 for file_name in sorted(os.listdir(args.amr_matrix_png_dir)):
+ − 816 file_path = os.path.abspath(os.path.join(args.amr_matrix_png_dir, file_name))
+ − 817 amr_matrix_files.append(file_path)
0
+ − 818 # Prepare the features BED files.
+ − 819 feature_bed_files = []
+ − 820 for file_name in sorted(os.listdir(args.feature_bed_dir)):
+ − 821 file_path = os.path.abspath(os.path.join(args.feature_bed_dir, file_name))
+ − 822 feature_bed_files.append(file_path)
+ − 823 # Prepare the features PNG files.
+ − 824 feature_png_files = []
+ − 825 for file_name in sorted(os.listdir(args.feature_png_dir)):
+ − 826 file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name))
+ − 827 feature_png_files.append(file_path)
+ − 828 # Prepare the mutation regions TSV files.
+ − 829 mutation_regions_files = []
+ − 830 for file_name in sorted(os.listdir(args.mutation_regions_dir)):
+ − 831 file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name))
+ − 832 mutation_regions_files.append(file_path)
+ − 833
+ − 834 markdown_report = PimaReport(args.analysis_name,
1
+ − 835 args.amr_deletions_file,
+ − 836 amr_matrix_files,
0
+ − 837 args.assembly_fasta_file,
+ − 838 args.assembly_name,
1
+ − 839 args.compute_sequence_length_file,
+ − 840 args.contig_coverage_file,
+ − 841 args.dbkey,
+ − 842 args.dnadiff_snps_file,
0
+ − 843 feature_bed_files,
+ − 844 feature_png_files,
1
+ − 845 args.flye_assembly_info_file,
+ − 846 args.flye_version,
+ − 847 args.genome_insertions_file,
0
+ − 848 args.gzipped,
+ − 849 args.illumina_fastq_file,
+ − 850 args.mutation_regions_bed_file,
+ − 851 mutation_regions_files,
1
+ − 852 args.pima_css,
+ − 853 args.plasmids_file,
+ − 854 args.reference_insertions_file)
0
+ − 855 markdown_report.make_report()
