Parse the metrics.tsv file from RNA-SeQC
764434
def parse_metrics(self, f):
        
        headers = None
        for l in f['f'].splitlines():
            s = l.strip().split("\t")
            if headers is None:
                headers = s
            else:
                s_name = s[ headers.index('Sample') ]
                data = dict()
                for idx, h in enumerate(headers):
                    try:
                        data[h] = float(s[idx])
                    except ValueError:
                        data[h] = s[idx]
                self.rna_seqc_metrics[s_name] = data
502417
def _parse_reactome_association_file(
            self, file, limit=None, subject_prefix=None, object_prefix=None):
        
        eco_map = Reactome.get_eco_map(Reactome.map_files['eco_map'])
        count = 0
        with open(file, 'r') as tsvfile:
            reader = csv.reader(tsvfile, delimiter="\t")
            for row in reader:
                (component, pathway_id, pathway_iri, pathway_label, go_ecode,
                species_name) = row
                count += 1
                self._add_component_pathway_association(
                    eco_map, component, subject_prefix, pathway_id, object_prefix,
                    pathway_label, go_ecode)

                if limit is not None and count >= limit:
                    break

        return
266347
def parse(handle, seqtype=None, robust=False):
    ''''''
    if seqtype is None:
        seqtype = _get_seqtype_from_ext(handle)

    if seqtype.startswith('gz-'):
        handle = _unzip_handle(handle)
        seqtype = seqtype[3:]

    # False positive from pylint, both handles are fileobj-like
    # pylint: disable=redefined-variable-type
    if robust:
        if seqtype == "embl":
            handle = sanity_check_embl(handle)
        elif seqtype == "genbank":
            handle = sanity_check_genbank(handle)
        elif seqtype == "fasta":
            handle = sanity_check_fasta(handle)
    # pylint: enable=redefined-variable-type

    return SeqIO.parse(handle, seqtype)
531990
def parse_reports(self):
    

    # Set up vars
    self.read_dist = dict()
    first_regexes = {
        'total_reads': r"Total Reads\s+(\d+)\s*",
        'total_tags': r"Total Tags\s+(\d+)\s*",
        'total_assigned_tags': r"Total Assigned Tags\s+(\d+)\s*",
    }
    second_regexes = {
        'cds_exons': r"CDS_Exons\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        '5_utr_exons': r"5'UTR_Exons\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        '3_utr_exons': r"3'UTR_Exons\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'introns': r"Introns\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tss_up_1kb': r"TSS_up_1kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tss_up_5kb': r"TSS_up_5kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tss_up_10kb': r"TSS_up_10kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tes_down_1kb': r"TES_down_1kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tes_down_5kb': r"TES_down_5kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
        'tes_down_10kb': r"TES_down_10kb\s+(\d+)\s+(\d+)\s+([\d\.]+)\s*",
    }

    # Go through files and parse data using regexes
    for f in self.find_log_files('rseqc/read_distribution'):
        d = dict()
        for k, r in first_regexes.items():
            r_search = re.search(r, f['f'], re.MULTILINE)
            if r_search:
                d[k] = int(r_search.group(1))
        for k, r in second_regexes.items():
            r_search = re.search(r, f['f'], re.MULTILINE)
            if r_search:
                d['{}_total_bases'.format(k)] = int(r_search.group(1))
                d['{}_tag_count'.format(k)] = int(r_search.group(2))
                d['{}_tags_kb'.format(k)] = float(r_search.group(2))

        d['other_intergenic_tag_count'] = d['total_tags']-d['total_assigned_tags']

        # Calculate some percentages for parsed file
        if 'total_tags' in d:
            t = float(d['total_tags'])
            pcts = dict()
            for k in d:
                if k.endswith('_tag_count'):
                    pk = '{}_tag_pct'.format(k[:-10])
                    pcts[pk] = (float(d[k]) / t)*100.0
            d.update(pcts)

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

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

    if len(self.read_dist) > 0:

        # Write to file
        self.write_data_file(self.read_dist, 'multiqc_rseqc_read_distribution')

        # Plot bar graph of groups
        keys = OrderedDict()
        keys['cds_exons_tag_count'] = {'name': "CDS_Exons"}
        keys['5_utr_exons_tag_count'] = {'name': "5'UTR_Exons"}
        keys['3_utr_exons_tag_count'] = {'name': "3'UTR_Exons"}
        keys['introns_tag_count'] = {'name': "Introns"}
        keys['tss_up_1kb_tag_count'] = {'name': "TSS_up_1kb"}
        keys['tss_up_5kb_tag_count'] = {'name': "TSS_up_5kb"}
        keys['tss_up_10kb_tag_count'] = {'name': "TSS_up_10kb"}
        keys['tes_down_1kb_tag_count'] = {'name': "TES_down_1kb"}
        keys['tes_down_5kb_tag_count'] = {'name': "TES_down_5kb"}
        keys['tes_down_10kb_tag_count'] = {'name': "TES_down_10kb"}
        keys['other_intergenic_tag_count'] = {'name': "Other_intergenic"}

        # Config for the plot
        pconfig = {
            'id': 'rseqc_read_distribution_plot',
            'title': 'RSeQC: Read Distribution',
            'ylab': '# Tags',
            'cpswitch_counts_label': 'Number of Tags',
            'cpswitch_c_active': False
        }

        self.add_section (
            name = 'Read Distribution',
            anchor = 'rseqc-read_distribution',
            description = '<a href="http://rseqc.sourceforge.net/#read-distribution-py" target="_blank">Read Distribution</a>' \
                " calculates how mapped reads are distributed over genome features.",
            plot = bargraph.plot(self.read_dist, keys, pconfig)
        )

    # Return number of samples found
    return len(self.read_dist)
894895
def _parse_novel(csv_file, sps="new"):
    """"""
    read = 0
    seen = set()
    safe_makedir("novel")
    with open("novel/hairpin.fa", "w") as fa_handle, open("novel/miRNA.str", "w") as str_handle:
        with open(csv_file) as in_handle:
            for line in in_handle:
                if line.startswith("mature miRBase miRNAs detected by miRDeep2"):
                    break
                if line.startswith("novel miRNAs predicted"):
                    read = 1
                    line = next(in_handle)
                    continue
                if read and line.strip():
                    cols = line.strip().split("\t")
                    name, start, score = cols[0], cols[16], cols[1]
                    if float(score) < 1:
                        continue
                    m5p, m3p, pre = cols[13], cols[14], cols[15].replace('u', 't').upper()
                    m5p_start = cols[15].find(m5p) + 1
                    m3p_start = cols[15].find(m3p) + 1
                    m5p_end = m5p_start + len(m5p) - 1
                    m3p_end = m3p_start + len(m3p) - 1
                    if m5p in seen:
                        continue
                    fa_handle.write(">{sps}-{name} {start}\n{pre}\n".format(**locals()))
                    str_handle.write(">{sps}-{name} ({score}) [{sps}-{name}-5p:{m5p_start}-{m5p_end}] [{sps}-{name}-3p:{m3p_start}-{m3p_end}]\n".format(**locals()))
                    seen.add(m5p)
    return op.abspath("novel")
1182586
def parse_result_line(cls, columns, line, metrics):
        
        line = IORMetricsExtractor.RE_MULTIPLE_SPACES.sub(' ', line)
        line = line.split(' ')
        operation = line[0]
        assert len(line) == len(columns)
        for i, col in enumerate(line):
            col = columns[i]
            if col in cls.METAS_IGNORED:
                continue
            desc = cls.METAS.get(col)
            if desc is None:
                LOGGER.warn('Unrecognized column: %s', col)
            meta = cls.get_meta_name(operation, desc.get('name') or col)
            value = desc['metric'].type(line[i])
            metrics[meta] = value
1207364
def _read_gtf(gtf):
    
    if not gtf:
        return gtf
    db = defaultdict(list)
    with open(gtf) as in_handle:
        for line in in_handle:
            if line.startswith("#"):
                continue
            cols = line.strip().split("\t")
            name = [n.split("=")[1] for n in cols[-1].split(";") if n.startswith("Name")]
            chrom, start, end, strand = cols[0], cols[3], cols[4], cols[6]
            if cols[2] == "miRNA_primary_transcript":
                db[name[0]].append([chrom, int(start), int(end), strand])
    return db
846635
def get_rna_counts(rna_file, transcript_name):
    
    # check if the RNA file exists
    if not os.path.exists(rna_file):
        msg = 'RNA-Seq BAM file "{}" does not exist'.format(rna_file)
        logging.error(msg)
        raise OSError(msg)
    rna_counts = {}

    cov_file = tempfile.NamedTemporaryFile(delete=False)
    try:
        subprocess.check_call(
            ['bedtools', 'genomecov', '-ibam', rna_file,
             '-bg'], stdout=cov_file)
    except subprocess.CalledProcessError as e:
        # needs testing
        raise ribocore.RNACountsError('Could not generate coverage for RNA BAM file. \n{}\n'.format(e))
    for line in open(cov_file.name):
        line = line.split()
        if line[0] == transcript_name:
            position, count = int(line[1]) + 1, int(line[3])
            rna_counts[position] = count
    cov_file.close()
    os.unlink(cov_file.name)
    return rna_counts
1362832
def parse_htseq_report (self, f):
        
        keys = [ '__no_feature', '__ambiguous', '__too_low_aQual', '__not_aligned', '__alignment_not_unique' ]
        parsed_data = dict()
        assigned_counts = 0
        for l in f['f']:
            s = l.split("\t")
            if s[0] in keys:
                parsed_data[s[0][2:]] = int(s[-1])
            else:
                try:
                    assigned_counts += int(s[-1])
                except (ValueError, IndexError):
                    pass
        if len(parsed_data) > 0:
            parsed_data['assigned'] = assigned_counts
            parsed_data['total_count'] = sum([v for v in parsed_data.values()])
            parsed_data['percent_assigned'] = (float(parsed_data['assigned']) / float(parsed_data['total_count'])) * 100.0
            return parsed_data
        return None
265222
def summary(args):
    
    from jcvi.formats.base import DictFile
    from jcvi.utils.cbook import percentage, Registry

    p = OptionParser(summary.__doc__)
    p.add_option("--extra", help="Cross with extra tsv file [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    frfile, statusfile = args
    status = DictFile(statusfile)
    fp = open(frfile)
    registry = Registry()  # keeps all the tags for any given gene
    for row in fp:
        seqid, gene, tag = row.split()
        if tag == '.':
            registry[gene].append("outside")
        else:
            registry[gene].append("inside")
            if tag[0] == '[':
                registry[gene].append("no_syntenic_model")
                if tag.startswith("[S]"):
                    registry[gene].append("[S]")
                    gstatus = status.get(gene, None)
                    if gstatus == 'complete':
                        registry[gene].append("complete")
                    elif gstatus == 'pseudogene':
                        registry[gene].append("pseudogene")
                    elif gstatus == 'partial':
                        registry[gene].append("partial")
                    else:
                        registry[gene].append("gmap_fail")
                elif tag.startswith("[NS]"):
                    registry[gene].append("[NS]")
                    if "random" in tag or "Scaffold" in tag:
                        registry[gene].append("random")
                    else:
                        registry[gene].append("real_ns")
                elif tag.startswith("[NF]"):
                    registry[gene].append("[NF]")
            else:
                registry[gene].append("syntenic_model")

    inside = registry.count("inside")
    outside = registry.count("outside")
    syntenic = registry.count("syntenic_model")
    non_syntenic = registry.count("no_syntenic_model")
    s = registry.count("[S]")
    ns = registry.count("[NS]")
    nf = registry.count("[NF]")
    complete = registry.count("complete")
    pseudogene = registry.count("pseudogene")
    partial = registry.count("partial")
    gmap_fail = registry.count("gmap_fail")
    random = registry.count("random")
    real_ns = registry.count("real_ns")

    complete_models = registry.get_tag("complete")
    pseudogenes = registry.get_tag("pseudogene")
    partial_deletions = registry.get_tag("partial")

    m = "{0} inside synteny blocks\n".format(inside)
    m += "{0} outside synteny blocks\n".format(outside)
    m += "{0} has syntenic gene\n".format(syntenic)
    m += "{0} lack syntenic gene\n".format(non_syntenic)
    m += "{0} has sequence match in syntenic location\n".format(s)
    m += "{0} has sequence match in non-syntenic location\n".format(ns)
    m += "{0} has sequence match in un-ordered scaffolds\n".format(random)
    m += "{0} has sequence match in real non-syntenic location\n".format(real_ns)
    m += "{0} has no sequence match\n".format(nf)
    m += "{0} syntenic sequence - complete model\n".format(percentage(complete, s))
    m += "{0} syntenic sequence - partial model\n".format(percentage(partial, s))
    m += "{0} syntenic sequence - pseudogene\n".format(percentage(pseudogene, s))
    m += "{0} syntenic sequence - gmap fail\n".format(percentage(gmap_fail, s))
    print(m, file=sys.stderr)

    aa = ["complete_models", "partial_deletions", "pseudogenes"]
    bb = [complete_models, partial_deletions, pseudogenes]
    for a, b in zip(aa, bb):
        fw = open(a, "w")
        print("\n".join(b), file=fw)
        fw.close()

    extra = opts.extra
    if extra:
        registry.update_from(extra)

    fp.seek(0)
    fw = open("registry", "w")
    for row in fp:
        seqid, gene, tag = row.split()
        ts = registry[gene]
        print("\t".join((seqid, gene, tag, "-".join(ts))), file=fw)
    fw.close()

    logging.debug("Registry written.")
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)
815213
def parse_reports(self):
    

    # Set up vars
    self.picard_pcrmetrics_data = dict()
    self.picard_pcrmetrics_samplestats = dict()

    # Go through logs and find Metrics
    for f in self.find_log_files('picard/pcr_metrics', filehandles=True):
        s_name = None
        for l in f['f']:
            # New log starting
            if 'TargetedPcrMetrics' in l and 'INPUT' in l:
                s_name = 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'])

            if s_name is not None:
                if 'TargetedPcrMetrics' in l and '## METRICS CLASS' in l:
                    keys = f['f'].readline().strip("\n").split("\t")
                    vals = f['f'].readline().strip("\n").split("\t")
                    if len(vals) == len(keys):
                        if s_name in self.picard_pcrmetrics_data:
                            log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
                        self.add_data_source(f, s_name, section='TargetedPcrMetrics')
                        self.picard_pcrmetrics_data[s_name] = dict()
                        for i, k in enumerate(keys):
                            try:
                                # Multiply percentages by 100
                                if k.startswith('PCT_'):
                                    vals[i] = float(vals[i]) * 100.0
                                self.picard_pcrmetrics_data[s_name][k] = float(vals[i])
                            except ValueError:
                                self.picard_pcrmetrics_data[s_name][k] = vals[i]

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

    if len(self.picard_pcrmetrics_data) > 0:

        # Write parsed data to a file
        self.write_data_file(self.picard_pcrmetrics_data, 'multiqc_picard_pcrmetrics')

        # Add to general stats table
        self.general_stats_headers['PCT_AMPLIFIED_BASES'] = {
            'title': '% Amplified Bases',
            'description': 'The fraction of aligned bases that mapped to or near an amplicon.',
            'min': 0,
            'max': 100,
            'suffix': '%',
            'scale': 'BrBG'
        }
        self.general_stats_headers['MEDIAN_TARGET_COVERAGE'] = {
            'title': 'Median Target Coverage',
            'description': 'The median coverage of reads that mapped to target regions of an experiment.',
            'min': 0,
            'suffix': 'X',
            'scale': 'GnBu',
        }

        for s_name in self.picard_pcrmetrics_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.picard_pcrmetrics_data[s_name] )

        # Bar plot of ignored bases
        keys = OrderedDict()
        keys['ON_AMPLICON_BASES'] = {'name': 'On-amplicon bases'}
        keys['NEAR_AMPLICON_BASES'] = {'name': 'Near-amplicon bases'}
        keys['OFF_AMPLICON_BASES'] = {'name': 'Off-amplicon bases', 'color': '#f28f43'}

        # Config for the plot
        pconfig = {
            'id': 'picard_pcr_metrics_bases',
            'title': 'Picard: PCR Amplicon Bases',
            'ylab': '# Bases',
            'cpswitch_counts_label': '# Bases',
            'hide_zero_cats': False
        }

        self.add_section (
            name = 'PCR Amplicon Bases',
            anchor = 'picard-pcrmetrics-bases',
            description = 'Metrics about reads obtained from targeted PCR experiments.',
            helptext = '''
            This plot shows the number of bases aligned on or near to amplified regions of the genome.

            * `ON_AMPLICON_BASES`: The number of `PF_BASES_ALIGNED` that mapped to an amplified region of the genome.
            * `NEAR_AMPLICON_BASES`: The number of `PF_BASES_ALIGNED` that mapped to within a fixed interval of an amplified region, but not on a baited region.
            * `OFF_AMPLICON_BASES`: The number of `PF_BASES_ALIGNED` that mapped neither on or near an amplicon.

            For more information see the [Picard documentation](https://broadinstitute.github.io/picard/picard-metric-definitions.html#TargetedPcrMetrics).''',
            plot = bargraph.plot(self.picard_pcrmetrics_data, keys, pconfig)
        )


    # Return the number of detected samples to the parent module
    return len(self.picard_pcrmetrics_data)
349172
def wr_tsv(fout_tsv, tsv_data, **kws):
    """"""
    items_str = "items" if "items" not in kws else kws["items"]
    if tsv_data:
        ifstrm = sys.stdout if fout_tsv is None else open(fout_tsv, 'w')
        num_items = prt_tsv(ifstrm, tsv_data, **kws)
        if fout_tsv is not None:
            sys.stdout.write("  {N:>5} {ITEMS} WROTE: {FOUT}\n".format(
                N=num_items, ITEMS=items_str, FOUT=fout_tsv))
            ifstrm.close()
    else:
        sys.stdout.write("      0 {ITEMS}. NOT WRITING {FOUT}\n".format(
            ITEMS=items_str, FOUT=fout_tsv))
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)
607031
def read_tsv(self):
        
        logging.info('Parsing MOB-recon outputs')
        for sample in self.metadata:
            if os.path.isfile(sample[self.analysistype].contig_report):
                # Read in the .tsv file with pandas. Skip the comment lines
                df = pd.read_csv(sample[self.analysistype].contig_report, delimiter='\t')
                for header in df:
                    # Remove any unwanted whitespace
                    clean_header = header.lstrip().rstrip()
                    # primary_key is the primary key, and value is the value of the cell for that key + header combo
                    for primary_key, value in df[header].items():
                        # Update the dictionary with the new data
                        try:
                            sample[self.analysistype].report_dict[primary_key].update({clean_header: value})
                        # Create the nested dictionary if it hasn't been created yet
                        except KeyError:
                            sample[self.analysistype].report_dict[primary_key] = dict()
                            sample[self.analysistype].report_dict[primary_key].update({clean_header: value})
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
354373
def parse_tstv_by_qual(self):
        

        self.vcftools_tstv_by_qual = dict()
        for f in self.find_log_files('vcftools/tstv_by_qual', filehandles=True):
            d = {}
            for line in f['f'].readlines()[1:]: # don't add the header line (first row)
                key = float(line.split()[0]) # taking the first column (QUAL_THRESHOLD) as key
                val = float(line.split()[6]) # taking Ts/Tv_GT_QUAL_THRESHOLD as value
                if (val == float('inf')) or (val == float('-inf')):
                    val = float('nan')
                d[key] = val
            self.vcftools_tstv_by_qual[f['s_name']] = d

        # Filter out ignored sample names
        self.vcftools_tstv_by_qual = self.ignore_samples(self.vcftools_tstv_by_qual)

        if len(self.vcftools_tstv_by_qual) == 0:
            return 0

        pconfig = {
            'id': 'vcftools_tstv_by_qual',
            'title': 'VCFTools: TsTv by Qual',
            'ylab': 'TsTv Ratio',
            'xlab': 'SNP Quality Threshold',
            'xmin': 0,
            'ymin': 0,
            'smooth_points': 400, # this limits huge filesizes and prevents browser crashing
            'smooth_points_sumcounts': False
        }

        helptext = '''
        `Transition` is a purine-to-purine or pyrimidine-to-pyrimidine point mutations.
        `Transversion` is a purine-to-pyrimidine or pyrimidine-to-purine point mutation.
        `Quality` here is the Phred-scaled quality score as given in the QUAL column of VCF.
        Note: only bi-allelic SNPs are used (multi-allelic sites and INDELs are skipped.)
        Refer to Vcftools's manual (https://vcftools.github.io/man_latest.html) on `--TsTv-by-qual`
        '''

        self.add_section(
            name = 'TsTv by Qual',
            anchor = 'vcftools-tstv-by-qual',
            description = "Plot of `TSTV-BY-QUAL` - the transition to transversion ratio as a function of SNP quality from the output of vcftools TsTv-by-qual.",
            helptext = helptext,
            plot = linegraph.plot(self.vcftools_tstv_by_qual,pconfig)
        )

        return len(self.vcftools_tstv_by_qual)