<commit_before>#!/usr/bin/env python

""" MultiQC module to parse output from Bismark """

from __future__ import print_function
from collections import defaultdict, OrderedDict
import io
import json
import logging
import os
import re

import multiqc
from multiqc import config

# Initialise the logger
log = logging.getLogger('MultiQC : {0:<14}'.format('Bismark'))

# Log parsing regexes
regexes = {
    'alignment': {
        'total_reads': r"^Sequence(?:s| pairs) analysed in total:\s+(\d+)$",
        'aligned_reads': r"^Number of(?: paired-end)? alignments with a unique best hit(?: from the different alignments)?:\s+(\d+)$",
        'no_alignments': r"^Sequence(?:s| pairs) with no alignments under any condition:\s+(\d+)$",
        'ambig_reads': r"^Sequence(?:s| pairs) did not map uniquely:\s+(\d+)$",
        'discarded_reads': r"^Sequence(?:s| pairs) which were discarded because genomic sequence could not be extracted:\s+(\d+)$",
        'total_c': r"^Total number of C's analysed:\s+(\d+)$",
        'meth_cpg': r"^Total methylated C's in CpG context:\s+(\d+)",
        'meth_chg': r"^Total methylated C's in CHG context:\s+(\d+)",
        'meth_chh': r"^Total methylated C's in CHH context:\s+(\d+)",
        'unmeth_cpg': r"^Total unmethylated C's in CpG context:\s+(\d+)",
        'unmeth_chg': r"^Total unmethylated C's in CHG context:\s+(\d+)",
        'unmeth_chh': r"^Total unmethylated C's in CHH context:\s+(\d+)",
        'percent_cpg_meth': r"^C methylated in CpG context:\s+([\d\.]+)%",
        'percent_chg_meth': r"^C methylated in CHG context:\s+([\d\.]+)%",
        'percent_chh_meth': r"^C methylated in CHH context:\s+([\d\.]+)%",
        'strand_ot': r"^CT(?:\/GA)?\/CT:\s+(\d+)\s+\(\(converted\) top strand\)$",
        'strand_ctot': r"^GA(?:\/CT)?\/CT:\s+(\d+)\s+\(complementary to \(converted\) top strand\)$",
        'strand_ctob': r"^GA(?:\/CT)?\/GA:\s+(\d+)\s+\(complementary to \(converted\) bottom strand\)$",
        'strand_ob': r"^CT(?:\/GA)?\/GA:\s+(\d+)\s+\(\(converted\) bottom strand\)$",
        'strand_directional': r"^Option '--(directional)' specified \(default mode\): alignments to complementary strands \(CTOT, CTOB\) were ignored \(i.e. not performed\)$"
    },
    'dedup': {
        # 'aligned_reads' overwrites previous, but I trust this more
        'aligned_reads': r"^Total number of alignments analysed in .+:\s+(\d+)$",
        'dup_reads': r"^Total number duplicated alignments removed:\s+(\d+)",
        'dup_reads_percent': r"^Total number duplicated alignments removed:\s+\d+\s+\(([\d\.]+)%\)",
        'dedup_reads': r"^Total count of deduplicated leftover sequences:\s+(\d+)",
        'dedup_reads_percent': r"^Total count of deduplicated leftover sequences:\s+\d+\s+\(([\d\.]+)% of total\)"
    },
    'methextract': {
        # These calls are typically done after deduplication.
        # Overwrites what was found in the alignment report, but now
        # after deduplication (if that was done)
        'total_c': r"^Total number of C's analysed:\s+(\d+)$",
        'meth_cpg': r"^Total methylated C's in CpG context:\s+(\d+)",
        'meth_chg': r"^Total methylated C's in CHG context:\s+(\d+)",
        'meth_chh': r"^Total methylated C's in CHH context:\s+(\d+)",
        'unmeth_cpg': r"^Total C to T conversions in CpG context:\s+(\d+)",
        'unmeth_chg': r"^Total C to T conversions in CHG context:\s+(\d+)",
        'unmeth_chh': r"^Total C to T conversions in CHH context:\s+(\d+)",
        'percent_cpg_meth': r"^C methylated in CpG context:\s+([\d\.]+)%",
        'percent_chg_meth': r"^C methylated in CHG context:\s+([\d\.]+)%",
        'percent_chh_meth': r"^C methylated in CHH context:\s+([\d\.]+)%"
    }
}

class MultiqcModule(multiqc.BaseMultiqcModule):

    def __init__(self):

        # Initialise the parent object
        super(MultiqcModule, self).__init__(log)

        # Static variables
        self.name = "Bismark"
        self.anchor = "bismark"
        self.intro = '<p><a href="http://www.bioinformatics.babraham.ac.uk/projects/bismark/" target="_blank">Bismark</a> \
            is a tool to map bisulfite converted sequence reads and determine cytosine methylation states.</p>'

        # Find and load any Bismark reports
        self.bismark_data = {
            'alignment': {},
            'dedup': {},
            'methextract': {},
            'merged': defaultdict(lambda: dict())
        }
        # Find and parse bismark alignment reports
        for f in self.find_log_files(['_PE_report.txt', '_SE_report.txt']):
            parsed_data = self.parse_bismark_report(f['f'], regexes['alignment'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['alignment']:
                    log.debug("Duplicate alignment sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['alignment'][f['s_name']] = parsed_data
        
        # Find and parse bismark deduplication reports
        for f in self.find_log_files('.deduplication_report.txt'):
            parsed_data = self.parse_bismark_report(f['f'], regexes['dedup'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['dedup']:
                    log.debug("Duplicate deduplication sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['dedup'][f['s_name']] = parsed_data
        
        # Find and parse bismark methylation extractor reports
        for f in self.find_log_files('_splitting_report.txt'):
            parsed_data = self.parse_bismark_report(f['f'], regexes['methextract'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['methextract']:
                    log.debug("Duplicate methylation extraction sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['methextract'][f['s_name']] = parsed_data
        
        # Merge the keys into a merged dict (in order)
        for t in ['alignment','dedup','methextract']:
            for sn in self.bismark_data[t].keys():
                for k, v in self.bismark_data[t][sn].items():
                    self.bismark_data['merged'][sn][k] = v
        
        if len(self.bismark_data['merged']) == 0:
            log.debug("Could not find any reports in {}".format(config.analysis_dir))
            raise UserWarning

        log.info("Found {} reports".format(len(self.bismark_data['merged'])))

        # Write parsed report data to a file
        self.write_csv_file(self.bismark_data['merged'], 'multiqc_bismark.txt')

        self.sections = list()
        
        # Basic Stats Table
        self.bismark_stats_table()

        # Section 1 - Column chart of alignment stats
        self.sections.append({
            'name': 'Alignment Rates',
            'anchor': 'bismark-alignment',
            'content': self.bismark_alignment_chart()
        })

        # Section 2 - Methylation percentages
        self.sections.append({
            'name': 'Cytosine Methylation',
            'anchor': 'bismark-methylation',
            'content': self.bismark_methlyation_chart()
        })

        # Section 3 - Strand Alignments
        self.sections.append({
            'name': 'Strand Alignment',
            'anchor': 'bismark-strands',
            'content': self.bismark_strand_chart()
        })

    def parse_bismark_report(self, report, regexes):
        """ Search a bismark report with a set of regexes """
        parsed_data = {}
        for k, r in regexes.items():
            r_search = re.search(r, report, re.MULTILINE)
            if r_search:
                try:
                    parsed_data[k] = float(r_search.group(1))
                except ValueError:
                    parsed_data[k] = r_search.group(1) # NaN
        if len(parsed_data) == 0: return None
        return parsed_data

    def bismark_stats_table(self):
        """ Take the parsed stats from the Bismark reports and add them to the
        basic stats table at the top of the report """

        # Use several try blocks in case one of the report types is missing
        # If data is missing, header rows won't be added
        
        # Alignment / Meth Extract data
        try:
            for sn, data in self.bismark_data['merged'].items():
                config.general_stats['rows'][sn]['percent_cpg_meth'] = '<td class="text-right">{:.1f}%</td>'.format(data['percent_cpg_meth'])
                config.general_stats['rows'][sn]['total_c'] = '<td class="text-right">{:.1f}</td>'.format(data['total_c']/1000000)
            config.general_stats['headers']['percent_cpg_meth'] = '<th class="chroma-col" data-chroma-scale="Greens" data-chroma-min="0"><span data-toggle="tooltip" title="Bismark: % Cytosines methylated in CpG context (alignment)">% Meth</span></th>'
            config.general_stats['headers']['total_c'] = '<th class="chroma-col" data-chroma-scale="Purples" data-chroma-min="0"><span data-toggle="tooltip" title="Bismark: Total number of C\'s analysed, in millions (alignment)">M C\'s</span></th>'
        except KeyError:
            pass
        
        # Deduplication data
        try:
            for sn, data in self.bismark_data['merged'].items():
                config.general_stats['rows'][sn]['bismark_dedup_reads_percent'] = '<td class="text-right">{:.1f}%</td>'.format(data['dup_reads_percent'])
                config.general_stats['rows'][sn]['bismark_dedup_reads'] = '<td class="text-right">{:.1f}</td>'.format(data['dedup_reads']/1000000)<commit_msg>General stats table tweak for bismark.<commit_after>#!/usr/bin/env python

""" MultiQC module to parse output from Bismark """

from __future__ import print_function
from collections import defaultdict, OrderedDict
import io
import json
import logging
import os
import re

import multiqc
from multiqc import config

# Initialise the logger
log = logging.getLogger('MultiQC : {0:<14}'.format('Bismark'))

# Log parsing regexes
regexes = {
    'alignment': {
        'total_reads': r"^Sequence(?:s| pairs) analysed in total:\s+(\d+)$",
        'aligned_reads': r"^Number of(?: paired-end)? alignments with a unique best hit(?: from the different alignments)?:\s+(\d+)$",
        'no_alignments': r"^Sequence(?:s| pairs) with no alignments under any condition:\s+(\d+)$",
        'ambig_reads': r"^Sequence(?:s| pairs) did not map uniquely:\s+(\d+)$",
        'discarded_reads': r"^Sequence(?:s| pairs) which were discarded because genomic sequence could not be extracted:\s+(\d+)$",
        'total_c': r"^Total number of C's analysed:\s+(\d+)$",
        'meth_cpg': r"^Total methylated C's in CpG context:\s+(\d+)",
        'meth_chg': r"^Total methylated C's in CHG context:\s+(\d+)",
        'meth_chh': r"^Total methylated C's in CHH context:\s+(\d+)",
        'unmeth_cpg': r"^Total unmethylated C's in CpG context:\s+(\d+)",
        'unmeth_chg': r"^Total unmethylated C's in CHG context:\s+(\d+)",
        'unmeth_chh': r"^Total unmethylated C's in CHH context:\s+(\d+)",
        'percent_cpg_meth': r"^C methylated in CpG context:\s+([\d\.]+)%",
        'percent_chg_meth': r"^C methylated in CHG context:\s+([\d\.]+)%",
        'percent_chh_meth': r"^C methylated in CHH context:\s+([\d\.]+)%",
        'strand_ot': r"^CT(?:\/GA)?\/CT:\s+(\d+)\s+\(\(converted\) top strand\)$",
        'strand_ctot': r"^GA(?:\/CT)?\/CT:\s+(\d+)\s+\(complementary to \(converted\) top strand\)$",
        'strand_ctob': r"^GA(?:\/CT)?\/GA:\s+(\d+)\s+\(complementary to \(converted\) bottom strand\)$",
        'strand_ob': r"^CT(?:\/GA)?\/GA:\s+(\d+)\s+\(\(converted\) bottom strand\)$",
        'strand_directional': r"^Option '--(directional)' specified \(default mode\): alignments to complementary strands \(CTOT, CTOB\) were ignored \(i.e. not performed\)$"
    },
    'dedup': {
        # 'aligned_reads' overwrites previous, but I trust this more
        'aligned_reads': r"^Total number of alignments analysed in .+:\s+(\d+)$",
        'dup_reads': r"^Total number duplicated alignments removed:\s+(\d+)",
        'dup_reads_percent': r"^Total number duplicated alignments removed:\s+\d+\s+\(([\d\.]+)%\)",
        'dedup_reads': r"^Total count of deduplicated leftover sequences:\s+(\d+)",
        'dedup_reads_percent': r"^Total count of deduplicated leftover sequences:\s+\d+\s+\(([\d\.]+)% of total\)"
    },
    'methextract': {
        # These calls are typically done after deduplication.
        # Overwrites what was found in the alignment report, but now
        # after deduplication (if that was done)
        'total_c': r"^Total number of C's analysed:\s+(\d+)$",
        'meth_cpg': r"^Total methylated C's in CpG context:\s+(\d+)",
        'meth_chg': r"^Total methylated C's in CHG context:\s+(\d+)",
        'meth_chh': r"^Total methylated C's in CHH context:\s+(\d+)",
        'unmeth_cpg': r"^Total C to T conversions in CpG context:\s+(\d+)",
        'unmeth_chg': r"^Total C to T conversions in CHG context:\s+(\d+)",
        'unmeth_chh': r"^Total C to T conversions in CHH context:\s+(\d+)",
        'percent_cpg_meth': r"^C methylated in CpG context:\s+([\d\.]+)%",
        'percent_chg_meth': r"^C methylated in CHG context:\s+([\d\.]+)%",
        'percent_chh_meth': r"^C methylated in CHH context:\s+([\d\.]+)%"
    }
}

class MultiqcModule(multiqc.BaseMultiqcModule):

    def __init__(self):

        # Initialise the parent object
        super(MultiqcModule, self).__init__(log)

        # Static variables
        self.name = "Bismark"
        self.anchor = "bismark"
        self.intro = '<p><a href="http://www.bioinformatics.babraham.ac.uk/projects/bismark/" target="_blank">Bismark</a> \
            is a tool to map bisulfite converted sequence reads and determine cytosine methylation states.</p>'

        # Find and load any Bismark reports
        self.bismark_data = {
            'alignment': {},
            'dedup': {},
            'methextract': {},
            'merged': defaultdict(lambda: dict())
        }
        # Find and parse bismark alignment reports
        for f in self.find_log_files(['_PE_report.txt', '_SE_report.txt']):
            parsed_data = self.parse_bismark_report(f['f'], regexes['alignment'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['alignment']:
                    log.debug("Duplicate alignment sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['alignment'][f['s_name']] = parsed_data
        
        # Find and parse bismark deduplication reports
        for f in self.find_log_files('.deduplication_report.txt'):
            parsed_data = self.parse_bismark_report(f['f'], regexes['dedup'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['dedup']:
                    log.debug("Duplicate deduplication sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['dedup'][f['s_name']] = parsed_data
        
        # Find and parse bismark methylation extractor reports
        for f in self.find_log_files('_splitting_report.txt'):
            parsed_data = self.parse_bismark_report(f['f'], regexes['methextract'])
            if parsed_data is not None:
                if f['s_name'] in self.bismark_data['methextract']:
                    log.debug("Duplicate methylation extraction sample log found! Overwriting: {}".format(f['s_name']))
                self.bismark_data['methextract'][f['s_name']] = parsed_data
        
        # Merge the keys into a merged dict (in order)
        for t in ['alignment','dedup','methextract']:
            for sn in self.bismark_data[t].keys():
                for k, v in self.bismark_data[t][sn].items():
                    self.bismark_data['merged'][sn][k] = v
        
        if len(self.bismark_data['merged']) == 0:
            log.debug("Could not find any reports in {}".format(config.analysis_dir))
            raise UserWarning

        log.info("Found {} reports".format(len(self.bismark_data['merged'])))

        # Write parsed report data to a file
        self.write_csv_file(self.bismark_data['merged'], 'multiqc_bismark.txt')

        self.sections = list()
        
        # Basic Stats Table
        self.bismark_stats_table()

        # Section 1 - Column chart of alignment stats
        self.sections.append({
            'name': 'Alignment Rates',
            'anchor': 'bismark-alignment',
            'content': self.bismark_alignment_chart()
        })

        # Section 2 - Methylation percentages
        self.sections.append({
            'name': 'Cytosine Methylation',
            'anchor': 'bismark-methylation',
            'content': self.bismark_methlyation_chart()
        })

        # Section 3 - Strand Alignments
        self.sections.append({
            'name': 'Strand Alignment',
            'anchor': 'bismark-strands',
            'content': self.bismark_strand_chart()
        })

    def parse_bismark_report(self, report, regexes):
        """ Search a bismark report with a set of regexes """
        parsed_data = {}
        for k, r in regexes.items():
            r_search = re.search(r, report, re.MULTILINE)
            if r_search:
                try:
                    parsed_data[k] = float(r_search.group(1))
                except ValueError:
                    parsed_data[k] = r_search.group(1) # NaN
        if len(parsed_data) == 0: return None
        return parsed_data

    def bismark_stats_table(self):
        """ Take the parsed stats from the Bismark reports and add them to the
        basic stats table at the top of the report """

        # Use several try blocks in case one of the report types is missing
        # If data is missing, header rows won't be added
        
        # Alignment / Meth Extract data
        try:
            for sn, data in self.bismark_data['merged'].items():
                config.general_stats['rows'][sn]['percent_cpg_meth'] = '<td class="text-right">{:.1f}%</td>'.format(data['percent_cpg_meth'])
                config.general_stats['rows'][sn]['total_c'] = '<td class="text-right">{:.1f}</td>'.format(data['total_c']/1000000)
            config.general_stats['headers']['percent_cpg_meth'] = '<th class="chroma-col" data-chroma-scale="Greens" data-chroma-max="100" data-chroma-min="0"><span data-toggle="tooltip" title="Bismark: % Cytosines methylated in CpG context (alignment)">% Meth</span></th>'
            config.general_stats['headers']['total_c'] = '<th class="chroma-col" data-chroma-scale="Purples" data-chroma-min="0"><span data-toggle="tooltip" title="Bismark: Total number of C\'s analysed, in millions (alignment)">M C\'s</span></th>'
        except KeyError:
            pass
        
        # Deduplication data
        try:
            for sn, data in self.bismark_data['merged'].items():
                config.general_stats['rows'][sn]['bismark_dedup_reads_percent'] = '<td class="text-right">{:.1f}%</td>'.format(data['dup_reads_percent'])
                config.general_stats['rows'][sn]['bismark_dedup_reads'] = '<td class="text-right">{:.1f}</td>'.format(data['dedup_reads']/1000000)
{'fraction_non_alphanumeric': 0.1604043983532232, 'fraction_numerical': 0.003283130960446089, 'mean_word_length': 2.7305598755832037, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4596931', 'n_tokens_mistral': 6453, 'n_tokens_neox': 5967, 'n_words': 1449}