Mercurial > repos > greg > pima_report
changeset 0:0a558f444c98 draft
Uploaded
| author | greg |
|---|---|
| date | Fri, 03 Mar 2023 22:06:23 +0000 |
| parents | |
| children | 67d0939b56b0 |
| files | .shed.yml macros.xml pima.css pima_report.py pima_report.xml |
| diffstat | 5 files changed, 987 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.shed.yml Fri Mar 03 22:06:23 2023 +0000 @@ -0,0 +1,9 @@ +name: pima_report +owner: greg +description: Generates the PIMA analysis summary report +long_description: Generates the PIMA analysis summary report +categories: +- Nanopore +remote_repository_url: https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/pima/pima_report +homepage_url: https://github.com/gregvonkuster/galaxy_tools +type: unrestricted
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros.xml Fri Mar 03 22:06:23 2023 +0000 @@ -0,0 +1,22 @@ +<macros> + <token name="@TOOL_VERSION@">1.0.0</token> + <token name="@VERSION_SUFFIX@">0</token> + <token name="@PROFILE@">21.01</token> + <xml name="requirements"> + <requirements> + <requirement type="package" version="1.79">biopython</requirement> + <requirement type="package" version="5.1.0">gawk</requirement> + <requirement type="package" version="3.4">grep</requirement> + <requirement type="package" version="1.5.0">mdutils</requirement> + <requirement type="package" version="1.5.2">pandas</requirement> + <requirement type="package" version="1.10">pypandoc</requirement> + <requirement type="package" version="57.1">weasyprint</requirement> + </requirements> + </xml> + <xml name="citations"> + <citations> + <citation type="doi">10.1101/2020.01.09.900589</citation> + </citations> + </xml> +</macros> +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pima.css Fri Mar 03 22:06:23 2023 +0000 @@ -0,0 +1,143 @@ +html { + line-height: 1.5; + font-family: Georgia, serif; + font-size: 20px; + color: #1a1a1a; + background-color: #fdfdfd; +} +body { + margin: 0 auto; + max-width: 50em; + padding-left: 20px; + padding-right: 20px; + padding-top: 20px; + padding-bottom: 20px; + hyphens: auto; + overflow-wrap: break-word; + font-kerning: normal; +} +@media print { + body { + background-color: transparent; + color: black; + font-size: 10pt; + } + p, h2, h3 { + orphans: 3; + widows: 3; + } + h2, h3, h4 { + page-break-after: avoid; + } +} +p { + margin: 1em 0; +} +a { + color: #3333FF; +} +a:visited { + color: #1a1a1a; +} +img { + max-width: 100%; +} +h1, h2, h3, h4, h5, h6 { + margin-top: 1.4em; +} +h5, h6 { + font-size: .8em; + font-style: italic; +} +h6 { + font-weight: normal; +} +ol, ul { + padding-left: 1.7em; + margin-top: 1em; +} +li > ol, li > ul { + margin-top: 0; +} +blockquote { + margin: 1em 0 1em 1.7em; + padding-left: 1em; + border-left: 2px solid #e6e6e6; + color: #606060; +} +code { + font-family: Menlo, Monaco, 'Lucida Console', Consolas, monospace; + font-size: 85%; + margin: 0; +} +pre { + margin: 1em 0; + overflow: auto; +} +pre code { + padding: 0; + overflow: visible; + overflow-wrap: normal; +} +.sourceCode { + background-color: transparent; + overflow: visible; +} +hr { + background-color: #1a1a1a; + border: none; + height: 1px; + margin: 1em 0; +} +table { + border-collapse: collapse; + font-variant-numeric: lining-nums tabular-nums; + /* + margin-left: auto; + margin-right: auto; */ +} + +table caption { + margin-bottom: 0.75em; +} +tbody { + margin-top: 0.5em; + border-top: 1px solid #1a1a1a; + border-bottom: 1px solid #1a1a1a; + text-align: center; +} +th { + border-top: 1px solid #1a1a1a; + padding: 0.25em 0.5em 0.25em 0.5em; +} +td { + padding: 0.125em 0.5em 0.25em 0.5em; + +} + +tr:nth-child(even) {background: #CCC} +tr:nth-child(odd) {background: #FFF} + +header { + margin-bottom: 4em; + text-align: center; +} +#TOC li { + list-style: none; +} +#TOC ul { + padding-left: 1.3em; +} +#TOC > ul { + padding-left: 0; +} +#TOC a:not(:hover) { + text-decoration: none; +} +code{white-space: pre-wrap;} +span.smallcaps{font-variant: small-caps;} +span.underline{text-decoration: underline;} +div.column{display: inline-block; vertical-align: top; width: 50%;} +div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;} +ul.task-list{list-style: none;} +.display.math{display: block; text-align: center; margin: 0.5rem auto;}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pima_report.py Fri Mar 03 22:06:23 2023 +0000 @@ -0,0 +1,733 @@ +import argparse +import os +import pandas +import pypandoc +import re +import subprocess +import sys + +from Bio import SeqIO +from datetime import date +from mdutils.mdutils import MdUtils + +CDC_ADVISORY = 'The analysis and report presented here should be treated as preliminary. Please contact the CDC/BDRD with any results regarding _Bacillus anthracis_.' + + +class PimaReport: + + def __init__(self, analysis_name=None, assembly_fasta_file=None, assembly_name=None, feature_bed_files=None, feature_png_files=None, + contig_coverage_file=None, dbkey=None, gzipped=None, illumina_fastq_file=None, mutation_regions_bed_file=None, + mutation_regions_tsv_files=None, pima_css=None): + self.ofh = open("process_log.txt", "w") + + self.ofh.write("analysis_name: %s\n" % str(analysis_name)) + self.ofh.write("assembly_fasta_file: %s\n" % str(assembly_fasta_file)) + self.ofh.write("assembly_name: %s\n" % str(assembly_name)) + self.ofh.write("feature_bed_files: %s\n" % str(feature_bed_files)) + self.ofh.write("feature_png_files: %s\n" % str(feature_png_files)) + self.ofh.write("contig_coverage_file: %s\n" % str(contig_coverage_file)) + self.ofh.write("dbkey: %s\n" % str(dbkey)) + self.ofh.write("gzipped: %s\n" % str(gzipped)) + self.ofh.write("illumina_fastq_file: %s\n" % str(illumina_fastq_file)) + self.ofh.write("mutation_regions_bed_file: %s\n" % str(mutation_regions_bed_file)) + self.ofh.write("mutation_regions_tsv_files: %s\n" % str(mutation_regions_tsv_files)) + self.ofh.write("pima_css: %s\n" % str(pima_css)) + + # General + self.doc = None + self.report_md = 'pima_report.md' + + # Inputs + self.analysis_name = analysis_name + self.assembly_fasta_file = assembly_fasta_file + self.assembly_name = assembly_name + self.feature_bed_files = feature_bed_files + self.feature_png_files = feature_png_files + self.contig_coverage_file = contig_coverage_file + self.dbkey = dbkey + self.gzipped = gzipped + self.illumina_fastq_file = illumina_fastq_file + self.mutation_regions_bed_file = mutation_regions_bed_file + self.mutation_regions_tsv_files = mutation_regions_tsv_files + self.read_type = 'Illumina' + self.ont_bases = None + self.ont_n50 = None + self.ont_read_count = None + self.pima_css = pima_css + + # Titles + self.alignment_title = 'Comparison with reference' + self.alignment_notes_title = 'Alignment notes' + self.amr_matrix_title = 'AMR matrix' + self.assembly_methods_title = 'Assembly' + self.assembly_notes_title = 'Assembly notes' + self.basecalling_title = 'Basecalling' + self.basecalling_methods_title = 'Basecalling' + self.contamination_methods_title = 'Contamination check' + self.contig_alignment_title = 'Alignment vs. reference contigs' + self.feature_title = 'Features found in the assembly' + self.feature_methods_title = 'Feature annotation' + self.feature_plot_title = 'Feature annotation plots' + self.large_indel_title = 'Large insertions & deletions' + self.methods_title = 'Methods' + self.mutation_title = 'Mutations found in the sample' + self.mutation_methods_title = 'Mutation screening' + self.plasmid_methods_title = 'Plasmid annotation' + self.reference_methods_title = 'Reference comparison' + self.snp_indel_title = 'SNPs and small indels' + #self.summary_title = 'Summary' + self.summary_title = 'Analysis of %s' % analysis_name + + # Methods + self.methods = pandas.Series(dtype='float64') + self.methods[self.contamination_methods_title] = pandas.Series(dtype='float64') + self.methods[self.assembly_methods_title] = pandas.Series(dtype='float64') + self.methods[self.reference_methods_title] = pandas.Series(dtype='float64') + self.methods[self.mutation_methods_title] = pandas.Series(dtype='float64') + self.methods[self.feature_methods_title] = pandas.Series(dtype='float64') + self.methods[self.plasmid_methods_title] = pandas.Series(dtype='float64') + + # Contamination + self.kraken_fracs = pandas.Series(dtype=object) + + # Notes + self.assembly_notes = pandas.Series(dtype=object) + self.alignment_notes = pandas.Series(dtype=object) + self.contig_alignment = pandas.Series(dtype=object) + + # Values + self.assembly_size = 0 + self.contig_info = None + self.did_flye_ont_fastq = False + self.did_medaka_ont_assembly = False + self.feature_hits = pandas.Series(dtype='float64') + self.illumina_length_mean = 0 + self.illumina_read_count = 0 + self.illumina_bases = 0 + self.mean_coverage = 0 + self.num_assembly_contigs = 0 + + # Actions + self.did_guppy_ont_fast5 = False + self.did_qcat_ont_fastq = False + self.info_illumina_fastq() + self.load_contig_info() + + def run_command(self, command): + self.ofh.write("\nXXXXXX In run_command, command:\n%s\n\n" % str(command)) + try: + return re.split('\\n', subprocess.check_output(command, shell=True).decode('utf-8')) + except Exception: + message = 'Command %s failed: exiting...' % command + sys.exit(message) + + def format_kmg(self, number, decimals=0): + self.ofh.write("\nXXXXXX In format_kmg, number:\n%s\n" % str(number)) + self.ofh.write("XXXXXX In format_kmg, decimals:\n%s\n\n" % str(decimals)) + if number == 0: + return '0' + magnitude_powers = [10**9, 10**6, 10**3, 1] + magnitude_units = ['G', 'M', 'K', ''] + for i in range(len(magnitude_units)): + if number >= magnitude_powers[i]: + magnitude_power = magnitude_powers[i] + magnitude_unit = magnitude_units[i] + return ('{:0.' + str(decimals) + 'f}').format(number / magnitude_power) + magnitude_unit + + def load_contig_info(self): + self.contig_info = pandas.Series(dtype=object) + 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) + self.contig_info[self.read_type].columns = ['contig', 'size', 'coverage'] + 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() + + def load_fasta(self, fasta): + sequence = pandas.Series(dtype=object) + for contig in SeqIO.parse(fasta, 'fasta'): + sequence[contig.id] = contig + return sequence + + def load_assembly(self): + self.assembly = self.load_fasta(self.assembly_fasta_file) + self.num_assembly_contigs = len(self.assembly) + for i in self.assembly: + self.assembly_size += len(i.seq) + self.assembly_size = self.format_kmg(self.assembly_size, decimals=1) + + def info_illumina_fastq(self): + self.ofh.write("\nXXXXXX In info_illumina_fastq\n\n") + if self.gzipped: + opener = 'gunzip -c' + else: + opener = 'cat' + command = ' '.join([opener, + self.illumina_fastq_file, + '| awk \'{getline;s += length($1);getline;getline;}END{print s/(NR/4)"\t"(NR/4)"\t"s}\'']) + output = self.run_command(command) + self.ofh.write("output:\n%s\n" % str(output)) + self.ofh.write("re.split('\\t', self.run_command(command)[0]:\n%s\n" % str(re.split('\\t', self.run_command(command)[0]))) + values = [] + for i in re.split('\\t', self.run_command(command)[0]): + if i == '': + values.append(float('nan')) + else: + values.append(float(i)) + self.ofh.write("values:\n%s\n" % str(values)) + self.ofh.write("values[0]:\n%s\n" % str(values[0])) + self.illumina_length_mean += values[0] + self.ofh.write("values[1]:\n%s\n" % str(values[1])) + self.illumina_read_count += int(values[1]) + self.ofh.write("values[2]:\n%s\n" % str(values[2])) + self.illumina_bases += int(values[2]) + # The original PIMA code inserts self.illumina_fastq into + # a list for no apparent reason. We don't do that here. + # self.illumina_length_mean /= len(self.illumina_fastq) + self.illumina_length_mean /= 1 + self.illumina_bases = self.format_kmg(self.illumina_bases, decimals=1) + + def start_doc(self): + #header_text = 'Analysis of %s' % self.analysis_name + self.doc = MdUtils(file_name=self.report_md, title='') + + def add_run_information(self): + self.ofh.write("\nXXXXXX In add_run_information\n\n") + self.doc.new_line() + self.doc.new_header(1, 'Run information') + # Tables in md.utils are implemented as a wrapping function. + Table_list = [ + "Category", + "Information", + "Date", + date.today(), + "ONT FAST5", + "N/A", + "ONT FASTQ", + "N/A", + "Illumina FASTQ", + self.wordwrap_markdown(self.analysis_name), + "Assembly", + self.wordwrap_markdown(self.assembly_name), + "Reference", + self.wordwrap_markdown(self.dbkey), + ] + self.doc.new_table(columns=2, rows=7, text=Table_list, text_align='left') + self.doc.new_line() + self.doc.new_line() + + def add_ont_library_information(self): + self.ofh.write("\nXXXXXX In add_ont_library_information\n\n") + if self.ont_n50 is None: + return + self.doc.new_line() + self.doc.new_header(2, 'ONT library statistics') + Table_List = [ + "Category", + "Quantity", + "ONT N50", + '{:,}'.format(self.ont_n50), + "ONT reads", + '{:,}'.format(self.ont_read_count), + "ONT bases", + '{:s}'.format(self.ont_bases), + "Illumina FASTQ", + self.wordwrap_markdown(self.illumina_fastq_file), + "Assembly", + self.wordwrap_markdown(self.assembly_name), + "Reference", + self.wordwrap_markdown(self.dbkey), + ] + self.doc.new_table(columns=2, rows=7, text=Table_List, text_align='left') + self.doc.new_line() + + def add_illumina_library_information(self): + self.ofh.write("\nXXXXXX In add_illumina_library_information\n\n") + if self.illumina_length_mean is None: + return + self.doc.new_line() + self.doc.new_header(2, 'Illumina library statistics') + Table_List = [ + "Illumina Info.", + "Quantity", + 'Illumina mean length', + '{:.1f}'.format(self.illumina_length_mean), + 'Illumina reads', + '{:,}'.format(self.illumina_read_count), + 'Illumina bases', + '{:s}'.format(self.illumina_bases) + ] + self.doc.new_table(columns=2, rows=4, text=Table_List, text_align='left') + + def add_assembly_information(self): + self.ofh.write("\nXXXXXX In add_assembly_information\n\n") + if self.assembly_fasta_file is None: + return + self.load_assembly() + self.doc.new_line() + self.doc.new_header(2, 'Assembly statistics') + Table_List = [ + "Category", + "Information", + "Contigs", + str(self.num_assembly_contigs), + "Assembly size", + str(self.assembly_size), + ] + self.doc.new_table(columns=2, rows=3, text=Table_List, text_align='left') + + def info_ont_fastq(self, fastq_file): + self.ofh.write("\nXXXXXX In info_ont_fastq, fastq_file:\n%s\n\n" % str(fastq_file)) + opener = 'cat' + if self.gzipped: + opener = 'gunzip -c' + else: + opener = 'cat' + command = ' '.join([opener, + fastq_file, + '| awk \'{getline;print length($0);s += length($1);getline;getline;}END{print "+"s}\'', + '| sort -gr', + '| awk \'BEGIN{bp = 0;f = 0}', + '{if(NR == 1){sub(/+/, "", $1);s=$1}else{bp += $1;if(bp > s / 2 && f == 0){n50 = $1;f = 1}}}', + 'END{printf "%d\\t%d\\t%d\\n", n50, (NR - 1), s;exit}\'']) + result = list(re.split('\\t', self.run_command(command)[0])) + if result[1] == '0': + self.error_out('No ONT reads found') + ont_n50, ont_read_count, ont_raw_bases = [int(i) for i in result] + + command = ' '.join([opener, + fastq_file, + '| awk \'{getline;print length($0);getline;getline;}\'']) + result = self.run_command(command) + result = list(filter(lambda x: x != '', result)) + ont_read_lengths = [int(i) for i in result] + + return ([ont_n50, ont_read_count, ont_raw_bases, ont_read_lengths]) + + def wordwrap_markdown(self, string): + if string: + if len(string) < 35: + return string + else: + if '/' in string: + adjust = string.split('/') + out = '' + max = 35 + for i in adjust: + out = out + '/' + i + if len(out) > max: + out += '<br>' + max += 35 + return out + else: + out = [string[i:i + 35] for i in range(0, len(string), 50)] + return '<br>'.join(out) + else: + return string + + def add_contig_info(self): + self.ofh.write("\nXXXXXX In add_contig_info\n\n") + if self.contig_info is None: + return + for method in ['ONT', 'Illumina']: + if method not in self.contig_info.index: + continue + self.doc.new_line() + self.doc.new_header(2, 'Assembly coverage by ' + method) + Table_List = ["Contig", "Length (bp)", "Coverage (X)"] + formatted = self.contig_info[method].copy() + formatted.iloc[:, 1] = formatted.iloc[:, 1].apply(lambda x: '{:,}'.format(x)) + for i in range(self.contig_info[method].shape[0]): + Table_List = Table_List + formatted.iloc[i, :].values.tolist() + row_count = int(len(Table_List) / 3) + self.doc.new_table(columns=3, rows=row_count, text=Table_List, text_align='left') + + def add_assembly_notes(self): + self.ofh.write("\nXXXXXX In add_assembly_notes\n\n") + if len(self.assembly_notes) == 0: + return + self.doc.new_line() + self.doc.new_line('<div style="page-break-after: always;"></div>') + self.doc.new_line() + self.doc.new_header(2, self.assembly_notes_title) + #for note in self.analysis.assembly_notes: + # self.doc.new_line(note) + + def add_contamination(self): + self.ofh.write("\nXXXXXX In add_contamination\n\n") + if len(self.kraken_fracs) == 0: + return + self.doc.new_line() + self.doc.new_header(2, 'Contamination check') + for read_type, kraken_fracs in self.kraken_fracs.iteritems(): + self.doc.new_line(read_type + ' classifications') + self.doc.new_line() + Table_List = ["Percent of Reads", "Reads", "Level", "Label"] + for index, row in kraken_fracs.iterrows(): + Table_List = Table_List + row.tolist() + row_count = int(len(Table_List) / 4) + self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left') + if self.contamination_methods_title not in self.methods: + self.methods[self.contamination_methods_title] = '' + method = 'Kraken2 was used to assign the raw reads into taxa.' + self.methods[self.contamination_methods_title] = self.methods[self.contamination_methods_title].append(pandas.Series(method)) + + def add_alignment(self): + self.ofh.write("\nXXXXXX In add_alignment\n\n") + # TODO: implement the draw_circos function for this. + if len(self.contig_alignment) > 0: + alignments = self.contig_alignment + else: + return + self.doc.new_line() + self.doc.new_header(level=2, title=self.alignment_title) + self.doc.new_line() + self.doc.new_header(level=3, title=self.snp_indel_title) + Table_1 = [ + "Category", + "Quantity", + 'SNPs', + #'{:,}'.format(self.analysis.quast_mismatches), + 'NA' + 'Small indels', + #'{:,}'.format(self.analysis.quast_indels) + 'NA' + ] + self.doc.new_table(columns=2, rows=3, text=Table_1, text_align='left') + self.doc.new_line('<div style="page-break-after: always;"></div>') + self.doc.new_line() + if len(self.alignment_notes) > 0: + self.doc.new_header(level=3, title=self.alignment_notes_title) + for note in self.alignment_notes: + self.doc.new_line(note) + for contig in alignments.index.tolist(): + contig_title = 'Alignment to %s' % contig + image_png = alignments[contig] + self.doc.new_line() + self.doc.new_header(level=3, title=contig_title) + self.doc.new_line(self.doc.new_inline_image(text='contig_title', path=os.path.abspath(image_png))) + self.doc.new_line('<div style="page-break-after: always;"></div>') + self.doc.new_line() + method = 'The genome assembly was aligned against the reference sequencing using dnadiff.' + self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method)) + + def add_features(self): + self.ofh.write("\nXXXXXX In add_features\n\n") + if len(self.feature_bed_files) == 0: + return + for bbf in self.feature_bed_files: + if os.path.getsize(bbf) > 0: + best = pandas.read_csv(filepath_or_buffer=bbf, sep='\t', header=None) + self.feature_hits[os.path.basename(bbf)] = best + if len(self.feature_hits) == 0: + return + self.ofh.write("self.feature_hits: %s\n" % str(self.feature_hits)) + self.doc.new_line() + self.doc.new_header(level=2, title=self.feature_title) + for feature_name in self.feature_hits.index.tolist(): + self.ofh.write("feature_name: %s\n" % str(feature_name)) + features = self.feature_hits[feature_name].copy() + self.ofh.write("features: %s\n" % str(features)) + if features.shape[0] == 0: + continue + features.iloc[:, 1] = features.iloc[:, 1].apply(lambda x: '{:,}'.format(x)) + features.iloc[:, 2] = features.iloc[:, 2].apply(lambda x: '{:,}'.format(x)) + self.doc.new_line() + self.doc.new_header(level=3, title=feature_name) + if (features.shape[0] == 0): + continue + for contig in pandas.unique(features.iloc[:, 0]): + self.ofh.write("contig: %s\n" % str(contig)) + self.doc.new_line(contig) + contig_features = features.loc[(features.iloc[:, 0] == contig), :] + self.ofh.write("contig_features: %s\n" % str(contig_features)) + Table_List = ['Start', 'Stop', 'Feature', 'Identity (%)', 'Strand'] + for i in range(contig_features.shape[0]): + self.ofh.write("i: %s\n" % str(i)) + feature = contig_features.iloc[i, :].copy(deep=True) + self.ofh.write("feature: %s\n" % str(feature)) + feature[4] = '{:.3f}'.format(feature[4]) + Table_List = Table_List + feature[1:].values.tolist() + self.ofh.write("Table_List: %s\n" % str(Table_List)) + row_count = int(len(Table_List) / 5) + self.ofh.write("row_count: %s\n" % str(row_count)) + self.doc.new_line() + self.doc.new_table(columns=7, rows=row_count, text=Table_List, text_align='left') + blastn_version = 'The genome assembly was queried for features using blastn.' + bedtools_version = 'Feature hits were clustered using bedtools and the highest scoring hit for each cluster was reported.' + method = '%s %s' % (blastn_version, bedtools_version) + self.methods[self.feature_methods_title] = self.methods[self.feature_methods_title].append(pandas.Series(method)) + + def add_feature_plots(self): + self.ofh.write("\nXXXXXX In add_feature_plots\n\n") + if len(self.feature_png_files) == 0: + return + self.doc.new_line() + self.doc.new_header(level=2, title='Feature Plots') + self.doc.new_paragraph('Only contigs with features are shown') + for feature_png_file in self.feature_png_files: + self.doc.new_line(self.doc.new_inline_image(text='Analysis', path=os.path.abspath(feature_png_file))) + + def add_mutations(self): + self.ofh.write("\nXXXXXX In add_mutations\n\n") + if len(self.mutation_regions_tsv_files) == 0: + return + try: + mutation_regions = pandas.read_csv(self.mutation_regions_bed_file, sep='\t', header=0, index_col=False) + except Exception: + # Likely an empty file. + return + amr_mutations = pandas.Series(dtype=object) + for region_i in range(mutation_regions.shape[0]): + region = mutation_regions.iloc[region_i, :] + region_name = str(region['name']) + self.ofh.write("Processing mutations for region %s\n" % region_name) + region_mutations_tsv_name = '%s_mutations.tsv' % region_name + if region_mutations_tsv_name not in self.mutation_regions_tsv_files: + continue + region_mutations_tsv = self.mutation_regions_tsv_files[region_mutations_tsv_name] + try: + region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False) + except Exception: + region_mutations = pandas.DataFrame() + if region_mutations.shape[0] == 0: + continue + # Figure out what kind of mutations are in this region + region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index) + region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel' + region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index) + region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index) + region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1) + region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']] + amr_mutations[region['name']] = region_mutations + # Report the mutations. + self.doc.new_line() + self.doc.new_header(level=2, title=self.mutation_title) + for region_name in amr_mutations.index.tolist(): + region_mutations = amr_mutations[region_name].copy() + self.doc.new_line() + self.doc.new_header(level=3, title=region_name) + if (region_mutations.shape[0] == 0): + self.doc.append('None') + continue + region_mutations.iloc[:, 1] = region_mutations.iloc[:, 1].apply(lambda x: '{:,}'.format(x)) + Table_List = ['Reference contig', 'Position', 'Reference', 'Alternate', 'Drug', 'Note'] + for i in range(region_mutations.shape[0]): + Table_List = Table_List + region_mutations.iloc[i, [0, 1, 3, 4, 5, 6]].values.tolist() + row_count = int(len(Table_List) / 6) + self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left') + method = '%s reads were mapped to the reference sequence using minimap2.' % self.read_type + self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method)) + method = 'Mutations were identified using samtools mpileup and varscan.' + self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method)) + + def add_amr_matrix(self): + self.ofh.write("\nXXXXXX In add_amr_matrix\n\n") + # Make sure that we have an AMR matrix to plot + #if not getattr(self.analysis, 'did_draw_amr_matrix', False): + # return + #amr_matrix_png = self.analysis.amr_matrix_png + #self.doc.new_line() + #self.doc.new_header(level=2, title=self.amr_matrix_title) + #self.doc.new_line('AMR genes and mutations with their corresponding drugs.') + #self.doc.new_line( + # self.doc.new_inline_image( + # text='AMR genes and mutations with their corresponding drugs', + # path=amr_matrix_png + # ) + #) + pass + + def add_large_indels(self): + self.ofh.write("\nXXXXXX In add_large_indels\n\n") + # Make sure we looked for mutations + #if len(self.analysis.large_indels) == 0: + # return + #large_indels = self.analysis.large_indels + #self.doc.new_line() + #self.doc.new_header(level=2, title=self.large_indel_title) + #for genome in ['Reference insertions', 'Query insertions']: + # genome_indels = large_indels[genome].copy() + # self.doc.new_line() + # self.doc.new_header(level=3, title=genome) + # if (genome_indels.shape[0] == 0): + # continue + # genome_indels.iloc[:, 1] = genome_indels.iloc[:, 1].apply(lambda x: '{:,}'.format(x)) + # genome_indels.iloc[:, 2] = genome_indels.iloc[:, 2].apply(lambda x: '{:,}'.format(x)) + # genome_indels.iloc[:, 3] = genome_indels.iloc[:, 3].apply(lambda x: '{:,}'.format(x)) + # Table_List = [ + # 'Reference contig', 'Start', 'Stop', 'Size (bp)' + # ] + # for i in range(genome_indels.shape[0]): + # Table_List = Table_List + genome_indels.iloc[i, :].values.tolist() + # row_count = int(len(Table_List) / 4) + # self.doc.new_table(columns=4, rows=row_count, text=Table_List, text_align='left') + #method = 'Large insertions or deletions were found as the complement of aligned regions using bedtools.' + #self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append( + # pandas.Series(method)) + #self.doc.new_line() + #self.doc.new_line('<div style="page-break-after: always;"></div>') + #self.doc.new_line() + pass + + def add_plasmids(self): + self.ofh.write("\nXXXXXX In add_plasmids\n\n") + """ + if not getattr(self.analysis, 'did_call_plasmids', False): + return + # Make sure we looked for mutations + #plasmids = self.analysis.plasmids + if plasmids is None: + return + plasmids = plasmids.copy() + self.doc.new_line() + #self.doc.new_header(level=2, title=self.analysis.plasmid_title) + if (plasmids.shape[0] == 0): + self.doc.new_line('None') + return + plasmids.iloc[:, 3] = plasmids.iloc[:, 3].apply(lambda x: '{:,}'.format(x)) + plasmids.iloc[:, 4] = plasmids.iloc[:, 4].apply(lambda x: '{:,}'.format(x)) + plasmids.iloc[:, 5] = plasmids.iloc[:, 5].apply(lambda x: '{:,}'.format(x)) + Table_List = [ + 'Genome contig', + 'Plasmid hit', + 'Plasmid acc.', + 'Contig size', + 'Aliged', + 'Plasmid size' + ] + for i in range(plasmids.shape[0]): + Table_List = Table_List + plasmids.iloc[i, 0:6].values.tolist() + row_count = int(len(Table_List) / 6) + self.doc.new_table(columns=6, rows=row_count, text=Table_List, text_align='left') + method = 'The plasmid reference database was queried against the genome assembly using minimap2.' + self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method)) + method = 'The resulting SAM was converted to a PSL using a custom version of sam2psl.' + self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method)) + method = 'Plasmid-to-genome hits were resolved using the pChunks algorithm.' + self.methods[self.plasmid_methods_title] = self.methods[self.plasmid_methods_title].append(pandas.Series(method)) + """ + pass + + def add_methods(self): + self.ofh.write("\nXXXXXX In add_methods\n\n") + self.doc.new_line('<div style="page-break-after: always;"></div>') + self.doc.new_line() + if len(self.methods) == 0: + return + self.doc.new_line() + self.doc.new_header(level=2, title=self.methods_title) + + for methods_section in self.methods.index.tolist(): + if self.methods[methods_section] is None or len(self.methods[methods_section]) == 0: + continue + self.doc.new_line() + self.doc.new_header(level=3, title=methods_section) + self.doc.new_paragraph(' '.join(self.methods[methods_section])) + + def add_summary(self): + self.ofh.write("\nXXXXXX In add_summary\n\n") + # Add summary title + self.doc.new_header(level=1, title=self.summary_title) + # First section of Summary + self.doc.new_header(level=1, title='CDC Advisory') + self.doc.new_paragraph(CDC_ADVISORY) + self.doc.new_line() + self.add_run_information() + self.add_ont_library_information() + methods = [] + if self.did_guppy_ont_fast5: + methods += ['ONT reads were basecalled using guppy'] + if self.did_qcat_ont_fastq: + methods += ['ONT reads were demultiplexed and trimmed using qcat'] + self.methods[self.basecalling_methods_title] = pandas.Series(methods) + self.add_illumina_library_information() + self.add_assembly_information() + self.add_contig_info() + self.add_assembly_notes() + if self.did_flye_ont_fastq: + method = 'ONT reads were assembled using Flye.' + self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.Series(method)) + if self.did_medaka_ont_assembly: + method = 'the genome assembly was polished using ont reads and medaka.' + self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.series(method)) + + def make_tex(self): + self.doc.new_table_of_contents(table_title='detailed run information', depth=2, marker="tableofcontents") + text = self.doc.file_data_text + text = text.replace("##--[", "") + text = text.replace("]--##", "") + self.doc.file_data_text = text + self.doc.create_md_file() + + def make_report(self): + self.ofh.write("\nXXXXXX In make_report\n\n") + self.start_doc() + self.add_summary() + self.add_contamination() + self.add_alignment() + self.add_features() + self.add_feature_plots() + self.add_mutations() + # TODO stuff working to here... + self.add_large_indels() + self.add_plasmids() + self.add_amr_matrix() + # self.add_snps() + self.add_methods() + self.make_tex() + # It took me quite a long time to find out that the value of the -t + # (implied) argument in the following command must be 'html' instead of + # the more logical 'pdf'. see the answer from snsn in this thread: + # https://github.com/jessicategner/pypandoc/issues/186 + self.ofh.write("\nXXXXX In make_report, calling pypandoc.convert_file...\n\n") + pypandoc.convert_file(self.report_md, + 'html', + extra_args=['--pdf-engine=weasyprint', '-V', '-css=%s' % self.pima_css], + outputfile='pima_report.pdf') + self.ofh.close() + + +parser = argparse.ArgumentParser() + +parser.add_argument('--analysis_name', action='store', dest='analysis_name', help='Sample identifier') +parser.add_argument('--assembly_fasta_file', action='store', dest='assembly_fasta_file', help='Assembly fasta file') +parser.add_argument('--assembly_name', action='store', dest='assembly_name', help='Assembly identifier') +parser.add_argument('--feature_bed_dir', action='store', dest='feature_bed_dir', help='Directory of best feature hits bed files') +parser.add_argument('--feature_png_dir', action='store', dest='feature_png_dir', help='Directory of best feature hits png files') +parser.add_argument('--contig_coverage_file', action='store', dest='contig_coverage_file', help='Contig coverage TSV file') +parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference genome') +parser.add_argument('--gzipped', action='store_true', dest='gzipped', default=False, help='Input sample is gzipped') +parser.add_argument('--illumina_fastq_file', action='store', dest='illumina_fastq_file', help='Input sample') +parser.add_argument('--mutation_regions_bed_file', action='store', dest='mutation_regions_bed_file', help='AMR mutation regions BRD file') +parser.add_argument('--mutation_regions_dir', action='store', dest='mutation_regions_dir', help='Directory of mutation regions TSV files') +parser.add_argument('--pima_css', action='store', dest='pima_css', help='PIMA css stypesheet') + +args = parser.parse_args() + +# Prepare the features BED files. +feature_bed_files = [] +for file_name in sorted(os.listdir(args.feature_bed_dir)): + file_path = os.path.abspath(os.path.join(args.feature_bed_dir, file_name)) + feature_bed_files.append(file_path) +# Prepare the features PNG files. +feature_png_files = [] +for file_name in sorted(os.listdir(args.feature_png_dir)): + file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name)) + feature_png_files.append(file_path) +# Prepare the mutation regions TSV files. +mutation_regions_files = [] +for file_name in sorted(os.listdir(args.mutation_regions_dir)): + file_path = os.path.abspath(os.path.join(args.feature_png_dir, file_name)) + mutation_regions_files.append(file_path) + +markdown_report = PimaReport(args.analysis_name, + args.assembly_fasta_file, + args.assembly_name, + feature_bed_files, + feature_png_files, + args.contig_coverage_file, + args.dbkey, + args.gzipped, + args.illumina_fastq_file, + args.mutation_regions_bed_file, + mutation_regions_files, + args.pima_css) +markdown_report.make_report()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pima_report.xml Fri Mar 03 22:06:23 2023 +0000 @@ -0,0 +1,80 @@ +<tool id="pima_report" name="PIMA: summary report" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="@PROFILE@"> + <description></description> + <macros> + <import>macros.xml</import> + </macros> + <expand macro="requirements"/> + <command detect_errors="exit_code"><![CDATA[ +#import re + +#set analysis_name = re.sub('[^\s\w\-]', '_', str($illumina_fastq_file.element_identifier)) +#set assembly_name = re.sub('[^\s\w\-]', '_', str($assembly_fasta_file.element_identifier)) + +mkdir feature_bed_dir && +mkdir feature_png_dir && +mkdir mutation_regions_dir && +touch 'pima_report.pdf' && + +#for $i in $features_bed: + #set file_name = $i.file_name + #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier)) + ln -s $i 'feature_bed_dir/$identifier' && +#end for +#for $i in $features_png: + #set file_name = $i.file_name + #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier)) + ln -s $i 'feature_png_dir/$identifier' && +#end for +#for $i in $mutation_regions: + #set file_name = $i.file_name + #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier)) + ln -s $i 'mutation_regions_dir/$identifier' && +#end for + +python '${__tool_directory__}/pima_report.py' +--analysis_name '$analysis_name' +--assembly_fasta_file '$assembly_fasta_file' +--assembly_name '$assembly_name' +--contig_coverage_file '$contig_coverage_file' +--dbkey '$aligned_sample.metadata.dbkey' +--feature_bed_dir 'feature_bed_dir' +--feature_png_dir 'feature_png_dir' +#if $illumina_fastq_file.ext.endswith(".gz"): + --gzipped +#end if +--illumina_fastq_file '$illumina_fastq_file' +--mutation_regions_dir 'mutation_regions_dir' +--mutation_regions_bed_file '$mutation_regions_bed_file' +--pima_css '${__tool_directory__}/pima.css' +&& mv 'pima_report.pdf' '$output' + ]]></command> + <inputs> + <param name="aligned_sample" type="data" format="bam" label="Aligned sample BAM file"/> + <param name="assembly_fasta_file" type="data" format="fasta" label="Assembly FASTA file"/> + <param name="features_bed" format="bed" type="data_collection" collection_type="list" label="Collection of best feature hits BED files"/> + <param name="features_png" format="png" type="data_collection" collection_type="list" label="Collection of best feature hits PNG files"/> + <param name="contig_coverage_file" type="data" format="tabular,tsv" label="Contig coverage tabular file"/> + <param name="illumina_fastq_file" type="data" format="fastqsanger,fastqsanger.gz" label="Fastq sample file"/> + <param name="mutation_regions" format="tabular,tsv" type="data_collection" collection_type="list" label="Collection of mutation regions tabular files"/> + <param name="mutation_regions_bed_file" type="data" format="mutations_regions,bed" label="Mutation regions BED file"/> + </inputs> + <outputs> + <data name="output" format="pdf"/> + </outputs> + <tests> + <test> + <param name="aligned_sample" value="aligned_sample.bam" ftype="bam"/> + <param name="assembly_fasta_file" value="assembly_fasta.fasta" ftype="fasta"/> + <param name="contig_coverage_file" value="contig_coverage.tabular" ftype="tabular"/> + <param name="illumina_fastq_file" value="illumina_fastq.fastq" ftype="fastq"/> + <output name="output" value="output.pdf" ftype="pdf"/> + </test> + </tests> + <help> +**What it does** + +Generates the PIMA analysis summary report. + </help> + <expand macro="citations"/> +</tool> +