Mercurial > repos > greg > pima_report

changeset 13:f03c80bb22e9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 16 Mar 2023 14:42:13 +0000
parents 99613333fd1f
children 95b1d1a9497d
files pima_report.py pima_report.xml
diffstat 2 files changed, 66 insertions(+), 50 deletions(-) [+]
line wrap: on
line diff
--- a/pima_report.py	Fri Mar 10 16:35:16 2023 +0000
+++ b/pima_report.py	Thu Mar 16 14:42:13 2023 +0000
@@ -16,13 +16,13 @@
 class PimaReport:
 
     def __init__(self, analysis_name=None, amr_deletions_file=None, amr_matrix_files=None, assembly_fasta_file=None,
-                 assembly_name=None, bedtools_version=None, blastn_version=None, compute_sequence_length_file=None,
-                 contig_coverage_file=None, dbkey=None, dnadiff_snps_file=None, dnadiff_version=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, kraken2_report_file=None,
-                 kraken2_version=None, minimap2_version=None, mutation_regions_bed_file=None,
-                 mutation_regions_tsv_files=None, pima_css=None, plasmids_file=None, reference_insertions_file=None,
-                 samtools_version=None, varscan_version=None):
+                 assembly_name=None, bedtools_version=None, blastn_version=None, circos_files=None,
+                 compute_sequence_length_file=None, contig_coverage_file=None, dbkey=None, dnadiff_snps_file=None,
+                 dnadiff_version=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,
+                 kraken2_report_file=None, kraken2_version=None, minimap2_version=None, mutation_regions_bed_file=None,
+                 mutation_regions_tsv_files=None, pima_css=None, plasmids_file=None, quast_report_file=None,
+                 reference_insertions_file=None, samtools_version=None, varscan_version=None):
         self.ofh = open("process_log.txt", "w")
 
         self.ofh.write("amr_deletions_file: %s\n" % str(amr_deletions_file))
@@ -32,6 +32,7 @@
         self.ofh.write("assembly_name: %s\n" % str(assembly_name))
         self.ofh.write("bedtools_version: %s\n" % str(bedtools_version))
         self.ofh.write("blastn_version: %s\n" % str(blastn_version))
+        self.ofh.write("circos_files: %s\n" % str(circos_files))
         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))
@@ -51,6 +52,7 @@
         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("quast_report_file: %s\n" % str(quast_report_file))
         self.ofh.write("reference_insertions_file: %s\n" % str(reference_insertions_file))
         self.ofh.write("samtools_version: %s\n" % str(samtools_version))
         self.ofh.write("varscan_version: %s\n" % str(varscan_version))
@@ -62,7 +64,8 @@
         # Inputs
         self.amr_deletions_file = amr_deletions_file
         self.amr_matrix_files = amr_matrix_files
-        self.analysis_name = re.sub('_', '.', analysis_name.rstrip(' _consensus_'))
+        self.analysis_name = analysis_name.split('_')[0]
+        self.ofh.write("self.analysis_name: %s\n" % str(self.analysis_name))
         self.assembly_fasta_file = assembly_fasta_file
         self.assembly_name = re.sub('_', '.', assembly_name.rstrip(' _consensus_'))
         if bedtools_version is None:
@@ -73,6 +76,7 @@
             self.blastn_version = 'blastn (version unknown)'
         else:
             self.blastn_version = re.sub('_', '.', blastn_version.rstrip(' _features_'))
+        self.circos_files = circos_files
         self.compute_sequence_length_file = compute_sequence_length_file
         self.contig_coverage_file = contig_coverage_file
         self.dbkey = dbkey
@@ -108,6 +112,8 @@
         self.ont_read_count = None
         self.pima_css = pima_css
         self.plasmids_file = plasmids_file
+        self.quast_report_file = quast_report_file
+        self.reference_insertions_file = reference_insertions_file
         self.reference_insertions_file = reference_insertions_file
         if samtools_version is None:
             self.samtools_version = 'samtools (version unknown)'
@@ -139,7 +145,7 @@
         self.plasmid_title = 'Plasmid annotation'
         self.reference_methods_title = 'Reference comparison'
         self.snp_indel_title = 'SNPs and small indels'
-        self.summary_title = 'Analysis of %s' % analysis_name
+        self.summary_title = 'Summary'
 
         # Methods
         self.methods = pandas.Series(dtype='float64')
@@ -163,13 +169,9 @@
         self.illumina_length_mean = 0
         self.illumina_read_count = 0
         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
@@ -265,7 +267,8 @@
         self.illumina_bases = self.format_kmg(self.illumina_bases, decimals=1)
 
     def start_doc(self):
-        self.doc = MdUtils(file_name=self.report_md, title='')
+        header_text = 'Analysis of ' + self.analysis_name
+        self.doc = MdUtils(file_name=self.report_md, title=header_text)
 
     def add_run_information(self):
         self.ofh.write("\nXXXXXX In add_run_information\n\n")
@@ -319,7 +322,7 @@
 
     def add_illumina_library_information(self):
         self.ofh.write("\nXXXXXX In add_illumina_library_information\n\n")
-        if self.illumina_length_mean is None:
+        if self.illumina_length_mean == 0:
             return
         self.doc.new_line()
         self.doc.new_header(2, 'Illumina library statistics')
@@ -475,38 +478,41 @@
 
     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.quast_mismatches),
-            'Small indels',
-            '{:,}'.format(self.quast_indels)
-        ]
+        if self.quast_report_file is not None:
+            # Process quast values.
+            quast_report = pandas.read_csv(self.quast_report_file, header=0, index_col=0, sep='\t')
+            quast_mismatches = int(float(quast_report.loc['# mismatches per 100 kbp', :][0]) * (float(quast_report.loc['Total length (>= 0 bp)', :][0]) / 100000.))
+            quast_indels = int(float(quast_report.loc['# indels per 100 kbp', :][0]) * (float(quast_report.loc['Total length (>= 0 bp)', :][0]) / 100000.))
+            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(quast_mismatches),
+                'Small indels',
+                '{:,}'.format(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>')
         self.doc.new_line()
+        # TODO: self.alignment_notes is not currently populated.
         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()
+        if len(self.circos_files) > 0:
+            # Add circos PNG files.
+            for circos_file in self.circos_files:
+                contig = os.path.basename(circos_file)
+                contig_title = 'Alignment to %s' % 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(circos_file)))
+                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 version %s.' % self.dnadiff_version
         self.methods[self.reference_methods_title] = self.methods[self.reference_methods_title].append(pandas.Series(method))
 
@@ -643,7 +649,7 @@
                 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 %s.' % (self.read_type, self.minimap2_version)
         self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
-        method = 'Mutations were identified using %s mpileup and %s.' % (self.samtools_version, self.varscan_version)
+        method = 'Mutations were identified using %s and %s.' % (self.samtools_version, self.varscan_version)
         self.methods[self.mutation_methods_title] = self.methods[self.mutation_methods_title].append(pandas.Series(method))
 
     def add_amr_matrix(self):
@@ -672,15 +678,6 @@
             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)
@@ -835,6 +832,7 @@
 parser.add_argument('--assembly_name', action='store', dest='assembly_name', help='Assembly identifier')
 parser.add_argument('--bedtools_version', action='store', dest='bedtools_version', default=None, help='Bedtools version string')
 parser.add_argument('--blastn_version', action='store', dest='blastn_version', default=None, help='Blastn version string')
+parser.add_argument('--circos_png_dir', action='store', dest='circos_png_dir', help='Directory of circos PNG files')
 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')
@@ -854,6 +852,7 @@
 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('--quast_report_file', action='store', dest='quast_report_file', help='Quast report tabular file')
 parser.add_argument('--reference_insertions_file', action='store', dest='reference_insertions_file', help='Reference insertions BED file')
 parser.add_argument('--samtools_version', action='store', dest='samtools_version', default=None, help='Samtools version string')
 parser.add_argument('--varscan_version', action='store', dest='varscan_version', default=None, help='Varscan version string')
@@ -865,6 +864,11 @@
 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 circos PNG files.
+circos_files = []
+for file_name in sorted(os.listdir(args.circos_png_dir)):
+    file_path = os.path.abspath(os.path.join(args.circos_png_dir, file_name))
+    circos_files.append(file_path)
 # Prepare the features BED files.
 feature_bed_files = []
 for file_name in sorted(os.listdir(args.feature_bed_dir)):
@@ -888,6 +892,7 @@
                              args.assembly_name,
                              args.bedtools_version,
                              args.blastn_version,
+                             circos_files,
                              args.compute_sequence_length_file,
                              args.contig_coverage_file,
                              args.dbkey,
@@ -907,6 +912,7 @@
                              mutation_regions_files,
                              args.pima_css,
                              args.plasmids_file,
+                             args.quast_report_file,
                              args.reference_insertions_file,
                              args.samtools_version,
                              args.varscan_version)
--- a/pima_report.xml	Fri Mar 10 16:35:16 2023 +0000
+++ b/pima_report.xml	Thu Mar 16 14:42:13 2023 +0000
@@ -36,6 +36,7 @@
 #end if
 
 mkdir amr_matrix_png_dir &&
+mkdir circos_png_dir &&
 mkdir feature_bed_dir &&
 mkdir feature_png_dir &&
 mkdir mutation_regions_dir &&
@@ -46,6 +47,11 @@
     #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier))
     ln -s $i 'amr_matrix_png_dir/$identifier' &&
 #end for
+#for $i in $circos_png:
+    #set file_name = $i.file_name
+    #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier))
+    ln -s $i 'circos_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))
@@ -74,6 +80,7 @@
 #if str($blastn_features) not in ['None', '']:
     --blastn_version '$blastn_version'
 #end if
+--circos_png_dir 'circos_png_dir'
 --compute_sequence_length_file '$compute_sequence_length_file'
 --contig_coverage_file '$contig_coverage_file'
 --dbkey '$aligned_sample.metadata.dbkey'
@@ -103,6 +110,7 @@
 --mutation_regions_bed_file '$mutation_regions_bed_file'
 --pima_css '${__tool_directory__}/pima.css'
 --plasmids_file '$plasmids_file'
+--quast_report_file '$quast_report_file'
 --reference_insertions_file '$reference_insertions_file'
 #if str($samtools_pileup_file) not in ['None', '']:
     --samtools_version '$samtools_version'
@@ -119,6 +127,7 @@
         <param name="assembly_fasta_file" type="data" format="fasta" label="Assembly FASTA file"/>
         <param name="bedtools_complementbed_file" type="data" format="bed" label="Bedtools ComplementBed BED file"/>
         <param name="blastn_features" format="tabular" type="data_collection" collection_type="list" label="Collection of blastn tabular files"/>
+        <param name="circos_png" format="png" type="data_collection" collection_type="list" label="Collection of circos PNG files"/>
         <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"/>
@@ -131,6 +140,7 @@
         <param name="minimap2_bam_file" type="data" format="bam" label="Minimap2 BAM 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="quast_report_file" type="data" format="tabular" label="Quast report tabular 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"/>
         <param name="samtools_pileup_file" type="data" format="pileup" label="Samtools pileup file"/>