Find RSeQC gene_body_coverage reports and parse their data
172530
def parse_reports(self):
    

    # Set up vars
    self.gene_body_cov_hist_counts = dict()
    self.gene_body_cov_hist_percent = dict()

    # TODO - Do separate parsing step to find skewness values
    # and add these to the general stats table?

    # Go through files and parse data
    for f in self.find_log_files('rseqc/gene_body_coverage'):

        # RSeQC >= v2.4
        if f['f'].startswith('Percentile'):
            keys = []
            nrows = 0
            for l in f['f'].splitlines():
                s = l.split()
                if len(keys) == 0:
                    keys = s[1:]
                else:
                    nrows += 1
                    s_name = self.clean_s_name(s[0], f['root'])
                    if s_name in self.gene_body_cov_hist_counts:
                        log.debug("Duplicate sample name found! Overwriting: {}".format(s_name))
                    self.add_data_source(f, s_name, section='gene_body_coverage')
                    self.gene_body_cov_hist_counts[s_name] = OrderedDict()
                    for k, var in enumerate(s[1:]):
                        self.gene_body_cov_hist_counts[s_name][int(keys[k])] = float(var)
            if nrows == 0:
                log.warning("Empty geneBodyCoverage file found: {}".format(f['fn']))


        # RSeQC < v2.4
        elif f['f'].startswith('Total reads'):
            if f['s_name'].endswith('.geneBodyCoverage'):
                f['s_name'] = f['s_name'][:-17]
            if f['s_name'] in self.gene_body_cov_hist_counts:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f, section='gene_body_coverage')
            self.gene_body_cov_hist_counts[f['s_name']] = OrderedDict()
            nrows = 0
            for l in f['f'].splitlines():
                s = l.split()
                try:
                    nrows += 1
                    self.gene_body_cov_hist_counts[f['s_name']][int(s[0])] = float(s[1])
                except:
                    pass
            if nrows == 0:
                del self.gene_body_cov_hist_counts[f['s_name']]
                log.warning("Empty geneBodyCoverage file found: {}".format(f['fn']))

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

    if len(self.gene_body_cov_hist_counts) > 0:

        # Make a normalised percentage version of the data
        for s_name in self.gene_body_cov_hist_counts:
            self.gene_body_cov_hist_percent[s_name] = OrderedDict()
            total = sum( self.gene_body_cov_hist_counts[s_name].values() )
            for k, v in self.gene_body_cov_hist_counts[s_name].items():
                self.gene_body_cov_hist_percent[s_name][k] = (v/total)*100

        # Add line graph to section
        pconfig = {
            'id': 'rseqc_gene_body_coverage_plot',
            'title': 'RSeQC: Gene Body Coverage',
            'ylab': '% Coverage',
            'xlab': "Gene Body Percentile (5' -> 3')",
            'xmin': 0,
            'xmax': 100,
            'tt_label': "<strong>{point.x}% from 5'</strong>: {point.y:.2f}",
            'data_labels': [
                {'name': 'Percentages', 'ylab': 'Percentage Coverage'},
                {'name': 'Counts', 'ylab': 'Coverage'}
            ]
        }
        self.add_section (
            name = 'Gene Body Coverage',
            anchor = 'rseqc-gene_body_coverage',
            description = '<a href="http://rseqc.sourceforge.net/#genebody-coverage-py" target="_blank">Gene Body Coverage</a>' \
                " calculates read coverage over gene bodies." \
                " This is used to check if reads coverage is uniform and" \
                " if there is any 5' or 3' bias.",
            plot = linegraph.plot([self.gene_body_cov_hist_percent, self.gene_body_cov_hist_counts], pconfig)
        )

    # Return number of samples found
    return len(self.gene_body_cov_hist_counts)
2688
def parse_reports(self):
    

    # Set up vars
    self.junction_saturation_all = dict()
    self.junction_saturation_known = dict()
    self.junction_saturation_novel = dict()

    # Go through files and parse data
    for f in self.find_log_files('rseqc/junction_saturation'):
        parsed = dict()
        for l in f['f'].splitlines():
            r = re.search(r"^([xyzw])=c\(([\d,]+)\)$", l)
            if r:
                parsed[r.group(1)] = [float(i) for i in r.group(2).split(',')]
        if len(parsed) == 4:
            if parsed['z'][-1] == 0:
                log.warn("Junction saturation data all zeroes, skipping: '{}'".format(f['s_name']))
            else:
                if f['s_name'] in self.junction_saturation_all:
                    log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
                self.add_data_source(f, section='junction_saturation')
                self.junction_saturation_all[f['s_name']] = OrderedDict()
                self.junction_saturation_known[f['s_name']] = OrderedDict()
                self.junction_saturation_novel[f['s_name']] = OrderedDict()
                for k, v in enumerate(parsed['x']):
                    self.junction_saturation_all[f['s_name']][v] = parsed['z'][k]
                    self.junction_saturation_known[f['s_name']][v] = parsed['y'][k]
                    self.junction_saturation_novel[f['s_name']][v] = parsed['w'][k]

    # Filter to strip out ignored sample names
    self.junction_saturation_all = self.ignore_samples(self.junction_saturation_all)
    self.junction_saturation_known = self.ignore_samples(self.junction_saturation_known)
    self.junction_saturation_novel = self.ignore_samples(self.junction_saturation_novel)

    if len(self.junction_saturation_all) > 0:

        # Add line graph to section
        pconfig = {
            'id': 'rseqc_junction_saturation_plot',
            'title': 'RSeQC: Junction Saturation',
            'ylab': 'Number of Junctions',
            'ymin': 0,
            'xlab': "Percent of reads",
            'xmin': 0,
            'xmax': 100,
            'tt_label': "<strong>{point.x}% of reads</strong>: {point.y:.2f}",
            'data_labels': [
                {'name': 'Known Junctions'},
                {'name': 'Novel Junctions'},
                {'name': 'All Junctions'}
            ],
            'cursor': 'pointer',
            'click_func': plot_single()
        }
        self.add_section (
            name = 'Junction Saturation',
            anchor = 'rseqc-junction_saturation',
            description = '''<a href="http://rseqc.sourceforge.net/#junction-saturation-py" target="_blank">Junction Saturation</a>
                counts the number of known splicing junctions that are observed
                in each dataset. If sequencing depth is sufficient, all (annotated) splice junctions should
                be rediscovered, resulting in a curve that reaches a plateau. Missing low abundance splice
                junctions can affect downstream analysis.</p>
                <div class="alert alert-info" id="rseqc-junction_sat_single_hint">
                  <span class="glyphicon glyphicon-hand-up"></span>
                  Click a line to see the data side by side (as in the original RSeQC plot).
                </div><p>''',
            plot = linegraph.plot([
                    self.junction_saturation_known,
                    self.junction_saturation_novel,
                    self.junction_saturation_all
                ], pconfig)
        )

    # Return number of samples found
    return len(self.junction_saturation_all)
78156
def parse_epcr(self):
        
        # Use the metadata object from the vtyper_object
        for sample in self.vtyper_object.metadata:
            # Initialise the dictionary
            sample[self.analysistype].result_dict = dict()
            # Read in the output file
            with open(sample[self.analysistype].resultsfile) as epcrresults:
                for result in epcrresults:
                    # Only the lines without a # contain results
                    if "#" not in result:
                        # Split on \t
                        # vtx2a_0_0 2014-SEQ-0121_127_length_1407_cov_50.7797_ID_10924  -   228   576 2   0 349/100-350
                        # primer_set: vtx2a_0_0, contig: 2014-SEQ-0121_127_length_1407_cov_50.7797_ID_10924, strand: -,
                        # start: 228, stop: 576, number of forward mismatches: 2, number of reverse mismatches: 2
                        # amplicon_combo: 349/100-350
                        primer_set, contig, strand, start, stop, total_mismatches, indels, amplicon_combo = \
                            result.rstrip().split('\t')
                        # Set the mismatches to be an int
                        total_mismatches = int(total_mismatches)
                        # Set the position of the amplicon on the contig. Ensure that the lower value is first
                        genome_pos = '{min}-{max}'.format(min=min([int(start), int(stop)]),
                                                          max=max([int(start), int(stop)]))
                        # Extract the gene name from the modified name used when creating the primer file: LMhlyA_0_0
                        # becomes LMhlyA
                        gene_re = re.search(r'([\w-]+)_(\d{1,3})_(\d{1,3})', primer_set)
                        gene = gene_re.groups()[0]
                        # Split the amplicon length from amplicon_combo: 349/100-350 -> 349
                        amplicon_length = amplicon_combo.split('/')[0]
                        # Populate the dictionary if the 'total_mismatches' key doesn't exist, or if the current number
                        # of mismatches is better than the previous 'best' number of mismatches
                        try:
                            if total_mismatches < sample[self.analysistype].result_dict[gene]['total_mismatches']:
                                self.populate_results_dict(sample=sample,
                                                           gene=gene,
                                                           total_mismatches=total_mismatches,
                                                           genome_pos=genome_pos,
                                                           amplicon_length=amplicon_length,
                                                           contig=contig,
                                                           primer_set=primer_set)
                        except KeyError:
                            self.populate_results_dict(sample=sample,
                                                       gene=gene,
                                                       total_mismatches=total_mismatches,
                                                       genome_pos=genome_pos,
                                                       amplicon_length=amplicon_length,
                                                       contig=contig,
                                                       primer_set=primer_set)
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]))
1109604
def _whole_gene_covered_by_nucmer_hits(nucmer_hits, ref_seq, percent_threshold, max_nt_extend):
        ''''''
        coords = AssemblyCompare.nucmer_hits_to_ref_coords(nucmer_hits)
        covered = []
        for coords_list in coords.values():
            covered.extend(coords_list)
        pyfastaq.intervals.merge_overlapping_in_list(covered)
        return (2 * max_nt_extend + pyfastaq.intervals.length_sum_from_list(covered)) / len(ref_seq) >= percent_threshold
325870
def parse_fastqc_report(self, file_contents, s_name=None, f=None):
        

        # Make the sample name from the input filename if we find it
        fn_search = re.search(r"Filename\s+(.+)", file_contents)
        if fn_search:
            s_name = self.clean_s_name(fn_search.group(1) , f['root'])

        if s_name in self.fastqc_data.keys():
            log.debug("Duplicate sample name found! Overwriting: {}".format(s_name))
        self.add_data_source(f, s_name)
        self.fastqc_data[s_name] = { 'statuses': dict() }

        # Parse the report
        section = None
        s_headers = None
        self.dup_keys = []
        for l in file_contents.splitlines():
            if l == '>>END_MODULE':
                section = None
                s_headers = None
            elif l.startswith('>>'):
                (section, status) = l[2:].split("\t", 1)
                section = section.lower().replace(' ', '_')
                self.fastqc_data[s_name]['statuses'][section] = status
            elif section is not None:
                if l.startswith('#'):
                    s_headers = l[1:].split("\t")
                    # Special case: Total Deduplicated Percentage header line
                    if s_headers[0] == 'Total Deduplicated Percentage':
                        self.fastqc_data[s_name]['basic_statistics'].append({
                            'measure': 'total_deduplicated_percentage',
                            'value': float(s_headers[1])
                        })
                    else:
                        # Special case: Rename dedup header in old versions of FastQC (v10)
                        if s_headers[1] == 'Relative count':
                            s_headers[1] = 'Percentage of total'
                        s_headers = [s.lower().replace(' ', '_') for s in s_headers]
                        self.fastqc_data[s_name][section] = list()

                elif s_headers is not None:
                    s = l.split("\t")
                    row = dict()
                    for (i, v) in enumerate(s):
                        v.replace('NaN','0')
                        try:
                            v = float(v)
                        except ValueError:
                            pass
                        row[s_headers[i]] = v
                    self.fastqc_data[s_name][section].append(row)
                    # Special case - need to remember order of duplication keys
                    if section == 'sequence_duplication_levels':
                        try:
                            self.dup_keys.append(float(s[0]))
                        except ValueError:
                            self.dup_keys.append(s[0])

        # Tidy up the Basic Stats
        self.fastqc_data[s_name]['basic_statistics'] = {d['measure']: d['value'] for d in self.fastqc_data[s_name]['basic_statistics']}

        # Calculate the average sequence length (Basic Statistics gives a range)
        length_bp = 0
        total_count = 0
        for d in self.fastqc_data[s_name].get('sequence_length_distribution', {}):
            length_bp += d['count'] * self.avg_bp_from_range(d['length'])
            total_count += d['count']
        if total_count > 0:
            self.fastqc_data[s_name]['basic_statistics']['avg_sequence_length'] = length_bp / total_count
368244
def make_scrnaseq_object(samples):
    
    local_sitelib = R_sitelib()
    counts_dir = os.path.dirname(dd.get_in_samples(samples, dd.get_combined_counts))
    gtf_file = dd.get_in_samples(samples, dd.get_transcriptome_gtf)
    if not gtf_file:
        gtf_file = dd.get_in_samples(samples, dd.get_gtf_file)
    rda_file = os.path.join(counts_dir, "se.rda")
    if not file_exists(rda_file):
        with file_transaction(rda_file) as tx_out_file:
            rcode = "%s-run.R" % os.path.splitext(rda_file)[0]
            rrna_file = "%s-rrna.txt" % os.path.splitext(rda_file)[0]
            rrna_file = _find_rRNA_genes(gtf_file, rrna_file)
            with open(rcode, "w") as out_handle:
                out_handle.write(_script.format(**locals()))
            rscript = Rscript_cmd()
            try:
                # do.run([rscript, "--no-environ", rcode],
                #        "SingleCellExperiment",
                #        log_error=False)
                rda_file = rcode
            except subprocess.CalledProcessError as msg:
                logger.exception()
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()
350590
def dump_rexobj_results(rexobj, options=None):
    '''
    
    '''
    print("-" * 60)
    print("Match count: ", rexobj.res_count)
    matches = rexobj.matches
    for match in matches:
        print("Loc:", match.loc, ":: ")
        for key in match.named_groups.keys():
            print("%s: %s" %
                  (key, match.named_groups[key]))
        print("")
107239
def parse_reports(self):
    

    # Set up vars
    self.bamtools_stats_data = dict()
    regexes = {
        'total_reads': r"Total reads:\s*(\d+)",
        'mapped_reads': r"Mapped reads:\s*(\d+)",
        'mapped_reads_pct': r"Mapped reads:\s*\d+\s+\(([\d\.]+)%\)",
        'forward_strand': r"Forward strand:\s*(\d+)",
        'forward_strand_pct': r"Forward strand:\s*\d+\s+\(([\d\.]+)%\)",
        'reverse_strand': r"Reverse strand:\s*(\d+)",
        'reverse_strand_pct': r"Reverse strand:\s*\d+\s+\(([\d\.]+)%\)",
        'failed_qc': r"Failed QC:\s*(\d+)",
        'failed_qc_pct': r"Failed QC:\s*\d+\s+\(([\d\.]+)%\)",
        'duplicates': r"Duplicates:\s*(\d+)",
        'duplicates_pct': r"Duplicates:\s*\d+\s+\(([\d\.]+)%\)",
        'paired_end': r"Paired-end reads:\s*(\d+)",
        'paired_end_pct': r"Paired-end reads:\s*\d+\s+\(([\d\.]+)%\)",
        'proper_pairs': r"'Proper-pairs'\s*(\d+)",
        'proper_pairs_pct': r"'Proper-pairs'\s*\d+\s+\(([\d\.]+)%\)",
        'both_mapped': r"Both pairs mapped:\s*(\d+)",
        'both_mapped_pct': r"Both pairs mapped:\s*\d+\s+\(([\d\.]+)%\)",
        'read_1': r"Read 1:\s*(\d+)",
        'read_2': r"Read 2:\s*(\d+)",
        'singletons': r"Singletons:\s*(\d+)",
        'singletons_pct': r"Singletons:\s*\d+\s+\(([\d\.]+)%\)",
    }

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

        if len(d) > 0:
            if f['s_name'] in self.bamtools_stats_data:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f, section='stats')
            self.bamtools_stats_data[f['s_name']] = d

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

    if len(self.bamtools_stats_data) > 0:

        # Write to file
        self.write_data_file(self.bamtools_stats_data, 'multiqc_bamtools_stats')

        # Add to general stats table
        self.general_stats_headers['duplicates_pct'] = {
            'title': '% Duplicates',
            'description': '% Duplicate Reads',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'OrRd'
        }
        self.general_stats_headers['mapped_reads_pct'] = {
            'title': '% Mapped',
            'description': '% Mapped Reads',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'RdYlGn'
        }
        for s_name in self.bamtools_stats_data:
            if s_name not in self.general_stats_data:
                self.general_stats_data[s_name] = dict()
            self.general_stats_data[s_name].update( self.bamtools_stats_data[s_name] )

        # Make dot plot of counts
        keys = OrderedDict()
        defaults = {
            'min': 0,
            'max': 100,
            'decimalPlaces': 2,
            'suffix': '%'
        }
        num_defaults = {
            'min': 0,
            'modify': lambda x: float(x) / 1000000.0,
            'decimalPlaces': 2
        }

        keys['total_reads'] = dict(num_defaults, **{'title': 'Total reads', 'description': 'Total reads (millions)' });
        keys['mapped_reads_pct'] = dict(defaults, **{'title': 'Mapped reads' })
        keys['forward_strand_pct'] = dict(defaults, **{'title': 'Forward strand' })
        keys['reverse_strand_pct'] = dict(defaults, **{'title': 'Reverse strand' })
        keys['failed_qc_pct'] = dict(defaults, **{'title': 'Failed QC' })
        keys['duplicates_pct'] = dict(defaults, **{'title': 'Duplicates' })
        keys['paired_end_pct'] = dict(defaults, **{'title': 'Paired-end', 'description': 'Paired-end reads' })
        keys['proper_pairs_pct'] = dict(defaults, **{'title': 'Proper-pairs' })
        keys['both_mapped_pct'] = dict(defaults, **{'title': 'Both mapped', 'description': 'Both pairs mapped' })
        keys['bt_read_1'] = dict(num_defaults, **{'title': 'Read 1', 'description': 'Read 1 (millions)' });
        keys['bt_read_2'] = dict(num_defaults, **{'title': 'Read 2', 'description': 'Read 2 (millions)' });
        keys['singletons_pct'] = dict(defaults, **{'title': 'Singletons' })

        self.add_section (
            name = 'Bamtools Stats',
            anchor = 'bamtools-stats',
            plot = beeswarm.plot(self.bamtools_stats_data, keys)
        )

    # Return number of samples found
    return len(self.bamtools_stats_data)
1057932
def print_table(seqs, id2name, name):
    
    itable = open('%s.itable' % (name.rsplit('.', 1)[0]), 'w')
    print('\t'.join(['#sequence', 'gene', 'model', 'insertion', 'gene position', 'model position', 'length', 'orf?', 'intron?', 'orf?intron?', 'insertion', 'orf', 'intron']), file=itable)
    for seq, info in list(seqs.items()):
        gene, model, insertions = info
        name = id2name[seq]
        for i, ins in enumerate(insertions, 1):
            gene_pos, model_pos, length, iseq, \
                    orfs, introns, orfs_b, introns_b, orf_annotations = ins
            # append annotation to orf header
            for orf in orfs:
                parts = orf[0].split()
                annotation = orf_annotations[parts[0].split('>')[1]]
                orf[0] = '%s %s %s' % (parts[0], annotation, ' '.join(parts[1:]))
            # get orf position
            gene_pos = '-'.join([str(j) for j in gene_pos])
            # check if orf, intron is present
            if orfs_b is True or introns_b is True:
                orfs_introns_b = True
            else:
                orfs_introns_b = False
            out = [name, gene, model, i, gene_pos, model_pos, length, orfs_b, introns_b, orfs_introns_b]
            out.append('|'.join(iseq))
            out.append('|'.join(['|'.join(orf) for orf in orfs]))
            out.append('|'.join(['|'.join(intron) for intron in introns]))
            print('\t'.join([str(i) for i in out]), file=itable)
    itable.close()
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())
923734
def parse_coverage(coverage_report, parser):
    
    if parser in PARSERS:
        if coverage_report:
            return PARSERS[parser].parse_coverage_report(coverage_report)
        return None
    return NotImplemented
1059720
def parse_coverage (self, f):
        
        data = dict()
        s_names = None
        j = 1
        for l in f['f'].splitlines():
            s = l.strip().split("\t")
            if s_names is None:
                s_names = s
                for s_name in s_names:
                    data[s_name] = dict()
            else:
                for i, v in enumerate(s):
                    data[s_names[i]][j] = float(v)
                j += 1
        if f['fn'] == 'meanCoverageNorm_high.txt':
            self.rna_seqc_norm_high_cov.update(data)
        elif f['fn'] == 'meanCoverageNorm_medium.txt':
            self.rna_seqc_norm_medium_cov.update(data)
        elif f['fn'] == 'meanCoverageNorm_low.txt':
            self.rna_seqc_norm_low_cov.update(data)
164845
def parse_genome_results(self, f):
    
    regexes = {
        'bam_file': r"bam file = (.+)",
        'total_reads': r"number of reads = ([\d,]+)",
        'mapped_reads': r"number of mapped reads = ([\d,]+)",
        'mapped_bases': r"number of mapped bases = ([\d,]+)",
        'sequenced_bases': r"number of sequenced bases = ([\d,]+)",
        'mean_insert_size': r"mean insert size = ([\d,\.]+)",
        'median_insert_size': r"median insert size = ([\d,\.]+)",
        'mean_mapping_quality': r"mean mapping quality = ([\d,\.]+)",
        'general_error_rate': r"general error rate = ([\d,\.]+)",
    }
    d = dict()
    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 ValueError:
                d[k] = r_search.group(1)
    # Check we have an input filename
    if 'bam_file' not in d:
        log.debug("Couldn't find an input filename in genome_results file {}".format(f['fn']))
        return None

    # Get a nice sample name
    s_name = self.clean_s_name(d['bam_file'], f['root'])

    # Add to general stats table & calculate a nice % aligned
    try:
        self.general_stats_data[s_name]['total_reads'] = d['total_reads']
        self.general_stats_data[s_name]['mapped_reads'] = d['mapped_reads']
        d['percentage_aligned'] = (d['mapped_reads'] / d['total_reads'])*100
        self.general_stats_data[s_name]['percentage_aligned'] = d['percentage_aligned']
        self.general_stats_data[s_name]['general_error_rate'] = d['general_error_rate']*100
    except KeyError:
        pass

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

    # get data
    data = collect_data(parent_module)

    # Filter to strip out ignored sample names
    data = parent_module.ignore_samples(data)

    # Reference data in parent module
    parent_module.picard_variantCalling_data = data

    if len(data) > 0:
        derive_data(data)

        # Write parsed data to a file
        parent_module.write_data_file(data, 'multiqc_picard_varientCalling')

        parent_module.general_stats_headers['DBSNP_TITV'] = {
            'title': 'TiTV ratio (known)',
            'description': "The Transition/Transversion ratio of the passing bi-allelic SNP calls made at SNP-database sites.",
            'min': 0,
            'scale': 'Blues',
            'shared_key': 'titv_ratio'
        }

        parent_module.general_stats_headers['NOVEL_TITV'] = {
            'title': 'TiTV ratio (novel)',
            'description': "The Transition/Transversion ratio of the passing bi-allelic SNP calls made at non-SNP-database sites.",
            'min': 0,
            'scale': 'Blues',
            'shared_key': 'titv_ratio'
        }

        parent_module.general_stats_headers['DBSNP_INS_DEL_RATIO'] = {
            'title': 'InDel ratio (known)',
            'description': "The Insertion / Deletion ratio of the passing bi-allelic SNP calls made at SNP-database sites.",
            'min': 0,
            'scale': 'Greens',
            'shared_key': 'indel_ratio',
            'hidden': True,
        }

        parent_module.general_stats_headers['NOVEL_INS_DEL_RATIO'] = {
            'title': 'InDel ratio (novel)',
            'description': "The Insertion / Deletion ratio of the passing bi-allelic SNP calls made at non-SNP-database sites.",
            'min': 0,
            'scale': 'Greens',
            'shared_key': 'indel_ratio',
            'hidden': True,
        }


        parent_module.general_stats_headers['total_called_variants_known'] = {
            'title': 'Called Variants (known)',
            'description': "Total counts of variants in SNP-database sites.",
            'shared_key': 'variant_count',
            'min': 0,
            'format': '{0:,.0f}',
            'hidden': True,
        }

        parent_module.general_stats_headers['total_called_variants_novel'] = {
            'title': 'Called Variants (novel)',
            'description': "Total counts of variants in non-SNP-database sites.",
            'shared_key': 'variant_count',
            'min': 0,
            'format': '{0:,.0f}',
            'hidden': True,
        }

        for s_name in data:
            if s_name not in parent_module.general_stats_data:
                parent_module.general_stats_data[s_name] = dict()
            parent_module.general_stats_data[s_name].update(data[s_name])


        # Variant Counts Bargraph
        parent_module.add_section(
            name='Variant Types',
            anchor='picard-variants-types',
            description='Variants that have been called, looking at variant types. Optionally filtered on label.',
            helptext="""
            Only passing variants are shown (i.e. non-filtered).\n
            SNPs are bi-allelic.\n
            Complex InDels are both an insertion and a deletion.
            """,
            plot=compare_variant_type_plot(data)
        )

        # Variant Counts Table
        parent_module.add_section(
            name='Variant Labels',
            anchor='picard-variants-labels',
            description='Variants that have been called, comparing with known variant sites.',
            helptext="""
            Only passing variants are shown (i.e. non-filtered).\n
            Variants contain bi-allelic SNPs, multi-allelic SNPs, simple and complex inserts and deletions.
            """,
            plot=compare_variants_label_plot(data)
        )

    return len(data)
471199
def get_protein_data_pgrouped(proteindata, p_acc, headerfields):
    """"""
    report = get_protein_data_base(proteindata, p_acc, headerfields)
    return get_cov_protnumbers(proteindata, p_acc, report)