changeset 1:67d0939b56b0 draft

Uploaded
author greg
date Tue, 07 Mar 2023 16:05:14 +0000
parents 0a558f444c98
children 9cb62054a87a
files pima_report.py pima_report.xml
diffstat 2 files changed, 248 insertions(+), 97 deletions(-) [+]
line wrap: on
line diff
--- a/pima_report.py	Fri Mar 03 22:06:23 2023 +0000
+++ b/pima_report.py	Tue Mar 07 16:05:14 2023 +0000
@@ -15,37 +15,57 @@
 
 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):
+    def __init__(self, analysis_name=None, amr_deletions_file=None, amr_matrix_files=None, assembly_fasta_file=None,
+                 assembly_name=None, compute_sequence_length_file=None, contig_coverage_file=None, dbkey=None,
+                 dnadiff_snps_file=None, feature_bed_files=None, feature_png_files=None, flye_assembly_info_file=None,
+                 flye_version=None, genome_insertions_file=None, gzipped=None, illumina_fastq_file=None,
+                 mutation_regions_bed_file=None, mutation_regions_tsv_files=None, pima_css=None, plasmids_file=None,
+                 reference_insertions_file=None):
         self.ofh = open("process_log.txt", "w")
 
+        self.ofh.write("amr_deletions_file: %s\n" % str(amr_deletions_file))
+        self.ofh.write("amr_matrix_files: %s\n" % str(amr_matrix_files))
         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("compute_sequence_length_file: %s\n" % str(compute_sequence_length_file))
+        self.ofh.write("contig_coverage_file: %s\n" % str(contig_coverage_file))
+        self.ofh.write("dbkey: %s\n" % str(dbkey))
+        self.ofh.write("dnadiff_snps_file: %s\n" % str(dnadiff_snps_file))
         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("flye_assembly_info_file: %s\n" % str(flye_assembly_info_file))
+        self.ofh.write("flye_version: %s\n" % str(flye_version))
         self.ofh.write("gzipped: %s\n" % str(gzipped))
+        self.ofh.write("genome_insertions_file: %s\n" % str(genome_insertions_file))
         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))
+        self.ofh.write("plasmids_file: %s\n" % str(plasmids_file))
+        # self.ofh.write("reference_genome: %s\n" % str(reference_genome))
+        self.ofh.write("reference_insertions_file: %s\n" % str(reference_insertions_file))
 
         # General
         self.doc = None
         self.report_md = 'pima_report.md'
 
         # Inputs
+        self.amr_deletions_file = amr_deletions_file
+        self.amr_matrix_files = amr_matrix_files
         self.analysis_name = analysis_name
         self.assembly_fasta_file = assembly_fasta_file
         self.assembly_name = assembly_name
+        self.compute_sequence_length_file = compute_sequence_length_file
+        self.contig_coverage_file = contig_coverage_file
+        self.dbkey = dbkey
+        self.dnadiff_snps_file = dnadiff_snps_file
         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.flye_assembly_info_file = flye_assembly_info_file
+        self.flye_version = flye_version
         self.gzipped = gzipped
+        self.genome_insertions_file = genome_insertions_file
         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
@@ -54,6 +74,9 @@
         self.ont_n50 = None
         self.ont_read_count = None
         self.pima_css = pima_css
+        self.plasmids_file = plasmids_file
+        # self.reference_genome = reference_genome
+        self.reference_insertions_file = reference_insertions_file
 
         # Titles
         self.alignment_title = 'Comparison with reference'
@@ -73,9 +96,9 @@
         self.mutation_title = 'Mutations found in the sample'
         self.mutation_methods_title = 'Mutation screening'
         self.plasmid_methods_title = 'Plasmid annotation'
+        self.plasmid_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
@@ -98,7 +121,6 @@
         # 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
@@ -106,6 +128,11 @@
         self.illumina_bases = 0
         self.mean_coverage = 0
         self.num_assembly_contigs = 0
+        # TODO: should the following 2 values  be passed as  parameters?
+        self.ont_n50_min = 2500
+        self.ont_coverage_min = 30
+        self.quast_indels = 0
+        self.quast_mismatches = 0
 
         # Actions
         self.did_guppy_ont_fast5 = False
@@ -139,6 +166,22 @@
         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()
+        if self.mean_coverage <= self.ont_coverage_min:
+            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
+            self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+        # Report if some contigs have low coverage.
+        low_coverage = self.contig_info[self.read_type].loc[self.contig_info[self.read_type]['coverage'] < self.ont_coverage_min, :]
+        if low_coverage.shape[0] >= 0:
+            for contig_i in range(low_coverage.shape[0]):
+                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
+                self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+        # See if some contigs have anolously low coverage.
+        fold_coverage = self.contig_info[self.read_type]['coverage'] / self.mean_coverage
+        low_coverage = self.contig_info[self.read_type].loc[fold_coverage < 1 / 5, :]
+        if low_coverage.shape[0] >= 0 :
+            for contig_i in range(low_coverage.shape[0]):
+                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
+                self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
 
     def load_fasta(self, fasta):
         sequence = pandas.Series(dtype=object)
@@ -185,7 +228,6 @@
         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):
@@ -256,6 +298,20 @@
         ]
         self.doc.new_table(columns=2, rows=4, text=Table_List, text_align='left')
 
+    def evaluate_assembly(self) :
+        assembly_info = pandas.read_csv(self.compute_sequence_length_file, sep='\t', header=None)
+        assembly_info.columns = ['contig', 'length']
+        self.contig_sizes = assembly_info
+        # Take a look at the number of contigs, their sizes,
+        # and circularity.  Warn if things don't look good.
+        if assembly_info.shape[0] > 4:
+            warning = 'Assembly produced {:d} contigs, more than ususally expected; assembly may be fragmented'.format(assembly_info.shape[0])
+            self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+        small_contigs = assembly_info.loc[assembly_info['length'] <= 3000, :]
+        if small_contigs.shape[0] > 0:
+            warning = 'Assembly produced {:d} small contigs ({:s}); assembly may include spurious sequences.'.format(small_contigs.shape[0], ', '.join(small_contigs['contig']))
+            self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
+
     def add_assembly_information(self):
         self.ofh.write("\nXXXXXX In add_assembly_information\n\n")
         if self.assembly_fasta_file is None:
@@ -291,15 +347,17 @@
         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])
+        # TODO: the following are not yet used...
+        # ont_read_lengths = [int(i) for i in result]
+        # ont_bases = self.format_kmg(ont_raw_bases, decimals=1)
+        if ont_n50 <= self.ont_n50_min:
+            warning = 'ONT N50 (%s) is less than the recommended minimum (%s)' % (str(ont_n50), str(self.ont_n50_min))
+            self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
 
     def wordwrap_markdown(self, string):
         if string:
@@ -347,8 +405,8 @@
         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)
+        for note in self.assembly_notes:
+            self.doc.new_line(note)
 
     def add_contamination(self):
         self.ofh.write("\nXXXXXX In add_contamination\n\n")
@@ -384,11 +442,9 @@
             "Category",
             "Quantity",
             'SNPs',
-            #'{:,}'.format(self.analysis.quast_mismatches),
-            'NA'
+            '{:,}'.format(self.quast_mismatches),
             'Small indels',
-            #'{:,}'.format(self.analysis.quast_indels)
-            'NA'
+            '{:,}'.format(self.quast_indels)
         ]
         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>')
@@ -444,12 +500,17 @@
                     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()
+                    # FIXME: Uncommenting the following line (which is
+                    # https://github.com/appliedbinf/pima_md/blob/main/MarkdownReport.py#L317)
+                    # will cause the job to fail with this exception:
+                    # ValueError: columns * rows is not equal to text length
+                    # 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')
+                self.ofh.write("Before new_table, len(Table_List):: %s\n" % str(len(Table_List)))
+                self.doc.new_table(columns=5, 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)
@@ -469,11 +530,34 @@
         self.ofh.write("\nXXXXXX In add_mutations\n\n")
         if len(self.mutation_regions_tsv_files) == 0:
             return
-        try:
+        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
+        # TODO: this is the only place where reference_genome is used,
+        # so I'm commenting it out for now.  We need to confirm if these
+        # errors that require the reference genmoe being passed are necessary.
+        # If so, we'll need to implement data tables in this tool.
+        # Make sure that the positions in the BED file fall within
+        # the chromosomes provided in the reference sequence.
+        """
+        for mutation_region in range(mutation_regions.shape[0]):
+            mutation_region = mutation_regions.iloc[mutation_region, :]
+            if not (mutation_region[0] in self.reference_genome):
+                self.ofh.write("\nMutation region: %s not found in reference genome.\n" % ' '.join(mutation_region.astype(str)))
+                continue
+            if not isinstance(mutation_region[1], int):
+                self.ofh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(mutation_region[1]))
+                break
+            elif not isinstance(mutation_region[2], int):
+                self.ofh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(mutation_region[2]))
+                break
+            if mutation_region[1] <= 0 or mutation_region[2] <= 0:
+                self.ofh.write("\nMutation region %s starts before the reference sequence.\n" % ' '.join(mutation_region.astype(str)))
+            if mutation_region[1] > len(self.reference_genome[mutation_region[0]].seq) or mutation_region[2] > len(self.reference_genome[mutation_region[0]].seq):
+                self.ofh.write("\nMutation region %s ends after the reference sequence.\n" % ' '.join(mutation_region.astype(str)))
+        """
         amr_mutations = pandas.Series(dtype=object)
         for region_i in range(mutation_regions.shape[0]):
             region = mutation_regions.iloc[region_i, :]
@@ -483,13 +567,13 @@
             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:
+            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
+            # 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)
@@ -521,78 +605,85 @@
     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
+        if len(self.amr_matrix_files) == 0:
+            return
+        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')
+        for amr_matrix_file in self.amr_matrix_files:
+            self.doc.new_line(self.doc.new_inline_image(text='AMR genes and mutations with their corresponding drugs',
+                                                        path=os.path.abspath(amr_matrix_file)))
 
     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
+        large_indels = pandas.Series(dtype='float64')
+        # Pull in insertions.
+        try:
+            reference_insertions = pandas.read_csv(filepath_or_buffer=self.reference_insertions_file, sep='\t', header=None)
+        except Exception:
+            reference_insertions = pandas.DataFrame()
+        try:
+            genome_insertions = pandas.read_csv(filepath_or_buffer=self.genome_insertions_file, sep='\t', header=None)
+        except Exception:
+            genome_insertions = pandas.DataFrame()
+        large_indels['Reference insertions'] = reference_insertions
+        large_indels['Query insertions'] = genome_insertions
+        # TODO: we don't seem to be reporting snps and deletions for some reason...
+        # Pull in the number of SNPs and small indels.
+        try:
+            snps = pandas.read_csv(filepath_or_buffer=self.dnadiff_snps_file, sep='\t', header=None)
+            # TODO: the following is not used...
+            # small_indels = snps.loc[(snps.iloc[:, 1] == '.') | (snps.iloc[:, 2] == '.'), :]
+            snps = snps.loc[(snps.iloc[:, 1] != '.') & (snps.iloc[:, 2] != '.'), :]
+        except Exception:
+            snps = pandas.DataFrame()
+        # Pull in deletions.
+        try:
+            amr_deletions = pandas.read_csv(filepath_or_buffer=self.amr_deletion_file, sep='\t', header=None)
+        except Exception:
+            amr_deletions = pandas.DataFrame()
+        if amr_deletions.shape[0] > 0:
+            amr_deletions.columns = ['contig', 'start', 'stop', 'name', 'type', 'drug', 'note']
+            amr_deletions = amr_deletions.loc[amr_deletions['type'].isin(['large-deletion', 'any']), :]
+        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()
 
     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:
+        try :
+            plasmids = pandas.read_csv(filepath_or_buffer=self.plasmids_file, sep='\t', header=0)
+        except Exception:
             return
         plasmids = plasmids.copy()
         self.doc.new_line()
-        #self.doc.new_header(level=2, title=self.analysis.plasmid_title)
+        self.doc.new_header(level=2, title=self.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'
-        ]
+        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)
@@ -603,8 +694,6 @@
         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")
@@ -614,7 +703,6 @@
             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
@@ -639,15 +727,25 @@
             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_contig_info()
+        self.evaluate_assembly()
         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.'
+        if self.flye_assembly_info_file is not None:
+            method = 'ONT reads were assembled using %s' % self.flye_version
             self.methods[self.assembly_methods_title] = self.methods[self.assembly_methods_title].append(pandas.Series(method))
+            # Pull in the assembly summary and look at the coverage.
+            assembly_info = pandas.read_csv(self.flye_assembly_info_file, header=0, index_col=0, sep='\t')
+            # Look for non-circular contigs.
+            open_contigs = assembly_info.loc[assembly_info['circ.'] == 'N', :]
+            if open_contigs.shape[0] > 0:
+                open_contig_ids = open_contigs.index.values
+                warning = 'Flye reported {:d} open contigs ({:s}); assembly may be incomplete.'.format(open_contigs.shape[0], ', '.join(open_contig_ids))
+                self.assembly_notes = self.assembly_notes.append(pandas.Series(warning))
         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))
+        self.info_ont_fastq(self.illumina_fastq_file)
+        self.add_assembly_notes()
 
     def make_tex(self):
         self.doc.new_table_of_contents(table_title='detailed run information', depth=2, marker="tableofcontents")
@@ -666,9 +764,9 @@
         self.add_features()
         self.add_feature_plots()
         self.add_mutations()
-        # TODO stuff working to here...
         self.add_large_indels()
         self.add_plasmids()
+        # TODO stuff working to here...
         self.add_amr_matrix()
         # self.add_snps()
         self.add_methods()
@@ -687,21 +785,36 @@
 
 parser = argparse.ArgumentParser()
 
+parser.add_argument('--amr_deletions_file', action='store', dest='amr_deletions_file', help='AMR deletions BED file')
+parser.add_argument('--amr_matrix_png_dir', action='store', dest='amr_matrix_png_dir', help='Directory of AMR matrix PNG files')
 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('--compute_sequence_length_file', action='store', dest='compute_sequence_length_file', help='Comnpute sequence length tabular file')
+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 identifier')
+parser.add_argument('--dnadiff_snps_file', action='store', dest='dnadiff_snps_file', help='DNAdiff snps tabular file')
 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('--flye_assembly_info_file', action='store', dest='flye_assembly_info_file', default=None, help='Flye assembly info tabular file')
+parser.add_argument('--flye_version', action='store', dest='flye_version', default=None, help='Flye version string')
+parser.add_argument('--genome_insertions_file', action='store', dest='genome_insertions_file', help='Genome insertions BED file')
 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')
+parser.add_argument('--plasmids_file', action='store', dest='plasmids_file', help='pChunks plasmids TSV file')
+parser.add_argument('--reference_insertions_file', action='store', dest='reference_insertions_file', help='Reference insertions BED file')
+# parser.add_argument('--reference_genome', action='store', dest='reference_genome', help='Reference genome fasta file')
 
 args = parser.parse_args()
 
+# Prepare the AMR matrix PNG files.
+amr_matrix_files = []
+for file_name in sorted(os.listdir(args.amr_matrix_png_dir)):
+    file_path = os.path.abspath(os.path.join(args.amr_matrix_png_dir, file_name))
+    amr_matrix_files.append(file_path)
 # Prepare the features BED files.
 feature_bed_files = []
 for file_name in sorted(os.listdir(args.feature_bed_dir)):
@@ -719,15 +832,24 @@
     mutation_regions_files.append(file_path)
 
 markdown_report = PimaReport(args.analysis_name,
+                             args.amr_deletions_file,
+                             amr_matrix_files,
                              args.assembly_fasta_file,
                              args.assembly_name,
+                             args.compute_sequence_length_file,
+                             args.contig_coverage_file,
+                             args.dbkey,
+                             args.dnadiff_snps_file,
                              feature_bed_files,
                              feature_png_files,
-                             args.contig_coverage_file,
-                             args.dbkey,
+                             args.flye_assembly_info_file,
+                             args.flye_version,
+                             args.genome_insertions_file,
                              args.gzipped,
                              args.illumina_fastq_file,
                              args.mutation_regions_bed_file,
                              mutation_regions_files,
-                             args.pima_css)
+                             args.pima_css,
+                             args.plasmids_file,
+                             args.reference_insertions_file)
 markdown_report.make_report()
--- a/pima_report.xml	Fri Mar 03 22:06:23 2023 +0000
+++ b/pima_report.xml	Tue Mar 07 16:05:14 2023 +0000
@@ -9,12 +9,22 @@
 
 #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))
+#if str($flye_assembly_info_file) not in ['None', '']:
+    #set flye_version = re.sub('[^\s\w\-]', '_', str($flye_assembly_info_file.element_identifier))
+    #set flye_version = $flye_version.rstrip('(assembly info)')
+#end if
 
+mkdir amr_matrix_png_dir &&
 mkdir feature_bed_dir &&
 mkdir feature_png_dir &&
 mkdir mutation_regions_dir &&
 touch 'pima_report.pdf' &&
 
+#for $i in $amr_matrices_png:
+    #set file_name = $i.file_name
+    #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier))
+    ln -s $i 'amr_matrix_png_dir/$identifier' &&
+#end for
 #for $i in $features_bed:
     #set file_name = $i.file_name
     #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier))
@@ -32,13 +42,22 @@
 #end for
 
 python '${__tool_directory__}/pima_report.py' 
+--amr_matrix_png_dir 'amr_matrix_png_dir'
+--amr_deletions_file '$amr_deletions_file'
 --analysis_name '$analysis_name'
 --assembly_fasta_file '$assembly_fasta_file'
 --assembly_name '$assembly_name'
+--compute_sequence_length_file '$compute_sequence_length_file'
 --contig_coverage_file '$contig_coverage_file'
 --dbkey '$aligned_sample.metadata.dbkey'
+--dnadiff_snps_file '$dnadiff_snps_file'
 --feature_bed_dir 'feature_bed_dir'
 --feature_png_dir 'feature_png_dir'
+#if str($flye_assembly_info_file) not in ['None', '']:
+    --flye_assembly_info_file '$flye_assembly_info_file'
+    --flye_version '$flye_version'
+#end if
+--genome_insertions_file '$genome_insertions_file'
 #if $illumina_fastq_file.ext.endswith(".gz"):
     --gzipped
 #end if
@@ -46,17 +65,27 @@
 --mutation_regions_dir 'mutation_regions_dir'
 --mutation_regions_bed_file '$mutation_regions_bed_file'
 --pima_css '${__tool_directory__}/pima.css'
+--plasmids_file '$plasmids_file'
+--reference_insertions_file '$reference_insertions_file'
 && mv 'pima_report.pdf' '$output'
     ]]></command>
     <inputs>
+        <param name="amr_matrices_png" format="png" type="data_collection" collection_type="list" label="Collection of AMR matrix PNG files"/>
         <param name="aligned_sample" type="data" format="bam" label="Aligned sample BAM file"/>
+        <param name="amr_deletions_file" type="data" format="bed" label="AMR deletions BED file"/>
         <param name="assembly_fasta_file" type="data" format="fasta" label="Assembly FASTA file"/>
+        <param name="compute_sequence_length_file" type="data" format="tabular,tsv" label="Compute sequence length tabular file"/>
+        <param name="contig_coverage_file" type="data" format="tabular,tsv" label="Contig coverage tabular file"/>
+        <param name="dnadiff_snps_file" type="data" format="tabular" label="DNAdiff snps tabular 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="flye_assembly_info_file" type="data" format="tabular,tsv" optional="true" label="Flye assembly info tabular file" help="Optional, ignored if not selected"/>
+        <param name="genome_insertions_file" type="data" format="bed" label="Genome insertions BED 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"/>
+        <param name="reference_insertions_file" type="data" format="bed" label="Reference insertions BED file"/>
+        <param name="plasmids_file" type="data" format="tsv" label="pChunks plasmids TSV file"/>
     </inputs>
     <outputs>
         <data name="output" format="pdf"/>