Find Qualimap RNASeq reports and parse their data
776694
def parse_reports(self):
    

    self.qualimap_rnaseq_genome_results = dict()
    regexes = {
        'reads_aligned': r"read(?:s| pairs) aligned\s*=\s*([\d,]+)",
        'total_alignments': r"total alignments\s*=\s*([\d,]+)",
        'non_unique_alignments': r"non-unique alignments\s*=\s*([\d,]+)",
        'reads_aligned_genes': r"aligned to genes\s*=\s*([\d,]+)",
        'ambiguous_alignments': r"ambiguous alignments\s*=\s*([\d,]+)",
        'not_aligned': r"not aligned\s*=\s*([\d,]+)",
        '5_3_bias': r"5'-3' bias\s*=\s*([\d,\.]+)$",
        'reads_aligned_exonic': r"exonic\s*=\s*([\d,]+)",
        'reads_aligned_intronic': r"intronic\s*=\s*([\d,]+)",
        'reads_aligned_intergenic': r"intergenic\s*=\s*([\d,]+)",
        'reads_aligned_overlapping_exon': r"overlapping exon\s*=\s*([\d,]+)",
    }
    for f in self.find_log_files('qualimap/rnaseq/rnaseq_results'):
        d = dict()

        # Get the sample name
        s_name_regex = re.search(r"bam file\s*=\s*(.+)", f['f'], re.MULTILINE)
        if s_name_regex:
            d['bam_file'] = s_name_regex.group(1)
            s_name = self.clean_s_name(d['bam_file'], f['root'])
        else:
            log.warn("Couldn't find an input filename in genome_results file {}/{}".format(f['root'], f['fn']))
            return None

        # Check for and 'fix' European style decimal places / thousand separators
        comma_regex = re.search(r"exonic\s*=\s*[\d\.]+ \(\d{1,3},\d+%\)", f['f'], re.MULTILINE)
        if comma_regex:
            log.debug("Trying to fix European comma style syntax in Qualimap report {}/{}".format(f['root'], f['fn']))
            f['f'] = f['f'].replace('.','')
            f['f'] = f['f'].replace(',','.')

        # Go through all numeric regexes
        for k, r in regexes.items():
            r_search = re.search(r, f['f'], re.MULTILINE)
            if r_search:
                try:
                    d[k] = float(r_search.group(1).replace(',',''))
                except UnicodeEncodeError:
                    # Qualimap reports infinity (\u221e) when 3' bias denominator is zero
                    pass
                except ValueError:
                    d[k] = r_search.group(1)

        # Add to general stats table
        for k in ['5_3_bias', 'reads_aligned']:
            try:
                self.general_stats_data[s_name][k] = d[k]
            except KeyError:
                pass

        # Save results
        if s_name in self.qualimap_rnaseq_genome_results:
            log.debug("Duplicate genome results sample name found! Overwriting: {}".format(s_name))
        self.qualimap_rnaseq_genome_results[s_name] = d
        self.add_data_source(f, s_name=s_name, section='rna_genome_results')


    #### Coverage profile
    self.qualimap_rnaseq_cov_hist = dict()
    for f in self.find_log_files('qualimap/rnaseq/coverage', filehandles=True):
        s_name = self.get_s_name(f)
        d = dict()
        for l in f['f']:
            if l.startswith('#'):
                continue
            coverage, count = l.split(None, 1)
            coverage = int(round(float(coverage)))
            count = float(count)
            d[coverage] = count

        if len(d) == 0:
            log.debug("Couldn't parse contents of coverage histogram file {}".format(f['fn']))
            return None

        # Save results
        if s_name in self.qualimap_rnaseq_cov_hist:
            log.debug("Duplicate coverage histogram sample name found! Overwriting: {}".format(s_name))
        self.qualimap_rnaseq_cov_hist[s_name] = d
        self.add_data_source(f, s_name=s_name, section='rna_coverage_histogram')

    # Filter to strip out ignored sample names
    self.qualimap_rnaseq_genome_results = self.ignore_samples(self.qualimap_rnaseq_genome_results)
    self.qualimap_rnaseq_cov_hist = self.ignore_samples(self.qualimap_rnaseq_cov_hist)

    #### Plots

    # Genomic Origin Bar Graph
    # NB: Ignore 'Overlapping Exon' in report - these make the numbers add up to > 100%
    if len(self.qualimap_rnaseq_genome_results) > 0:
        gorigin_cats = OrderedDict()
        gorigin_cats['reads_aligned_exonic'] = {'name': 'Exonic'}
        gorigin_cats['reads_aligned_intronic'] = {'name': 'Intronic'}
        gorigin_cats['reads_aligned_intergenic'] = {'name': 'Intergenic'}
        gorigin_pconfig = {
            'id': 'qualimap_genomic_origin',
            'title': 'Qualimap RNAseq: Genomic Origin',
            'ylab': 'Number of reads',
            'cpswitch_c_active': False
        }
        genomic_origin_helptext = '''
        There are currently three main approaches to map reads to transcripts in an
        RNA-seq experiment: mapping reads to a reference genome to identify expressed
        transcripts that are annotated (and discover those that are unknown), mapping
        reads to a reference transcriptome, and <i>de novo</i> assembly of transcript
        sequences (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).

        For RNA-seq QC analysis, QualiMap can be used to assess alignments produced by
        the first of these approaches. For input, it requires a GTF annotation file
        along with a reference genome, which can be used to reconstruct the exon
        structure of known transcripts. This allows mapped reads to be grouped by
        whether they originate in an exonic region (for QualiMap, this may include
        5&#8242; and 3&#8242; UTR regions as well as protein-coding exons), an intron,
        or an intergenic region (see the <a href="http://qualimap.bioinfo.cipf.es/doc_html/index.html"
        target="_blank">Qualimap 2 documentation</a>).

        The inferred genomic origins of RNA-seq reads are presented here as a bar graph
        showing either the number or percentage of mapped reads in each read dataset
        that have been assigned to each type of genomic region. This graph can be used
        to assess the proportion of useful reads in an RNA-seq experiment. That
        proportion can be reduced by the presence of intron sequences, especially if
        depletion of ribosomal RNA was used during sample preparation (<a href="https://doi.org/10.1038/nrg3642"
        target="_blank">Sims et al. 2014</a>). It can also be reduced by off-target
        transcripts, which are detected in greater numbers at the sequencing depths
        needed to detect poorly-expressed transcripts (<a href="https://doi.org/10.1101/gr.124321.111"
        target="_blank">Tarazona et al. 2011</a>).'''
        self.add_section (
            name = 'Genomic origin of reads',
            anchor = 'qualimap-reads-genomic-origin',
            description = 'Classification of mapped reads as originating in exonic, intronic or intergenic regions. These can be displayed as either the number or percentage of mapped reads.',
            helptext = genomic_origin_helptext,
            plot = bargraph.plot(self.qualimap_rnaseq_genome_results, gorigin_cats, gorigin_pconfig)
        )

    if len(self.qualimap_rnaseq_cov_hist) > 0:
        coverage_profile_helptext = '''
        There are currently three main approaches to map reads to transcripts in an
        RNA-seq experiment: mapping reads to a reference genome to identify expressed
        transcripts that are annotated (and discover those that are unknown), mapping
        reads to a reference transcriptome, and <i>de novo</i> assembly of transcript
        sequences (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).

        For RNA-seq QC analysis, QualiMap can be used to assess alignments produced by
        the first of these approaches. For input, it requires a GTF annotation file
        along with a reference genome, which can be used to reconstruct the exon
        structure of known transcripts. QualiMap uses this information to calculate the
        depth of coverage along the length of each annotated transcript. For a set of
        reads mapped to a transcript, the depth of coverage at a given base position is
        the number of high-quality reads that map to the transcript at that position
        (<a href="https://doi.org/10.1038/nrg3642" target="_blank">Sims et al. 2014</a>).

        QualiMap calculates coverage depth at every base position of each annotated
        transcript. To enable meaningful comparison between transcripts, base positions
        are rescaled to relative positions expressed as percentage distance along each
        transcript (*0%, 1%, &#8230;, 99%*). For the set of transcripts with at least
        one mapped read, QualiMap plots the cumulative mapped-read depth (y-axis) at
        each relative transcript position (x-axis). This plot shows the gene coverage
        profile across all mapped transcripts for each read dataset. It provides a
        visual way to assess positional biases, such as an accumulation of mapped reads
        at the 3&#8242; end of transcripts, which may indicate poor RNA quality in the
        original sample (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).'''
        self.add_section (
            name = 'Gene Coverage Profile',
            anchor = 'qualimap-genome-fraction-coverage',
            description = 'Mean distribution of coverage depth across the length of all mapped transcripts.',
            helptext = coverage_profile_helptext,
            plot = linegraph.plot(self.qualimap_rnaseq_cov_hist, {
                'id': 'qualimap_gene_coverage_profile',
                'title': 'Qualimap RNAseq: Coverage Profile Along Genes (total)',
                'ylab': 'Coverage',
                'xlab': 'Transcript Position (%)',
                'ymin': 0,
                'xmin': 0,
                'xmax': 100,
                'tt_label': '<b>{point.x} bp</b>: {point.y:.0f}%',
            })
        )


    #### General Stats
    self.general_stats_headers['5_3_bias'] = {
        'title': "5'-3' bias",
        'format': '{:,.2f}',
    }
    self.general_stats_headers['reads_aligned'] = {
        'title': '{} Aligned'.format(config.read_count_prefix),
        'description': 'Reads Aligned ({})'.format(config.read_count_desc),
        'min': 0,
        'scale': 'RdBu',
        'shared_key': 'read_count',
        'modify': lambda x: x * config.read_count_multiplier
    }

    # Return the number of reports we found
    return len(self.qualimap_rnaseq_genome_results.keys())
927362
def Parse(self):
    """"""
    for data in self.Query(self.EVENTS_QUERY):
      (timestamp, agent_bundle_identifier, agent_name, url, sender,
       sender_address, type_number, title, referrer, referrer_alias) = data
      yield [
          timestamp, "OSX_QUARANTINE", url, referrer, title, agent_name,
          agent_bundle_identifier, sender, sender_address, type_number,
          referrer_alias
      ]
966584
def get_summary_metrics(self, align_metrics, dup_metrics,
            insert_metrics=None, hybrid_metrics=None, vrn_vals=None,
            rnaseq_metrics=None):
        
        with open(align_metrics) as in_handle:
            align_vals = self._parse_align_metrics(in_handle)
        if dup_metrics:
            with open(dup_metrics) as in_handle:
                dup_vals = self._parse_dup_metrics(in_handle)
        else:
            dup_vals = {}
        (insert_vals, hybrid_vals, rnaseq_vals) = (None, None, None)
        if insert_metrics and file_exists(insert_metrics):
            with open(insert_metrics) as in_handle:
                insert_vals = self._parse_insert_metrics(in_handle)
        if hybrid_metrics and file_exists(hybrid_metrics):
            with open(hybrid_metrics) as in_handle:
                hybrid_vals = self._parse_hybrid_metrics(in_handle)
        if rnaseq_metrics and file_exists(rnaseq_metrics):
            with open(rnaseq_metrics) as in_handle:
                rnaseq_vals = self._parse_rnaseq_metrics(in_handle)

        return self._tabularize_metrics(align_vals, dup_vals, insert_vals,
                hybrid_vals, vrn_vals, rnaseq_vals)
1005747
def get_stats(self):
    """"""
    ostr = ''
    errtotal = self.deletions['total']+self.insertions['total']+self.mismatches
    ostr += "ALIGNMENT_COUNT\t"+str(self.alignment_count)+"\n"
    ostr += "ALIGNMENT_BASES\t"+str(self.alignment_length)+"\n"
    ostr += "ANY_ERROR\t"+str(errtotal)+"\n"
    ostr += "MISMATCHES\t"+str(self.mismatches)+"\n"
    ostr += "ANY_DELETION\t"+str(self.deletions['total'])+"\n"
    ostr += "COMPLETE_DELETION\t"+str(self.deletions['specific'])+"\n"
    ostr += "HOMOPOLYMER_DELETION\t"+str(self.deletions['homopolymer'])+"\n"
    ostr += "ANY_INSERTION\t"+str(self.insertions['total'])+"\n"
    ostr += "COMPLETE_INSERTION\t"+str(self.insertions['specific'])+"\n"
    ostr += "HOMOPOLYMER_INSERTION\t"+str(self.insertions['homopolymer'])+"\n"
    return ostr
429033
def read_RNAplfold(tmpdir, maxlen=None, seq_align="start", pad_with="E"):
    
    assert pad_with in {"P", "H", "I", "M", "E"}

    def read_profile(tmpdir, P):
        return [values.strip().split("\t")
                for seq_name, values in iter_fasta("{tmp}/{P}_profile.fa".format(tmp=tmpdir, P=P))]

    def nelem(P, pad_width):
        """get the right neutral element
        """
        return 1 if P is pad_with else 0

    arr_hime = np.array([pad_sequences(read_profile(tmpdir, P),
                                       value=[nelem(P, pad_with)],
                                       align=seq_align,
                                       maxlen=maxlen)
                         for P in RNAplfold_PROFILES_EXECUTE], dtype="float32")

    # add the pairness column
    arr_p = 1 - arr_hime.sum(axis=0)[np.newaxis]
    arr = np.concatenate((arr_p, arr_hime))

    # reshape to: seq, seq_length, num_channels
    arr = np.moveaxis(arr, 0, 2)
    return arr
303571
def parse_reports(self):
    

    # Set up vars
    self.read_dups = dict()

    # Go through files and parse data
    for f in self.find_log_files('rseqc/read_duplication_pos'):
        if f['f'].startswith('Occurrence	UniqReadNumber'):
            if f['s_name'] in self.read_dups:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f, section='read_duplication')
            self.read_dups[f['s_name']] = OrderedDict()
            for l in f['f'].splitlines():
                s = l.split()
                try:
                    if int(s[0]) <= 500:
                        self.read_dups[f['s_name']][int(s[0])] = int(s[1])
                except:
                    pass

    # Filter to strip out ignored sample names
    self.read_dups = self.ignore_samples(self.read_dups)

    if len(self.read_dups) > 0:

        # Add line graph to section
        pconfig = {
            'smooth_points': 200,
            'id': 'rseqc_read_dups_plot',
            'title': 'RSeQC: Read Duplication',
            'ylab': 'Number of Reads (log10)',
            'xlab': "Occurrence of read",
            'yLog': True,
            'tt_label': "<strong>{point.x} occurrences</strong>: {point.y} reads",
        }

        self.add_section (
            name = 'Read Duplication',
            anchor = 'rseqc-read_dups',
            description = '<a href="http://rseqc.sourceforge.net/#read-duplication-py" target="_blank">read_duplication.py</a>' \
                " calculates how many alignment positions have a certain number of exact duplicates."\
                " Note - plot truncated at 500 occurrences and binned.</p>",
            plot = linegraph.plot(self.read_dups, pconfig)
        )

    # Return number of samples found
    return len(self.read_dups)
791048
def parse_reports(self):
    

    # Set up vars
    self.picard_rrbs_metrics = dict()

    # Go through logs and find Metrics
    for f in self.find_log_files('picard/rrbs_metrics', filehandles=True):
        parsed_data = dict()
        s_name = None
        keys = None
        for l in f['f']:
            # New log starting
            if 'CollectRrbsMetrics' in l and 'INPUT' in l:
                s_name = None
                keys = None
                # Pull sample name from input
                fn_search = re.search(r"INPUT(?:=|\s+)(\[?[^\s]+\]?)", l, flags=re.IGNORECASE)
                if fn_search:
                    s_name = os.path.basename(fn_search.group(1).strip('[]'))
                    s_name = self.clean_s_name(s_name, f['root'])
                    parsed_data[s_name] = dict()

            if s_name is not None:
                if 'RrbsSummaryMetrics' in l and '## METRICS CLASS' in l:
                    keys = f['f'].readline().strip("\n").split("\t")
                elif keys:
                    vals = l.strip("\n").split("\t")
                    if len(vals) == len(keys):
                        for i, k in enumerate(keys):
                            try:
                                parsed_data[s_name][k] = float(vals[i])
                            except ValueError:
                                parsed_data[s_name][k] = vals[i]
                    else:
                        s_name = None
                        keys = None

        # Remove empty dictionaries
        for s_name in list(parsed_data.keys()):
            if len(parsed_data[s_name]) == 0:
                parsed_data.pop(s_name, None)

        # Collect parsed data
        for s_name in parsed_data.keys():
            if s_name in self.picard_rrbs_metrics:
                log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
            self.add_data_source(f, s_name, section='RrbsSummaryMetrics')
            self.picard_rrbs_metrics[s_name] = parsed_data[s_name]

    # Filter to strip out ignored sample names
    self.picard_rrbs_metrics = self.ignore_samples(self.picard_rrbs_metrics)

    if len(self.picard_rrbs_metrics) > 0:

        # Write parsed data to a file
        self.write_data_file(self.picard_rrbs_metrics, 'multiqc_picard_RrbsSummaryMetrics')

        # Add to general stats table
        self.general_stats_headers['PCT_CPG_BASES_CONVERTED'] = {
            'title': '% CpG Methylated',
            'description': 'Percentage of times a CpG cytosine was converted',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'format': '{:,.0f}',
            'scale': 'RdYlGn-rev',
            'modify': lambda x: 100 - self.multiply_hundred(x)
        }
        self.general_stats_headers['PCT_NON_CPG_BASES_CONVERTED'] = {
            'title': '% Non-CpG Methylated',
            'description': 'Percentage of times a non-CpG cytosine was converted',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'format': '{:,.0f}',
            'scale': 'RdYlGn',
            'modify': lambda x: 100 - self.multiply_hundred(x)
        }
        self.general_stats_headers['MEDIAN_CPG_COVERAGE'] = {
            'title': 'Median CpG Cov',
            'description': 'Median coverage of CpG sites',
            'min': 0
        }
        for s_name in self.picard_rrbs_metrics:
            if s_name not in self.general_stats_data:
                self.general_stats_data[s_name] = dict()
            self.general_stats_data[s_name].update( self.picard_rrbs_metrics[s_name] )



        # Make the bar plot of converted bases
        pdata_cpg = dict()
        pdata_noncpg = dict()
        for s_name in self.picard_rrbs_metrics.keys():
            pdata_cpg[s_name] = dict()
            pdata_cpg[s_name]['converted'] = self.picard_rrbs_metrics[s_name]['CPG_BASES_CONVERTED']
            pdata_cpg[s_name]['not_converted'] = self.picard_rrbs_metrics[s_name]['CPG_BASES_SEEN'] - self.picard_rrbs_metrics[s_name]['CPG_BASES_CONVERTED']
            pdata_noncpg[s_name] = dict()
            pdata_noncpg[s_name]['converted'] = self.picard_rrbs_metrics[s_name]['NON_CPG_BASES']
            pdata_noncpg[s_name]['not_converted'] = self.picard_rrbs_metrics[s_name]['NON_CPG_BASES'] - self.picard_rrbs_metrics[s_name]['NON_CPG_CONVERTED_BASES']

        keys = OrderedDict()
        keys['not_converted'] = {'name': 'Unconverted Bases (Methylated)'}
        keys['converted'] = {'name': 'Converted Bases (Unmethylated)'}

        # Config for the plot
        pconfig = {
            'id': 'picard_rrbs_converted_bases_plot',
            'title': 'Picard: Converted Bases',
            'ylab': '# CpG Bases',
            'cpswitch_counts_label': 'Number of Bases',
            'data_labels': [
                {'name': 'CpG', 'ylab': '# CpG Bases'},
                {'name': 'Non-CpG', 'ylab': '# Non-CpG Bases'}
            ]
        }

        self.add_section (
            name = 'RRBS Converted Bases',
            anchor = 'picard-rrbssummary-convertedbases',
            plot = bargraph.plot([pdata_cpg, pdata_noncpg], [keys, keys], pconfig)
        )


        # Make the bar plot of processed reads
        pdata = dict()
        for s_name in self.picard_rrbs_metrics.keys():
            pdata[s_name] = dict()
            pdata[s_name]['with_no_cpg'] = self.picard_rrbs_metrics[s_name]['READS_WITH_NO_CPG']
            pdata[s_name]['ignored_short'] = self.picard_rrbs_metrics[s_name]['READS_IGNORED_SHORT']
            pdata[s_name]['ignored_mismatches'] = self.picard_rrbs_metrics[s_name]['READS_IGNORED_MISMATCHES']
            pdata[s_name]['not_ignored'] = (
                self.picard_rrbs_metrics[s_name]['READS_ALIGNED'] -
                pdata[s_name]['with_no_cpg'] -
                pdata[s_name]['ignored_short'] -
                pdata[s_name]['ignored_mismatches']
            )

        keys = OrderedDict()
        keys['not_ignored'] = {'name': 'Utilised reads'}
        keys['with_no_cpg'] = {'name': 'Ignored (no CpG sites)'}
        keys['ignored_short'] = {'name': 'Ignored (too short)'}
        keys['ignored_mismatches'] = {'name': 'Ignored (exceeded mismatch threshold)'}

        # Config for the plot
        pconfig = {
            'id': 'picard_rrbs_ignored_reads_plot',
            'title': 'Picard: RRBS Read Counts',
            'ylab': '# Reads',
            'cpswitch_counts_label': 'Number of Reads'
        }

        self.add_section (
            name = 'RRBS Read Counts',
            anchor = 'picard-rrbssummary-readcounts',
            plot = bargraph.plot(pdata, keys, pconfig)
        )

    # Return the number of detected samples to the parent module
    return len(self.picard_rrbs_metrics)
188837
def parse_reports(self):
    

    # Set up vars
    self.bam_stat_data = dict()
    regexes = {
        'total_records': r"Total records:\s*(\d+)",
        'qc_failed': r"QC failed:\s*(\d+)",
        'optical_pcr_duplicate': r"Optical/PCR duplicate:\s*(\d+)",
        'non_primary_hits': r"Non primary hits\s*(\d+)",
        'unmapped_reads': r"Unmapped reads:\s*(\d+)",
        'mapq_lt_mapq_cut_non-unique': r"mapq < mapq_cut \(non-unique\):\s*(\d+)",
        'mapq_gte_mapq_cut_unique': r"mapq >= mapq_cut \(unique\):\s*(\d+)",
        'read_1': r"Read-1:\s*(\d+)",
        'read_2': r"Read-2:\s*(\d+)",
        'reads_map_to_sense': r"Reads map to '\+':\s*(\d+)",
        'reads_map_to_antisense': r"Reads map to '-':\s*(\d+)",
        'non-splice_reads': r"Non-splice reads:\s*(\d+)",
        'splice_reads': r"Splice reads:\s*(\d+)",
        'reads_mapped_in_proper_pairs': r"Reads mapped in proper pairs:\s*(\d+)",
        'proper-paired_reads_map_to_different_chrom': r"Proper-paired reads map to different chrom:\s*(\d+)",
    }

    #intiate PE check
    is_paired_end=False

    # Go through files and parse data using regexes
    for f in self.find_log_files('rseqc/bam_stat'):
        d = dict()
        for k, r in regexes.items():
            r_search = re.search(r, f['f'], re.MULTILINE)
            if r_search:
                d[k] = int(r_search.group(1))

        # Calculate some percentages
        if 'total_records' in d:
            t = float(d['total_records'])
            if 'mapq_gte_mapq_cut_unique' in d:
                d['unique_percent'] = (float(d['mapq_gte_mapq_cut_unique']) / t)*100.0
            if 'reads_mapped_in_proper_pairs' in d:
                d['proper_pairs_percent'] = (float(d['reads_mapped_in_proper_pairs']) / t)*100.0

        if len(d) > 0:
            if f['s_name'] in self.bam_stat_data:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f, section='bam_stat')
            #Check if SE or PE
            if d['read_2'] != 0:
                is_paired_end = True
            self.bam_stat_data[f['s_name']] = d

    # Filter to strip out ignored sample names
    self.bam_stat_data = self.ignore_samples(self.bam_stat_data)

    if len(self.bam_stat_data) > 0:
        # Write to file
        self.write_data_file(self.bam_stat_data, 'multiqc_rseqc_bam_stat')

        # Add to general stats table
        self.general_stats_headers['proper_pairs_percent'] = {
            'title': '% Proper Pairs',
            'description': '% Reads mapped in proper pairs',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'RdYlGn'
        }
        for s_name in self.bam_stat_data:
            if s_name not in self.general_stats_data:
                self.general_stats_data[s_name] = dict()

                #Only write if PE, i.e. there is something to write
                if is_paired_end:
                    self.general_stats_data[s_name].update( self.bam_stat_data[s_name] )

        # Make dot plot of counts
        pconfig = {
            'id': 'rseqc_bam_stat'
        }
        keys = OrderedDict()
        defaults = {
            'min': 0,
            'shared_key': 'read_count',
            'decimalPlaces': 2,
            'modify': lambda x: float(x) / 1000000.0,
        }
        keys['total_records'] = dict(defaults, **{'title': 'Total records' })
        keys['qc_failed'] = dict(defaults, **{'title': 'QC failed' })
        keys['optical_pcr_duplicate'] = dict(defaults, **{'title': 'Duplicates', 'description': 'Optical/PCR duplicate' })
        keys['non_primary_hits'] = dict(defaults, **{'title': 'Non primary hit' })
        keys['unmapped_reads'] = dict(defaults, **{'title': 'Unmapped', 'description': 'Unmapped reads' })
        keys['mapq_lt_mapq_cut_non'] = dict(defaults, **{'title': 'Non-unique', 'description': 'mapq < mapq_cut (non-unique)' })
        keys['mapq_gte_mapq_cut_unique'] = dict(defaults, **{'title': 'Unique', 'description': 'mapq >= mapq_cut (unique)' })
        if is_paired_end:
            keys['read_1'] = dict(defaults, **{'title': 'Read-1' })
            keys['read_2'] = dict(defaults, **{'title': 'Read-2' })
        keys['reads_map_to_sense'] = dict(defaults, **{'title': '+ve strand', 'description': "Reads map to '+'" })
        keys['reads_map_to_antisense'] = dict(defaults, **{'title': '-ve strand', 'description': "Reads map to '-'" })
        keys['non-splice_reads'] = dict(defaults, **{'title': 'Non-splice reads' })
        keys['splice_reads'] = dict(defaults, **{'title': 'Splice reads' })
        if is_paired_end:
            keys['reads_mapped_in_proper_pairs'] = dict(defaults, **{
                        'title': 'Proper pairs', 'description':'Reads mapped in proper pairs' })
            keys['proper-paired_reads_map_to_different_chrom'] = dict(defaults, **{
                        'title': 'Different chrom', 'description': 'Proper-paired reads map to different chrom' })

        self.add_section (
            name = 'Bam Stat',
            anchor = 'rseqc-bam_stat',
            description = 'All numbers reported in millions.',
            plot = beeswarm.plot(self.bam_stat_data, keys, pconfig)
        )

    # Return number of samples found
    return len(self.bam_stat_data)
845983
def parse_seqres(self, pdb):
    """"""
    
    seqresre = re.compile("SEQRES")
    
    seqreslines = [line for line in pdb.lines if seqresre.match(line)]
    
    for line in seqreslines:
      chain = line[11]
      resnames = line[19:70].strip()
      self.setdefault(chain, [])
      self[chain] += resnames.split()
159065
def _correct_qualimap_genome_results(samples):
    
    for s in samples:
        if verify_file(s.qualimap_genome_results_fpath):
            correction_is_needed = False
            with open(s.qualimap_genome_results_fpath, 'r') as f:
                content = f.readlines()
                metrics_started = False
                for line in content:
                    if ">> Reference" in line:
                        metrics_started = True
                    if metrics_started:
                        if line.find(',') != -1:
                            correction_is_needed = True
                            break
            if correction_is_needed:
                with open(s.qualimap_genome_results_fpath, 'w') as f:
                    metrics_started = False
                    for line in content:
                        if ">> Reference" in line:
                            metrics_started = True
                        if metrics_started:
                            if line.find(',') != -1:
                                line = line.replace(',', '')
                        f.write(line)
916854
def _parse_reports_by_type(self):
    

    data = dict()

    for file_meta in self.find_log_files('picard/sam_file_validation', filehandles=True):
        sample = file_meta['s_name']

        if sample in data:
            log.debug("Duplicate sample name found! Overwriting: {}".format(sample))

        filehandle = file_meta['f']
        first_line = filehandle.readline().rstrip()
        filehandle.seek(0)  # Rewind reading of the file

        if 'No errors found' in first_line:
            sample_data = _parse_no_error_report()
        elif first_line.startswith('ERROR') or first_line.startswith('WARNING'):
            sample_data = _parse_verbose_report(filehandle)
        else:
            sample_data = _parse_summary_report(filehandle)

        data[sample] = sample_data

    return data
372312
def _parse_qualimap_coverage(table):
    
    out = {}
    for row in table.find_all("tr"):
        col, val = [x.text for x in row.find_all("td")]
        if col == "Mean":
            out["Coverage (Mean)"] = val
    return out
990086
def parse_itms_services(url_data):
    """"""
    query = url_data.urlparts[3]
    for k, v, sep in urlutil.parse_qsl(query, keep_blank_values=True):
        if k == "url":
            url_data.add_url(v)
            break
799979
def parse_qc(self, qc_file):
        
        import pandas as pd
        series = pd.Series()
        try:
            with open(qc_file) as handle:
                line = handle.readlines()[0].strip().split("\t")  # list of strings per line
            series["NSC"] = line[-3]
            series["RSC"] = line[-2]
            series["qualityTag"] = line[-1]
        except:
            pass
        return series
388769
def parse_star_report (self, raw_data):
        

        regexes = {
            'total_reads':                  r"Number of input reads \|\s+(\d+)",
            'avg_input_read_length':        r"Average input read length \|\s+([\d\.]+)",
            'uniquely_mapped':              r"Uniquely mapped reads number \|\s+(\d+)",
            'uniquely_mapped_percent':      r"Uniquely mapped reads % \|\s+([\d\.]+)",
            'avg_mapped_read_length':       r"Average mapped length \|\s+([\d\.]+)",
            'num_splices':                  r"Number of splices: Total \|\s+(\d+)",
            'num_annotated_splices':        r"Number of splices: Annotated \(sjdb\) \|\s+(\d+)",
            'num_GTAG_splices':             r"Number of splices: GT/AG \|\s+(\d+)",
            'num_GCAG_splices':             r"Number of splices: GC/AG \|\s+(\d+)",
            'num_ATAC_splices':             r"Number of splices: AT/AC \|\s+(\d+)",
            'num_noncanonical_splices':     r"Number of splices: Non-canonical \|\s+(\d+)",
            'mismatch_rate':                r"Mismatch rate per base, % \|\s+([\d\.]+)",
            'deletion_rate':                r"Deletion rate per base \|\s+([\d\.]+)",
            'deletion_length':              r"Deletion average length \|\s+([\d\.]+)",
            'insertion_rate':               r"Insertion rate per base \|\s+([\d\.]+)",
            'insertion_length':             r"Insertion average length \|\s+([\d\.]+)",
            'multimapped':                  r"Number of reads mapped to multiple loci \|\s+(\d+)",
            'multimapped_percent':          r"% of reads mapped to multiple loci \|\s+([\d\.]+)",
            'multimapped_toomany':          r"Number of reads mapped to too many loci \|\s+(\d+)",
            'multimapped_toomany_percent':  r"% of reads mapped to too many loci \|\s+([\d\.]+)",
            'unmapped_mismatches_percent':  r"% of reads unmapped: too many mismatches \|\s+([\d\.]+)",
            'unmapped_tooshort_percent':    r"% of reads unmapped: too short \|\s+([\d\.]+)",
            'unmapped_other_percent':       r"% of reads unmapped: other \|\s+([\d\.]+)",
        }
        parsed_data = {}
        for k, r in regexes.items():
            r_search = re.search(r, raw_data, re.MULTILINE)
            if r_search:
                parsed_data[k] = float(r_search.group(1))
        # Figure out the numbers for unmapped as for some reason only the percentages are given
        try:
            total_mapped = parsed_data['uniquely_mapped'] + parsed_data['multimapped'] + parsed_data['multimapped_toomany']
            unmapped_count = parsed_data['total_reads'] - total_mapped
            total_unmapped_percent = parsed_data['unmapped_mismatches_percent'] + parsed_data['unmapped_tooshort_percent'] + parsed_data['unmapped_other_percent']
            try:
                parsed_data['unmapped_mismatches'] = int(round(unmapped_count * (parsed_data['unmapped_mismatches_percent'] / total_unmapped_percent), 0))
                parsed_data['unmapped_tooshort'] = int(round(unmapped_count * (parsed_data['unmapped_tooshort_percent'] / total_unmapped_percent), 0))
                parsed_data['unmapped_other'] = int(round(unmapped_count * (parsed_data['unmapped_other_percent'] / total_unmapped_percent), 0))
            except ZeroDivisionError:
                parsed_data['unmapped_mismatches'] = 0
                parsed_data['unmapped_tooshort'] = 0
                parsed_data['unmapped_other'] = 0
        except KeyError:
            pass

        if len(parsed_data) == 0: return None
        return parsed_data
1221617
def reporter(self):
        
        # Create a folder in which amplicons and raw BLAST results are to be stored
        detailed_reports = os.path.join(self.path, 'detailed_reports')
        make_path(detailed_reports)
        # Create the report path if necessary
        make_path(self.reportpath)
        with open(self.report, 'w') as report:
            # Initialise the header
            data = 'Sample,Gene,GenomeLocation,AmpliconSize,Contig,ForwardPrimers,ReversePrimers,' \
                   'ForwardMismatches,ReverseMismatches\n'
            for sample in self.metadata:
                # Initialise variables to convert attributes from set to list
                sample[self.analysistype].genes = dict()
                sample[self.analysistype].ntrange = dict()
                try:
                    # Check to ensure that there were amplicons present
                    if sample[self.analysistype].range:
                        # Iterate through each contig with genes present, and the list of genes on the contig
                        for contig, genes in sorted(sample[self.analysistype].genespresent.items()):
                            sample[self.analysistype].genes[contig] = list(genes)
                            # Iterate through the set of genes
                            for gene in sorted(genes):
                                # Determine which primers had the best lowest number of mismatches
                                # Create variables to store primer names, and number of mismatches
                                forward = list()
                                reverse = list()
                                forwardmismatches = int()
                                reversemismatches = int()
                                # Sort the dictionary of mismatches based on the number of mismatches
                                sorteddict = sorted(sample[self.analysistype].mismatches[contig][gene].items(),
                                                    key=operator.itemgetter(1))
                                # For every primer in the sorted dictionary, check to see if it is a forward or
                                # reverse primer, and determine the number of mismatches; if it is the first primer
                                # encountered for a gene, add it to the list, as because the dictionary is sorted,
                                # the number of mismatches should be the best. Additionally, if subsequent primer
                                # hits have the same number of mismatches, also add them to the list
                                for entry in sorteddict:
                                    # Add primer to the list if it is the first primer encountered
                                    if '-F' in entry[0] and not forward:
                                        forward.append(entry[0])
                                        forwardmismatches = entry[1]
                                    # Add the primer to the list if it has the same number of mismatches as a previous
                                    # primer in the list
                                    elif '-F' in entry[0] and entry[1] == forwardmismatches:
                                        forward.append(entry[0])
                                        forwardmismatches = entry[1]
                                    elif '-R' in entry[0] and not reverse:
                                        reverse.append(entry[0])
                                        reversemismatches = entry[1]
                                    elif '-R' in entry[0] and entry[1] == reversemismatches:
                                        reverse.append(entry[0])
                                        reversemismatches = entry[1]
                                # Make a variable to prevent writing out this long attribute name multiple times
                                ntrange = list(sample[self.analysistype].range[contig][gene])
                                sample[self.analysistype].ntrange[gene] = ntrange
                                # Extract the amplicons from the sequence file
                                self.ampliconfile(sample, contig, sorted(ntrange), forward, reverse)
                                # Copy the amplicons and raw BLAST outputs from FASTQ-formatted files to the
                                # detailed_reports folder
                                if sample[self.analysistype].filetype == 'fastq':
                                    try:
                                        shutil.copyfile(src=sample[self.analysistype].ampliconfile,
                                                        dst=os.path.join(detailed_reports,
                                                                         '{sn}_amplicons.fa'.format(sn=sample.name)))
                                        shutil.copyfile(src=sample[self.analysistype].report,
                                                        dst=os.path.join(detailed_reports,
                                                                         '{sn}_rawresults.csv'.format(sn=sample.name)))
                                    except FileNotFoundError:
                                        pass
                                # This first gene for a sample gets the sample name printed
                                data += '{sn},'.format(sn=sample.name)
                                # Populate the string with the gene name, properly formatted range, the length of
                                # the amplicon, and the name of the contig on which the gene was found
                                data += '{},{},{},{},{},{},{},{}\n'\
                                    .format(gene,
                                            '-'.join(str(x) for x in sorted(ntrange)),
                                            max(ntrange) - min(ntrange) + 1,
                                            contig,
                                            ';'.join(sorted(forward)),
                                            ';'.join(sorted(reverse)),
                                            forwardmismatches,
                                            reversemismatches)
                    # If there were no amplicons, add the sample name and nothing else
                    else:
                        data += '{sn}\n'.format(sn=sample.name)
                # If there were no BLAST hits, add the sample name, and nothing else
                except AttributeError:
                    data += '{sn}\n'.format(sn=sample.name)
                # Remove attributes that either take up too much room in the .json output, or are not JSON serializable
                try:
                    delattr(sample[self.analysistype], "blastresults")
                except AttributeError:
                    pass
                try:
                    delattr(sample[self.analysistype], "genespresent")
                except AttributeError:
                    pass
                try:
                    delattr(sample[self.analysistype], "contigs")
                except AttributeError:
                    pass
                try:
                    delattr(sample[self.analysistype], "range")
                except AttributeError:
                    pass
            # Write the string to the report
            report.write(data)
        # Clean up the BLAST database files
        db = os.path.splitext(self.formattedprimers)[0]
        # A list of all the file extensions associated with the BLASTdb
        dbextensions = ['.nhr', '.nin', '.nog', '.nsd', '.nsi', '.nsq']
        # Iterate through all the files, and delete each one - pass on IO errors
        for dbfile in zip(itertools.repeat(db), dbextensions):
            try:
                os.remove(''.join(dbfile))
            except IOError:
                pass
        try:
            os.remove(self.formattedprimers)  # Maybe want to keep this file?
        except IOError:
            pass
4955
def copies(mapping, s2bins, rna, min_rna = 800, mismatches = 0):
    
    cov = {} # cov[scaffold] = [bases, length]
    s2bins, bins2s = parse_s2bins(s2bins)
    rna_cov = parse_rna(rna, s2bins, min_rna)
    s2bins, bins2s = filter_missing_rna(s2bins, bins2s, rna_cov)
    # count bases mapped to scaffolds and rRNA gene regions
    for line in mapping:
        line = line.strip().split()
        # get scaffold lengths
        if line[0].startswith('@'):
            if line[0].startswith('@SQ') is False:
                continue
            s = line[1].split(':')[1]
            l = int(line[2].split(':')[1])
            # check if scaffold is binned
            if s not in s2bins:
                continue
            if s not in cov:
                cov[s] = [0, l]
        # check mismatch threshold
        mm = count_mismatches(line)
        if mm is False or mm > mismatches:
            continue
        # check that scaffold is in bin
        s, bases = line[2], len(line[9])
        if s not in cov:
            continue
        cov[s][0] += bases 
        rna_cov = rna_bases(rna_cov, s, bases, line) 
    print('# mismatches threshold: %s' % (mismatches))
    header = ['#rRNA scaffold', 'rRNA genes >=%sbp on scaffold' % (min_rna), \
            'rRNA coverage', \
            'bin', 'bin info', 'bin coverage', \
            'rRNAs >=%sbp in bin' % (min_rna), \
            'rRNA coverage/bin coverage', \
            'estimated number of copies']
    print('\t'.join(header))
    for bin, scaffolds in list(bins2s.items()):
        rna_count = sum([len(rna_cov[s][2]) for s in scaffolds if s in rna_cov])
        for s in scaffolds:
            if s not in rna_cov:
                continue
            out = []
            counts = rna_cov[s]
            bin_cov = calc_bin_cov(bins2s[bin], cov)
            num_genes = len(counts[2])
            rna_coverage = float(float(counts[0])/float(counts[1]))
            if bin_cov == 0:
                rna_div_bin = 0
            else:
                rna_div_bin = float(rna_coverage/bin_cov)
            est = int(max([rna_count, counts, rna_div_bin]))
            out = [s, num_genes, rna_coverage, bin, bin_cov, rna_count, rna_div_bin, est]
            print('\t'.join([str(i) for i in out]))