Find RSeQC bam_stat reports and parse their data
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)
238309
def _get_evidence_bam(work_dir, data):
    
    evidence_bam = glob.glob(os.path.join(work_dir, "results", "evidence",
                                            "evidence_*.%s*.bam" % (dd.get_sample_name(data))))
    if evidence_bam:
        return evidence_bam[0]
869399
def parse_reports(self):
    

    # Get data
    data = _parse_reports_by_type(self)

    if data:
        #  Filter to strip out ignored sample names (REQUIRED)
        data = self.ignore_samples(data)

        # Populate the general stats table
        _add_data_to_general_stats(self, data)

        # Add any found data to the report
        _add_section_to_report(self, data)

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

    self.picard_ValidateSamFile_data = data  # Seems like the right thing to do
    return len(data)
1178128
def make_checksum_validation_script(stats_list):
    """"""
    if not os.path.exists('./hash_check'):
        os.mkdir('./hash_check')

    with open('./hash_check/curl.sh', 'w') as curl_f, open(
        './hash_check/md5.txt', 'w'
    ) as md5_f, open('./hash_check/sha1.txt', 'w') as sha1_f:

        curl_f.write('#!/usr/bin/env bash\n\n')

        for stats_dict in stats_list:
            for sysmeta_xml in stats_dict['largest_sysmeta_xml']:
                print(sysmeta_xml)
                sysmeta_pyxb = d1_common.types.dataoneTypes_v1_2.CreateFromDocument(
                    sysmeta_xml
                )

                pid = sysmeta_pyxb.identifier.value().encode('utf-8')
                file_name = re.sub('\W+', '_', pid)
                size = sysmeta_pyxb.size
                base_url = stats_dict['gmn_dict']['base_url']

                if size > 100 * 1024 * 1024:
                    logging.info('Ignored large object. size={} pid={}')

                curl_f.write('# {} {}\n'.format(size, pid))
                curl_f.write(
                    'curl -o obj/{} {}/v1/object/{}\n'.format(
                        file_name, base_url, d1_common.url.encodePathElement(pid)
                    )
                )

                if sysmeta_pyxb.checksum.algorithm == 'MD5':
                    md5_f.write(
                        '{} obj/{}\n'.format(sysmeta_pyxb.checksum.value(), file_name)
                    )
                else:
                    sha1_f.write(
                        '{} obj/{}\n'.format(sysmeta_pyxb.checksum.value(), file_name)
                    )

    with open('./hash_check/check.sh', 'w') as f:
        f.write('#!/usr/bin/env bash\n\n')
        f.write('mkdir -p obj\n')
        f.write('./curl.sh\n')
        f.write('sha1sum -c sha1.txt\n')
        f.write('md5sum -c md5.txt\n')
1043567
def run(bam_file, data, out_dir):
    
    # Qualimap results should be saved to a directory named after sample.
    # MultiQC (for parsing additional data) picks the sample name after the dir as follows:
    #   <sample name>/raw_data_qualimapReport/insert_size_histogram.txt
    results_dir = os.path.join(out_dir, dd.get_sample_name(data))
    resources = config_utils.get_resources("qualimap", data["config"])
    options = " ".join(resources.get("options", ""))
    results_file = os.path.join(results_dir, "genome_results.txt")
    report_file = os.path.join(results_dir, "qualimapReport.html")
    utils.safe_makedir(results_dir)
    pdf_file = "qualimapReport.pdf"
    if not utils.file_exists(results_file) and not utils.file_exists(os.path.join(results_dir, pdf_file)):
        if "qualimap_full" in tz.get_in(("config", "algorithm", "tools_on"), data, []):
            logger.info("Full qualimap analysis for %s may be slow." % bam_file)
            ds_bam = bam_file
        else:
            ds_bam = bam.downsample(bam_file, data, 1e7, work_dir=out_dir)
            bam_file = ds_bam if ds_bam else bam_file
        if options.find("PDF") > -1:
            options = "%s -outfile %s" % (options, pdf_file)
        num_cores = data["config"]["algorithm"].get("num_cores", 1)
        qualimap = config_utils.get_program("qualimap", data["config"])
        max_mem = config_utils.adjust_memory(resources.get("memory", "1G"),
                                             num_cores)

        with file_transaction(data, results_dir) as tx_results_dir:
            utils.safe_makedir(tx_results_dir)

            export = "%s%s export JAVA_OPTS='-Xms32m -Xmx%s -Djava.io.tmpdir=%s' && " % (
                utils.java_freetype_fix(), utils.local_path_export(), max_mem, tx_results_dir)
            cmd = ("unset DISPLAY && {export} {qualimap} bamqc -bam {bam_file} -outdir {tx_results_dir} "
                   "--skip-duplicated --skip-dup-mode 0 "
                   "-nt {num_cores} {options}")
            species = None
            if (tz.get_in(("genome_resources", "aliases", "human"), data, "")
                  or dd.get_genome_build(data).startswith(("hg", "GRCh"))):
                species = "HUMAN"
            elif dd.get_genome_build(data).startswith(("mm", "GRCm")):
                species = "MOUSE"
            if species in ["HUMAN", "MOUSE"]:
                cmd += " -gd {species}"
            regions = (dd.get_coverage(data) if dd.get_coverage(data) not in [None, False, "None"]
                       else dd.get_variant_regions_merged(data))
            if regions:
                regions = bedutils.merge_overlaps(bedutils.clean_file(regions, data), data)
                bed6_regions = _bed_to_bed6(regions, out_dir)
                cmd += " -gff {bed6_regions}"
            bcbio_env = utils.get_bcbio_env()
            do.run(cmd.format(**locals()), "Qualimap: %s" % dd.get_sample_name(data), env=bcbio_env)
            tx_results_file = os.path.join(tx_results_dir, "genome_results.txt")
            cmd = "sed -i 's/bam file = .*/bam file = %s.bam/' %s" % (dd.get_sample_name(data), tx_results_file)
            do.run(cmd, "Fix Name Qualimap for {}".format(dd.get_sample_name(data)))
    # Qualimap output folder (results_dir) needs to be named after the sample (see comments above). However, in order
    # to keep its name after upload, we need to put  the base QC file (results_file) into the root directory (out_dir):
    base_results_file = os.path.join(out_dir, os.path.basename(results_file))
    shutil.copyfile(results_file, base_results_file)
    return {"base": base_results_file,
            "secondary": _find_qualimap_secondary_files(results_dir, base_results_file)}
277739
def get_stats(self, stat_name):
        

        return [self.get_stat(r, stat_name) for r in self.statistics.keys()]
701388
def parse_reports(self):
    

    # Set up vars
    self.picard_insertSize_data = dict()
    self.picard_insertSize_histogram = dict()
    self.picard_insertSize_samplestats = dict()

    # Go through logs and find Metrics
    for f in self.find_log_files('picard/insertsize', filehandles=True):
        s_name = None
        in_hist = False
        for l in f['f']:

            # Catch the histogram values
            if s_name is not None and in_hist is True:
                try:
                    sections = l.split("\t")
                    ins = int(sections[0])
                    tot_count = sum( [int(x) for x in sections[1:]] )
                    self.picard_insertSize_histogram[s_name][ins] = tot_count
                    self.picard_insertSize_samplestats[s_name]['total_count'] += tot_count
                except ValueError:
                    # Reset in case we have more in this log file
                    s_name = None
                    in_hist = False

            # New log starting
            if 'InsertSizeMetrics' 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 'InsertSizeMetrics' in l and '## METRICS CLASS' in l:
                    if s_name in self.picard_insertSize_data:
                        log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
                    self.add_data_source(f, s_name, section='InsertSizeMetrics')
                    keys = f['f'].readline().strip("\n").split("\t")
                    vals = f['f'].readline().strip("\n").split("\t")
                    self.picard_insertSize_samplestats[s_name] = {'total_count': 0, 'meansum':0, 'total_pairs':0 }
                    orientation_idx = keys.index('PAIR_ORIENTATION')
                    while len(vals) == len(keys):
                        pair_orientation = vals[orientation_idx]
                        rowkey = '{}_{}'.format(s_name, pair_orientation)
                        self.picard_insertSize_data[rowkey] = OrderedDict()
                        self.picard_insertSize_data[rowkey]['SAMPLE_NAME'] = s_name
                        for i, k in enumerate(keys):
                            try:
                                self.picard_insertSize_data[rowkey][k] = float(vals[i])
                            except ValueError:
                                try:
                                    self.picard_insertSize_data[rowkey][k] = float(vals[i].replace(',','.'))
                                    log.debug("Switching commas for points in '{}': {} - {}".format(f['fn'], vals[i], vals[i].replace(',','.')))
                                except ValueError:
                                    self.picard_insertSize_data[rowkey][k] = vals[i]
                            except IndexError:
                                pass # missing data
                        # Add to mean sums
                        rp = self.picard_insertSize_data[rowkey]['READ_PAIRS']
                        mis = self.picard_insertSize_data[rowkey]['MEAN_INSERT_SIZE']
                        self.picard_insertSize_samplestats[s_name]['meansum'] += (rp * mis)
                        self.picard_insertSize_samplestats[s_name]['total_pairs'] += rp

                        vals = f['f'].readline().strip("\n").split("\t")

                    # Skip lines on to histogram
                    l = f['f'].readline().strip("\n")
                    l = f['f'].readline().strip("\n")

                    self.picard_insertSize_histogram[s_name] = OrderedDict()
                    in_hist = True

        for key in list(self.picard_insertSize_data.keys()):
            if len(self.picard_insertSize_data[key]) == 0:
                self.picard_insertSize_data.pop(key, None)
        for s_name in list(self.picard_insertSize_histogram.keys()):
            if len(self.picard_insertSize_histogram[s_name]) == 0:
                self.picard_insertSize_histogram.pop(s_name, None)
                log.debug("Ignoring '{}' histogram as no data parsed".format(s_name))

    # Calculate summed mean values for all read orientations
    for s_name, v in self.picard_insertSize_samplestats.items():
        self.picard_insertSize_samplestats[s_name]['summed_mean'] = v['meansum'] / v['total_pairs']

    # Calculate summed median values for all read orientations
    for s_name in self.picard_insertSize_histogram:
        j = 0
        for idx, c in self.picard_insertSize_histogram[s_name].items():
            j += c
            if j > (self.picard_insertSize_samplestats[s_name]['total_count'] / 2):
                self.picard_insertSize_samplestats[s_name]['summed_median'] = idx
                break


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

    if len(self.picard_insertSize_data) > 0:

        # Write parsed data to a file
        self.write_data_file(self.picard_insertSize_data, 'multiqc_picard_insertSize')

        # Do we have median insert sizes?
        missing_medians = False
        for v in self.picard_insertSize_samplestats.values():
            if 'summed_median' not in v:
                missing_medians = True

        # Add to general stats table
        self.general_stats_headers['summed_median'] = {
            'title': 'Insert Size',
            'description': 'Median Insert Size, all read orientations (bp)',
            'min': 0,
            'suffix': ' bp',
            'format': '{:,.0f}',
            'scale': 'GnBu',
        }
        self.general_stats_headers['summed_mean'] = {
            'title': 'Mean Insert Size',
            'description': 'Mean Insert Size, all read orientations (bp)',
            'min': 0,
            'suffix': ' bp',
            'format': '{:,.0f}',
            'scale': 'GnBu',
            'hidden': False if missing_medians else True
        }
        for s_name in self.picard_insertSize_samplestats:
            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_insertSize_samplestats[s_name] )

        # Section with histogram plot
        if len(self.picard_insertSize_histogram) > 0:
            # Make a normalised percentage version of the data
            data_percent = {}
            for s_name, data in self.picard_insertSize_histogram.items():
                data_percent[s_name] = OrderedDict()
                total = float( sum( data.values() ) )
                for k, v in data.items():
                    data_percent[s_name][k] = (v/total)*100

            # Allow customisation of how smooth the the plot is
            try:
                insertsize_smooth_points = int(config.picard_config['insertsize_smooth_points'])
                log.debug("Custom Picard insert size smoothing: {}".format(insertsize_smooth_points))
            except (AttributeError, KeyError, ValueError):
                insertsize_smooth_points = 500

            # Plot the data and add section
            pconfig = {
                'smooth_points': insertsize_smooth_points,
                'smooth_points_sumcounts': [True, False],
                'id': 'picard_insert_size',
                'title': 'Picard: Insert Size',
                'ylab': 'Count',
                'xlab': 'Insert Size (bp)',
                'xDecimals': False,
                'tt_label': '<b>{point.x} bp</b>: {point.y:.0f}',
                'ymin': 0,
                'data_labels': [
                    {'name': 'Counts', 'ylab': 'Coverage'},
                    {'name': 'Percentages', 'ylab': 'Percentage of Counts'}
                ]
            }
            self.add_section (
                name = 'Insert Size',
                anchor = 'picard-insertsize',
                description = 'Plot shows the number of reads at a given insert size. Reads with different orientations are summed.',
                plot = linegraph.plot([self.picard_insertSize_histogram, data_percent], pconfig)
            )


    # Return the number of detected samples to the parent module
    return len(self.picard_insertSize_data)
411973
def parse_reports(self,
    log_key='picard/markdups',
    section_name='Mark Duplicates',
    section_anchor='picard-markduplicates',
    plot_title='Picard: Deduplication Stats',
    plot_id='picard_deduplication',
    data_filename='multiqc_picard_dups'):
    

    # Set up vars
    self.picard_dupMetrics_data = dict()

    # Go through logs and find Metrics
    for f in self.find_log_files(log_key, filehandles=True):
        s_name = f['s_name']
        for l in f['f']:
            # New log starting
            if 'markduplicates' in l.lower() and 'input' in l.lower():
                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 'UNPAIRED_READ_DUPLICATES' in l:
                    if s_name in self.picard_dupMetrics_data:
                        log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
                    self.add_data_source(f, s_name, section='DuplicationMetrics')
                    self.picard_dupMetrics_data[s_name] = dict()
                    keys = l.rstrip("\n").split("\t")
                    vals = f['f'].readline().rstrip("\n").split("\t")
                    for i, k in enumerate(keys):
                        try:
                            self.picard_dupMetrics_data[s_name][k] = float(vals[i])
                        except ValueError:
                            self.picard_dupMetrics_data[s_name][k] = vals[i]
                    # Check that this sample had some reads
                    if self.picard_dupMetrics_data[s_name].get('READ_PAIRS_EXAMINED', 0) == 0 and \
                       self.picard_dupMetrics_data[s_name].get('UNPAIRED_READS_EXAMINED', 0) == 0:
                        self.picard_dupMetrics_data.pop(s_name, None)
                        log.warn("Skipping MarkDuplicates sample '{}' as log contained no reads".format(s_name))
                    s_name = None

        for s_name in list(self.picard_dupMetrics_data.keys()):
            if len(self.picard_dupMetrics_data[s_name]) == 0:
                self.picard_dupMetrics_data.pop(s_name, None)
                log.debug("Removing {} as no data parsed".format(s_name))


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

    if len(self.picard_dupMetrics_data) > 0:

        # Write parsed data to a file
        self.write_data_file(self.picard_dupMetrics_data, data_filename)

        # Add to general stats table
        self.general_stats_headers['PERCENT_DUPLICATION'] = {
            'title': '% Dups',
            'description': '{} - Percent Duplication'.format(section_name),
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'OrRd',
            'modify': lambda x: self.multiply_hundred(x)
        }
        for s_name in self.picard_dupMetrics_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_dupMetrics_data[s_name] )


        # Make the bar plot and add to the MarkDuplicates section
        # NOTE: I had a hard time getting these numbers to add up as expected.
        # If you think I've done something wrong, let me know! Please add an
        # issue here: https://github.com/ewels/MultiQC/issues
        for sn in self.picard_dupMetrics_data.keys():
            self.picard_dupMetrics_data[sn]['UNPAIRED_READ_UNIQUE'] = self.picard_dupMetrics_data[sn]['UNPAIRED_READS_EXAMINED'] - self.picard_dupMetrics_data[sn]['UNPAIRED_READ_DUPLICATES']
            self.picard_dupMetrics_data[sn]['READ_PAIR_NOT_OPTICAL_DUPLICATES'] = self.picard_dupMetrics_data[sn]['READ_PAIR_DUPLICATES'] - self.picard_dupMetrics_data[sn]['READ_PAIR_OPTICAL_DUPLICATES']
            self.picard_dupMetrics_data[sn]['READ_PAIR_UNIQUE'] = self.picard_dupMetrics_data[sn]['READ_PAIRS_EXAMINED'] - self.picard_dupMetrics_data[sn]['READ_PAIR_DUPLICATES']

        keys = OrderedDict()
        keys_r = ['READ_PAIR_UNIQUE', 'UNPAIRED_READ_UNIQUE', 'READ_PAIR_NOT_OPTICAL_DUPLICATES',
                'READ_PAIR_OPTICAL_DUPLICATES', 'UNPAIRED_READ_DUPLICATES', 'UNMAPPED_READS']
        for k in keys_r:
            keys[k] = {'name': k.replace('_',' ').title()}

        # Config for the plot
        pconfig = {
            'id': plot_id,
            'title': plot_title,
            'ylab': '# Reads',
            'cpswitch_counts_label': 'Number of Reads',
            'cpswitch_c_active': False
        }

        self.add_section (
            name = section_name,
            anchor = section_anchor,
            plot = bargraph.plot(self.picard_dupMetrics_data, keys, pconfig)
        )

    # Return the number of detected samples to the parent module
    return len(self.picard_dupMetrics_data)
686150
def check_bam(bam, o):
    
    try:
        p = sp.Popen(['samtools', 'view', bam], stdout=sp.PIPE)
        # Count paired alignments
        paired = 0
        read_lengths = defaultdict(int)
        while o > 0:  # Count down number of lines
            line = p.stdout.readline().decode().split("\t")
            flag = int(line[1])
            read_lengths[len(line[9])] += 1
            if 1 & flag:  # check decimal flag contains 1 (paired)
                paired += 1
            o -= 1
        p.kill()
    except OSError:
        reason = "Note (samtools not in path): For NGS inputs, " \
                 "pep needs samtools to auto-populate " \
                 "'read_length' and 'read_type' attributes; " \
                 "these attributes were not populated."
        raise OSError(reason)

    _LOGGER.debug("Read lengths: {}".format(read_lengths))
    _LOGGER.debug("paired: {}".format(paired))
    return read_lengths, paired
991145
def bamtag(sam):
    ''' 
    '''
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    if is_python3():
        queryalignment = next(track)
    else:
        queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
1024059
def picard_idxstats(picard, align_bam):
    
    opts = [("INPUT", align_bam)]
    stdout = picard.run("BamIndexStats", opts, get_stdout=True)
    out = []
    AlignInfo = collections.namedtuple("AlignInfo", ["contig", "length", "aligned", "unaligned"])
    for line in stdout.split("\n"):
        if line:
            parts = line.split()
            if len(parts) == 2:
                _, unaligned = parts
                out.append(AlignInfo("nocontig", 0, 0, int(unaligned)))
            elif len(parts) == 7:
                contig, _, length, _, aligned, _, unaligned = parts
                out.append(AlignInfo(contig, int(length), int(aligned), int(unaligned)))
            else:
                raise ValueError("Unexpected output from BamIndexStats: %s" % line)
    return out
861267
def _run_qc_tools(bam_file, data):
    
    from bcbio.qc import (atropos, contamination, coverage, damage, fastqc, kraken,
                          qsignature, qualimap, samtools, picard, srna, umi, variant,
                          viral, preseq, chipseq)
    tools = {"fastqc": fastqc.run,
             "atropos": atropos.run,
             "small-rna": srna.run,
             "samtools": samtools.run,
             "qualimap": qualimap.run,
             "qualimap_rnaseq": qualimap.run_rnaseq,
             "qsignature": qsignature.run,
             "contamination": contamination.run,
             "coverage": coverage.run,
             "damage": damage.run,
             "variants": variant.run,
             "peddy": peddy.run_qc,
             "kraken": kraken.run,
             "picard": picard.run,
             "umi": umi.run,
             "viral": viral.run,
             "preseq": preseq.run,
             "chipqc": chipseq.run
             }
    qc_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "qc", data["description"]))
    metrics = {}
    qc_out = utils.deepish_copy(dd.get_summary_qc(data))
    for program_name in dd.get_algorithm_qc(data):
        if not bam_file and program_name != "kraken":  # kraken doesn't need bam
            continue
        if dd.get_phenotype(data) == "germline" and program_name != "variants":
            continue
        qc_fn = tools[program_name]
        cur_qc_dir = os.path.join(qc_dir, program_name)
        out = qc_fn(bam_file, data, cur_qc_dir)
        qc_files = None
        if out and isinstance(out, dict):
            # Check for metrics output, two cases:
            # 1. output with {"metrics"} and files ("base")
            if "metrics" in out:
                metrics.update(out.pop("metrics"))
            # 2. a dictionary of metrics
            elif "base" not in out:
                metrics.update(out)
            # Check for files only output
            if "base" in out:
                qc_files = out
        elif out and isinstance(out, six.string_types) and os.path.exists(out):
            qc_files = {"base": out, "secondary": []}
        if not qc_files:
            qc_files = _organize_qc_files(program_name, cur_qc_dir)
        if qc_files:
            qc_out[program_name] = qc_files

    metrics["Name"] = dd.get_sample_name(data)
    metrics["Quality format"] = dd.get_quality_format(data).lower()
    return {"qc": qc_out, "metrics": metrics}
1213390
def parse_reports(self):
    

    # Set up vars
    self.junction_annotation_data = dict()
    regexes = {
        'total_splicing_events': r"^Total splicing  Events:\s*(\d+)$",
        'known_splicing_events': r"^Known Splicing Events:\s*(\d+)$",
        'partial_novel_splicing_events': r"^Partial Novel Splicing Events:\s*(\d+)$",
        'novel_splicing_events': r"^Novel Splicing Events:\s*(\d+)$",
        'total_splicing_junctions': r"^Total splicing  Junctions:\s*(\d+)$",
        'known_splicing_junctions': r"^Known Splicing Junctions:\s*(\d+)$",
        'partial_novel_splicing_junctions': r"^Partial Novel Splicing Junctions:\s*(\d+)$",
        'novel_splicing_junctions': r"^Novel Splicing Junctions:\s*(\d+)$",
    }

    # Go through files and parse data using regexes
    for f in self.find_log_files('rseqc/junction_annotation'):
        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_splicing_events' in d:
            t = float(d['total_splicing_events'])
            if 'known_splicing_events' in d:
                d['known_splicing_events_pct'] = (float(d['known_splicing_events']) / t)*100.0
            if 'partial_novel_splicing_events' in d:
                d['partial_novel_splicing_events_pct'] = (float(d['partial_novel_splicing_events']) / t)*100.0
            if 'novel_splicing_events' in d:
                d['novel_splicing_events_pct'] = (float(d['novel_splicing_events']) / t)*100.0
        if 'total_splicing_junctions' in d:
            t = float(d['total_splicing_junctions'])
            if 'known_splicing_junctions' in d:
                d['known_splicing_junctions_pct'] = (float(d['known_splicing_junctions']) / t)*100.0
            if 'partial_novel_splicing_junctions' in d:
                d['partial_novel_splicing_junctions_pct'] = (float(d['partial_novel_splicing_junctions']) / t)*100.0
            if 'novel_splicing_junctions' in d:
                d['novel_splicing_junctions_pct'] = (float(d['novel_splicing_junctions']) / t)*100.0

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

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

    if len(self.junction_annotation_data) > 0:

        # Write to file
        self.write_data_file(self.junction_annotation_data, 'multiqc_rseqc_junction_annotation')

        # Plot junction annotations
        keys = [OrderedDict(), OrderedDict()]
        keys[0]['known_splicing_junctions'] = { 'name': 'Known Splicing Junctions' }
        keys[0]['partial_novel_splicing_junctions'] = { 'name': 'Partial Novel Splicing Junctions' }
        keys[0]['novel_splicing_junctions'] = { 'name': 'Novel Splicing Junctions' }
        keys[1]['known_splicing_events'] = { 'name': 'Known Splicing Events' }
        keys[1]['partial_novel_splicing_events'] = { 'name': 'Partial Novel Splicing Events' }
        keys[1]['novel_splicing_events'] = { 'name': 'Novel Splicing Events' }

        pconfig = {
            'id': 'rseqc_junction_annotation_junctions_plot',
            'title': 'RSeQC: Splicing Junctions',
            'ylab': '% Junctions',
            'cpswitch_c_active': False,
            'data_labels': [ 'Junctions', 'Events' ]
        }
        self.add_section (
            name = 'Junction Annotation',
            anchor = 'rseqc_junction_annotation',
            description = '<a href="http://rseqc.sourceforge.net/#junction-annotation-py" target="_blank">Junction annotation</a>' \
                " compares detected splice junctions to" \
                " a reference gene model. An RNA read can be spliced 2" \
                " or more times, each time is called a splicing event.",
            plot = bargraph.plot([self.junction_annotation_data, self.junction_annotation_data], keys, pconfig)
        )

    # Return number of samples found
    return len(self.junction_annotation_data)
299119
def detect_bam_features(bamfile, n_entries=1000):
    '''  '''

    inbam = pysam.Samfile(bamfile)
    inbam = inbam.fetch(until_eof=True)

    tags = ["NH", "X0", "XT"]
    available_tags = {x: 1 for x in tags}

    for n, read in enumerate(inbam):
        if n > n_entries:
            break

        if read.is_unmapped:
            continue

        else:
            for tag in tags:
                if not read.has_tag(tag):
                    available_tags[tag] = 0

    return available_tags
903611
def parse_reports(self):
    

    # Set up vars
    self.picard_GCbias_data = dict()
    self.picard_GCbiasSummary_data = dict()

    # Go through logs and find Metrics
    for f in self.find_log_files('picard/gcbias', filehandles=True):
        s_name = None
        gc_col = None
        cov_col = None
        for l in f['f']:
            # New log starting
            if 'GcBiasMetrics' 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 gc_col is not None and cov_col is not None :
                    try:
                        # Note that GC isn't always the first column.
                        s = l.strip("\n").split("\t")
                        self.picard_GCbias_data[s_name][ int(s[gc_col]) ] = float(s[cov_col])
                    except IndexError:
                        s_name = None
                        gc_col = None
                        cov_col = None

                if 'GcBiasDetailMetrics' in l and '## METRICS CLASS' in l:
                    if s_name in self.picard_GCbias_data:
                        log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
                    self.add_data_source(f, s_name, section='GcBiasDetailMetrics')
                    self.picard_GCbias_data[s_name] = dict()
                    # Get header - find columns with the data we want
                    l = f['f'].readline()
                    s = l.strip("\n").split("\t")
                    gc_col = s.index('GC')
                    cov_col = s.index('NORMALIZED_COVERAGE')

                if 'GcBiasSummaryMetrics' in l and '## METRICS CLASS' in l:
                    if s_name in self.picard_GCbias_data:
                        log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
                    self.add_data_source(f, s_name, section='GcBiasSummaryMetrics')
                    self.picard_GCbiasSummary_data[s_name] = dict()

                    keys = f['f'].readline().rstrip("\n").split("\t")
                    vals = f['f'].readline().rstrip("\n").split("\t")
                    for i, k in enumerate(keys):
                        try:
                            self.picard_GCbiasSummary_data[s_name][k] = float(vals[i])
                        except ValueError:
                            self.picard_GCbiasSummary_data[s_name][k] = vals[i]


        for s_name in list(self.picard_GCbias_data.keys()):
            if len(self.picard_GCbias_data[s_name]) == 0:
                self.picard_GCbias_data.pop(s_name, None)
                log.debug("Removing {} as no data parsed".format(s_name))

        for s_name in list(self.picard_GCbiasSummary_data.keys()):
            if len(self.picard_GCbiasSummary_data[s_name]) == 0:
                self.picard_GCbiasSummary_data.pop(s_name, None)
                log.debug("Removing {} as no data parsed".format(s_name))


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

    if len(self.picard_GCbias_data) > 0:

        # Plot the graph

        pconfig = {
            'id': 'picard_gcbias_plot',
            'title': 'Picard: GC Coverage Bias',
            'ylab': 'Normalized Coverage',
            'xlab': '% GC',
            'xmin': 0,
            'xmax': 100,
            'xDecimals': False,
            'ymin': 0,
            'yCeiling': 10,
            'tt_label': '<b>{point.x} %GC</b>: {point.y:.2f}',
            'yPlotLines': [
                {'value': 1, 'color': '#999999', 'width': 2, 'dashStyle': 'LongDash'},
            ]
        }
        self.add_section (
            name = 'GC Coverage Bias',
            anchor = 'picard-gcbias',
            description = 'This plot shows bias in coverage across regions of the genome with varying GC content.'\
                ' A perfect library would be a flat line at <code>y = 1</code>.',
            plot = linegraph.plot(self.picard_GCbias_data, pconfig)
        )

    if len(self.picard_GCbiasSummary_data) > 0:
        # Write parsed summary data to a file
        self.write_data_file(self.picard_GCbiasSummary_data, 'multiqc_picard_gcbias')


    # Return the number of detected samples to the parent module
    return len(self.picard_GCbias_data)
909802
def run(bam, chrom, pos1, pos2, reffa, chr_reffa, parameters):
    """"""

    # check for chr ref
    is_chr_query = chrom.startswith('chr')
    if is_chr_query and chr_reffa is None:
        chr_reffa = reffa

    # check bam ref type
    bam_header = subprocess.check_output("samtools view -H {}".format(bam), shell=True)
    is_chr_bam = bam_header.find('SN:chr') != -1
    if is_chr_bam:
        reffa = chr_reffa
    if not is_chr_query and is_chr_bam:
        chrom = 'chr' + chrom
    if is_chr_query and not is_chr_bam:
        chrom = re.sub(r'^chr', '', chrom)

    posmin = min(pos1, pos2)
    posmax = max(pos1, pos2)
    cmd = "samtools view -bh {bam} {chrom}:{pos1}-{pos2} " \
        "| samtools mpileup {parameters} -f {reffa} -".format(bam=bam, chrom=chrom,
                                                              pos1=posmin, pos2=posmax,
                                                              reffa=reffa, parameters=parameters)
    if pos1 == pos2:
        cmd += " | awk '$2 == {pos}'".format(pos=pos1)
    else:
        cmd += " | tail -n +2 | awk '$2 >= {posmin} && $2 <= {posmax}'".format(posmin=posmin, posmax=posmax)
    sys.stderr.write("Running:\n{}\n".format(cmd))
    child = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE)
    stdout, stderr = child.communicate()
    if child.returncode != 0:
        if len(stdout) == 0 and stderr is None:
            warnings.warn("Command:\n{cmd}\n did not exit with zero exit code. "
                          "Possibly no coverage for sample.".format(cmd=cmd))
        else:
            raise(Exception("Command:\n{cmd}\n did not exit with zero exit code. "
                            "Check command.".format(cmd=cmd)))
    else:
        return stdout
229054
def parse_hicpro_stats(self, f, rsection):
        
        s_name = self.clean_s_name(os.path.basename(f['root']), os.path.dirname(f['root']))

        if s_name not in self.hicpro_data.keys():
            self.hicpro_data[s_name] = {}

        self.add_data_source(f, s_name, section=rsection)
        for l in f['f'].splitlines():
            if not l.startswith('#'):
                s = l.split("\t")
                if s[0] in self.hicpro_data[s_name]:
                    log.debug("Duplicated keys found! Overwriting: {}".format(s[0]))
                self.hicpro_data[s_name][s[0]] = int(s[1])