Parse all kraken-report data files into sample counts dict
    and store global taxon id -> taxonomy data
498134
def process_samples(kraken_reports_fp, max_rank, min_rank):
    
    taxa = OrderedDict()
    sample_counts = OrderedDict()
    for krep_fp in kraken_reports_fp:
        if not osp.isfile(krep_fp):
            raise RuntimeError("ERROR: File '{}' not found.".format(krep_fp))

        # use the kraken report filename as the sample ID
        sample_id = osp.splitext(osp.split(krep_fp)[1])[0]

        with open(krep_fp, "rt") as kf:
            try:
                kdr = csv.DictReader(kf, fieldnames=field_names, 
                                     delimiter="\t")
                kdata = [entry for entry in kdr][1:]
            except OSError as oe:
                raise RuntimeError("ERROR: {}".format(oe))

        scounts, staxa = parse_kraken_report(kdata, max_rank=max_rank, 
                                             min_rank=min_rank)

        # update master records
        taxa.update(staxa)
        sample_counts[sample_id] = scounts

    return sample_counts, taxa
1134885
def parse_pointfinder(self):
        
        # Create the nested dictionary that stores the necessary values for creating summary reports
        self.populate_summary_dict()
        # Clear out any previous reports
        for organism in self.summary_dict:
            for report in self.summary_dict[organism]:
                try:
                    os.remove(self.summary_dict[organism][report]['summary'])
                except FileNotFoundError:
                    pass
        for sample in self.runmetadata.samples:
            # Find the PointFinder outputs. If the outputs don't exist, create the appropriate entries in the
            # summary dictionary as required
            try:
                self.summary_dict[sample.general.referencegenus]['prediction']['output'] = \
                    glob(os.path.join(sample[self.analysistype].pointfinder_outputs, '{seq}*prediction.txt'
                                      .format(seq=sample.name)))[0]
            except IndexError:
                try:
                    self.summary_dict[sample.general.referencegenus]['prediction']['output'] = str()
                except KeyError:
                    self.populate_summary_dict(genus=sample.general.referencegenus,
                                               key='prediction')
            try:
                self.summary_dict[sample.general.referencegenus]['table']['output'] = \
                    glob(os.path.join(sample[self.analysistype].pointfinder_outputs, '{seq}*table.txt'
                                      .format(seq=sample.name)))[0]
            except IndexError:
                try:
                    self.summary_dict[sample.general.referencegenus]['table']['output'] = str()
                except KeyError:
                    self.populate_summary_dict(genus=sample.general.referencegenus,
                                               key='table')
            try:
                self.summary_dict[sample.general.referencegenus]['results']['output'] = \
                    glob(os.path.join(sample[self.analysistype].pointfinder_outputs, '{seq}*results.tsv'
                                      .format(seq=sample.name)))[0]
            except IndexError:
                try:
                    self.summary_dict[sample.general.referencegenus]['results']['output'] = str()
                except KeyError:
                    self.populate_summary_dict(genus=sample.general.referencegenus,
                                               key='results')
            # Process the predictions
            self.write_report(summary_dict=self.summary_dict,
                              seqid=sample.name,
                              genus=sample.general.referencegenus,
                              key='prediction')
            # Process the results summary
            self.write_report(summary_dict=self.summary_dict,
                              seqid=sample.name,
                              genus=sample.general.referencegenus,
                              key='results')

            # Process the table summary
            self.write_table_report(summary_dict=self.summary_dict,
                                    seqid=sample.name,
                                    genus=sample.general.referencegenus)
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)
845600
def collect(self, file=sys.stderr):
        
        cur = gcstats()
        Ncur = len(cur)
        if self.stats is not None and file is not None:
            prev = self.stats
            Nprev = self.ntypes  # may be less than len(prev)

            if Ncur != Nprev:
                print("# Types %d -> %d" % (Nprev, Ncur), file=file)

            Scur, Sprev, first = set(cur), set(prev), True
            for T in Scur - Sprev:  # new types
                if first:
                    print('New Types', file=file)
                    first = False
                print(' ', T, cur[T], file=file)

            first = True
            for T in Sprev - Scur:  # collected types
                if first:
                    print('Cleaned Types', file=file)
                    first = False
                print(' ', T, -prev[T], file=file)

            first = True
            for T in Scur & Sprev:
                if cur[T] == prev[T]:
                    continue
                if first:
                    print('Known Types', file=file)
                    first = False
                print(' ', T, cur[T], 'delta', cur[T] - prev[T], file=file)

        else:  # first call
            print("All Types", file=file)
            for T, C in cur.items():
                print(' ', T, C, file=file)

        self.stats, self.ntypes = cur, len(cur)
458334
def reporter(self, analysistype='genesippr'):
        
        logging.info('Creating {} report'.format(analysistype))
        # Create a dictionary to link all the genera with their genes
        genusgenes = dict()
        # The organism-specific targets are in .tfa files in the target path
        targetpath = str()
        for sample in self.runmetadata.samples:
            if sample.general.bestassemblyfile != 'NA':
                targetpath = sample[analysistype].targetpath
        for organismfile in glob(os.path.join(targetpath, '*.tfa')):
            organism = os.path.splitext(os.path.basename(organismfile))[0]
            # Use BioPython to extract all the gene names from the file
            for record in SeqIO.parse(open(organismfile), 'fasta'):
                # Append the gene names to the genus-specific list
                try:
                    genusgenes[organism].add(record.id.split('_')[0])
                except (KeyError, IndexError):
                    genusgenes[organism] = set()
                    genusgenes[organism].add(record.id.split('_')[0])
        # Determine from which genera the gene hits were sourced
        for sample in self.runmetadata.samples:
            # Initialise the list to store the genera
            sample[analysistype].targetgenera = list()
            if sample.general.bestassemblyfile != 'NA':
                for organism in genusgenes:
                    # Iterate through all the genesippr hits and attribute each gene to the appropriate genus
                    for gene in sample[analysistype].results:
                        # If the gene name is in the genes from that organism, add the genus name to the list of
                        # genera found in the sample
                        if gene.split('_')[0] in genusgenes[organism]:
                            if organism not in sample[analysistype].targetgenera:
                                sample[analysistype].targetgenera.append(organism)
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        # The report will have every gene for all genera in the header
        header = 'Strain,Genus,{},\n'.format(','.join(self.genelist))
        data = str()
        with open(os.path.join(self.reportpath, analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                sample[analysistype].report_output = list()
                if sample.general.bestassemblyfile != 'NA':
                    # Add the genus/genera found in the sample
                    data += '{},{},'.format(sample.name, ';'.join(sample[analysistype].targetgenera))
                    best_dict = dict()
                    if sample[analysistype].results:
                        gene_check = list()
                        # Find the best match for all the hits
                        for target, pid in sample[analysistype].results.items():
                            gene_name = target.split('_')[0]
                            for gene in self.genelist:
                                # If the key matches a gene in the list of genes
                                if gene == gene_name:
                                    # If the percent identity is better, update the dictionary
                                    try:
                                        if float(pid) > best_dict[gene]:
                                            best_dict[gene] = float(pid)
                                    except KeyError:
                                        best_dict[gene] = float(pid)
                        for gene in self.genelist:
                            # If the gene was not found in the sample, print an empty cell in the report
                            try:
                                best_dict[gene]
                            except KeyError:
                                data += ','
                            # Print the required information for the gene
                            for name, identity in sample[analysistype].results.items():
                                if name.split('_')[0] == gene and gene not in gene_check:
                                    data += '{pid}%'.format(pid=best_dict[gene])
                                    try:
                                        if not sample.general.trimmedcorrectedfastqfiles[0].endswith('.fasta'):
                                            data += ' ({avgd} +/- {std}),'\
                                                .format(avgd=sample[analysistype].avgdepth[name],
                                                        std=sample[analysistype].standarddev[name])
                                        else:
                                            data += ','
                                    except IndexError:
                                        data += ','
                                    gene_check.append(gene)
                                    # Add the simplified results to the object - used in the assembly pipeline report
                                    sample[analysistype].report_output.append(gene)
                        # Add a newline after each sample
                        data += '\n'
                    # Add a newline if the sample did not have any gene hits
                    else:
                        data += '\n'
            # Write the header and data to file
            report.write(header)
            report.write(data)
1051162
def _extract_number_of_taxa(self):
        
        n_taxa = dict()
        for i in self.seq_records:
            if i.gene_code not in n_taxa:
                n_taxa[i.gene_code] = 0
            n_taxa[i.gene_code] += 1
        number_taxa = sorted([i for i in n_taxa.values()], reverse=True)[0]
        self.number_taxa = str(number_taxa)
93135
def _set_global_vars(metadata):
    
    fnames = collections.defaultdict(list)
    for sample in metadata.keys():
        for k, v in metadata[sample].items():
            if isinstance(v, six.string_types) and os.path.isfile(v):
                v = _expand_file(v)
                metadata[sample][k] = v
                fnames[v].append(k)
    global_vars = {}
    # Skip global vars -- more confusing than useful
    # loc_counts = collections.defaultdict(int)
    # global_var_sub = {}
    # for fname, locs in fnames.items():
    #     if len(locs) > 1:
    #         loc_counts[locs[0]] += 1
    #         name = "%s%s" % (locs[0], loc_counts[locs[0]])
    #         global_var_sub[fname] = name
    #         global_vars[name] = fname
    # for sample in metadata.keys():
    #     for k, v in metadata[sample].items():
    #         if isinstance(v, six.string_types) and v in global_var_sub:
    #             metadata[sample][k] = global_var_sub[v]
    return metadata, global_vars
406143
def _get_chrbands(self, limit, taxon):
        
        model = Model(self.graph)
        line_counter = 0
        myfile = '/'.join((self.rawdir, self.files[taxon]['file']))
        LOG.info("Processing Chr bands from FILE: %s", myfile)
        geno = Genotype(self.graph)

        # build the organism's genome from the taxon
        genome_label = self.files[taxon]['genome_label']
        taxon_id = 'NCBITaxon:' + taxon

        # add the taxon as a class.  adding the class label elsewhere
        model.addClassToGraph(taxon_id, None)
        model.addSynonym(taxon_id, genome_label)

        genome_id = geno.makeGenomeID(taxon_id)
        geno.addGenome(taxon_id, genome_label)
        model.addOWLPropertyClassRestriction(
            genome_id, self.globaltt['in taxon'], taxon_id)

        placed_scaffold_pattern = r'chr(\d+|X|Y|Z|W|MT|M)'
        # currently unused patterns
        # unlocalized_scaffold_pattern = placed_scaffold_pattern + r'_(\w+)_random'
        # unplaced_scaffold_pattern = r'chrUn_(\w+)'

        col = ['chrom', 'start', 'stop', 'band', 'rtype']
        with gzip.open(myfile, 'rb') as reader:
            for line in reader:
                line_counter += 1
                # skip comments
                line = line.decode().strip()
                if line[0] == '#':
                    continue
                # chr13	4500000	10000000	p12	stalk
                row = line.split('\t')
                chrom = row[col.index('chrom')]
                band = row[col.index('band')]
                rtype = row[col.index('rtype')]
                # NOTE
                # some less-finished genomes have placed and unplaced scaffolds
                # * Placed scaffolds:
                #    Scaffold has an oriented location within a chromosome.
                # * Unlocalized scaffolds:
                #     scaffold 's chromosome  is known,
                #     scaffold's position, orientation or both is not known.
                # *Unplaced scaffolds:
                #   it is not known which chromosome the scaffold belongs to.
                # find out if the thing is a full on chromosome, or a scaffold:
                # ex: unlocalized scaffold: chr10_KL568008v1_random
                # ex: unplaced scaffold: chrUn_AABR07022428v1

                mch = re.match(placed_scaffold_pattern+r'$', chrom)
                if mch is not None and len(mch.groups()) == 1:
                    # the chromosome is the first match of the pattern
                    # chrom = m.group(1)  # TODO unused
                    pass
                else:
                    # let's skip over anything that isn't a placed_scaffold
                    LOG.info("Skipping non-placed chromosome %s", chrom)
                    continue
                # the chrom class, taxon as the reference
                cclassid = makeChromID(chrom, taxon, 'CHR')

                # add the chromosome as a class
                geno.addChromosomeClass(chrom, taxon_id, genome_label)
                model.addOWLPropertyClassRestriction(
                    cclassid, self.globaltt['member of'], genome_id)

                # add the band(region) as a class
                maplocclass_id = cclassid+band
                maplocclass_label = makeChromLabel(chrom+band, genome_label)
                if band is not None and band.strip() != '':
                    region_type_id = self.map_type_of_region(rtype)
                    model.addClassToGraph(
                        maplocclass_id, maplocclass_label,
                        region_type_id)
                else:
                    region_type_id = self.globaltt['chromosome']
                # add the staining intensity of the band
                if re.match(r'g(neg|pos|var)', rtype):
                    if region_type_id in [
                            self.globaltt['chromosome_band'],
                            self.globaltt['chromosome_subband']]:
                        stain_type = self.resolve(rtype)
                        if stain_type is not None:
                            model.addOWLPropertyClassRestriction(
                                maplocclass_id,
                                self.globaltt['has_sequence_attribute'],
                                self.resolve(rtype))
                    else:
                        # usually happens if it's a chromosome because
                        # they don't actually have banding info
                        LOG.info("feature type %s != chr band", region_type_id)
                else:
                    LOG.warning('staining type not found: %s', rtype)

                # get the parent bands, and make them unique
                parents = list(self.make_parent_bands(band, set()))
                # alphabetical sort will put them in smallest to biggest
                parents.sort(reverse=True)

                # print("PARENTS of", maplocclass_id, "=", parents)
                # add the parents to the graph, in hierarchical order
                # TODO this is somewhat inefficient due to
                # re-adding upper-level nodes when iterating over the file
                for prnt in parents:
                    parent = prnt.strip()
                    if parent is None or parent == "":
                        continue
                    pclassid = cclassid + parent  # class chr parts
                    pclass_label = makeChromLabel(chrom + parent, genome_label)
                    rti = getChrPartTypeByNotation(parent, self.graph)
                    model.addClassToGraph(pclassid, pclass_label, rti)

                    # for canonical chromosomes,
                    # then the subbands are subsequences of the full band
                    # add the subsequence stuff as restrictions

                    if prnt != parents[-1]:
                        grandparent = 1 + parents.index(prnt)
                        pid = cclassid + parents[grandparent]   # the instance
                        model.addOWLPropertyClassRestriction(
                            pclassid, self.globaltt['is subsequence of'], pid)
                        model.addOWLPropertyClassRestriction(
                            pid, self.globaltt['has subsequence'], pclassid)
                    else:
                        # add the last one (p or q usually)
                        # as attached to the chromosome
                        model.addOWLPropertyClassRestriction(
                            pclassid, self.globaltt['is subsequence of'], cclassid)
                        model.addOWLPropertyClassRestriction(
                            cclassid, self.globaltt['has subsequence'], pclassid)

                # connect the band here to the first one in the parent list
                if len(parents) > 0:
                    model.addOWLPropertyClassRestriction(
                        maplocclass_id, self.globaltt['is subsequence of'],
                        cclassid + parents[0])
                    model.addOWLPropertyClassRestriction(
                        cclassid + parents[0], self.globaltt['has subsequence'],
                        maplocclass_id)

                if limit is not None and line_counter > limit:
                    break

        # TODO figure out the staining intensities for the encompassing bands

        return
916854
def _parse_reports_by_type(self):
    

    data = dict()

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

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

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

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

        data[sample] = sample_data

    return data
791048
def parse_reports(self):
    

    # Set up vars
    self.picard_rrbs_metrics = dict()

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

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

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

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

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

    if len(self.picard_rrbs_metrics) > 0:

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

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



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

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

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

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


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

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

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

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

    # Return the number of detected samples to the parent module
    return len(self.picard_rrbs_metrics)
863915
def main():
    
    # Get commandline args		
    args = get_args()
	
    # Instantiate the ete NCBI taxa object
    ncbi = NCBITaxa(dbfile=args.database)
    ## dbfile location
    if args.verbose > 1:
        sys.stderr.write('Taxa database is stored at {}\n'.format(ncbi.dbfile))

    # Update the database if required.
    if args.update is True:
        if args.verbose > 1:
            msg = 'Updating the taxonomy database. This may take several minutes...\n'
            sys.stderr.write(msg)
        ncbi.update_taxonomy_database()
            
    # If names were provided in taxid list, convert to taxids
    args.taxid = args.taxid.replace('"', '').replace("'", '').split(',')
    args.taxid = name2taxid(args.taxid, ncbi)

    # Output
    if args.outfile is None:
        outFH = sys.stdout
    else:
        outFH = open(args.outfile, 'w')
    ## header
    if args.taxon_info:
        outFH.write('\t'.join(['name', 'taxid', 'rank', 'lineage']) + '\n')
    elif not args.just_taxids:
        outFH.write('\t'.join(['parent_taxid',
                               'descendent_taxid',
                               'descendent_name']) + '\n')
    ## body
    for taxid in args.taxid:
        if args.taxon_info:
            taxon_info(taxid, ncbi, outFH)
        else:
            desc_taxa(taxid, ncbi,  outFH, args.just_taxids)
            
    outFH.close()
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)
529817
def parse_single_report(f):
    
    # Fixme: Separate GATKReport parsing and data subsetting. A GATKReport parser now available from the GATK MultiqcModel.

    data = dict()
    in_CompOverlap = False
    in_CountVariants = False
    in_TiTv = False
    for l in f:
        # Detect section headers
        if '#:GATKTable:CompOverlap' in l:
            in_CompOverlap = True
        elif '#:GATKTable:CountVariants' in l:
            in_CountVariants = True
        elif '#:GATKTable:TiTvVariantEvaluator' in l:
            in_TiTv = True
        else:
            # Parse contents using nested loops
            if in_CompOverlap:
                headers = l.split()
                while in_CompOverlap:
                    l = f.readline().strip("\n")
                    d = dict()
                    try:
                        for i, s in enumerate(l.split()):
                            d[headers[i]] = s
                        if d['Novelty'] == 'all':
                            data['reference'] = d['CompRod']
                            data['comp_rate'] = float(d['compRate'])
                            data['concordant_rate'] = float(d['concordantRate'])
                            data['eval_variants'] = int(d['nEvalVariants'])
                            data['novel_sites'] = int(d['novelSites'])
                        elif d['Novelty'] == 'known':
                            data['known_sites'] = int(d['nEvalVariants'])
                    except KeyError:
                        in_CompOverlap = False
            elif in_CountVariants:
                headers = l.split()
                while in_CountVariants:
                    l = f.readline().strip("\n")
                    d = dict()
                    try:
                        for i, s in enumerate(l.split()):
                            d[headers[i]] = s
                        if d['Novelty'] == 'all':
                            data['snps'] = int(d['nSNPs'])
                            data['mnps'] = int(d['nMNPs'])
                            data['insertions'] = int(d['nInsertions'])
                            data['deletions'] = int(d['nDeletions'])
                            data['complex'] = int(d['nComplex'])
                            data['symbolic'] = int(d['nSymbolic'])
                            data['mixed'] = int(d['nMixed'])
                            data['nocalls'] = int(d['nNoCalls'])
                    except KeyError:
                        in_CountVariants = False
            elif in_TiTv:
                headers = l.split()
                data['titv_reference'] = 'unknown'
                while in_TiTv:
                    l = f.readline().strip("\n")
                    d = dict()
                    try:
                        for i, s in enumerate(l.split()):
                            d[headers[i]] = s
                        if d['Novelty'] == 'known':
                            data['titv_reference'] = d['CompRod']
                            data['known_titv'] = float(d['tiTvRatio'])
                        elif d['Novelty'] == 'novel':
                            data['novel_titv'] = float(d['tiTvRatio'])
                    except KeyError:
                        in_TiTv = False

    return data
1203407
def parse_sargasso_logs(self, f):
        
        species_name = list()
        items = list()
        header = list()
        is_first_line = True
        for l in f['f'].splitlines():
            s = l.split(",")
            # Check that this actually is a Sargasso file
            if is_first_line and s[0]!='Sample':
                return None

            if len(s) < 7:
                continue
            if is_first_line:
                #prepare header
                is_first_line = False
                header = s
                for i in header[1:]:
                    #find out what species included
                    sname = i.split('-')[-1]
                    if sname not in species_name:
                        species_name.append(sname)
                    #find out what is being counted
                    kname = ("-".join(i.split('-')[-3:-1]))
                    if kname not in items:
                        items.append(kname)
            else:
                #start sample lines.
                sample_name = s.pop(0)

                chunk_by_species = [s[i:i + len(items)] for i in range(0, len(s), len(items))];
                for idx,v in enumerate(chunk_by_species):
                    #adding species name to the sample name for easy interpretation
                    new_sample_name = '_'.join([sample_name,species_name[idx]])

                    # Clean up sample name
                    new_sample_name = self.clean_s_name(new_sample_name, f['root'])

                    if new_sample_name in self.sargasso_data.keys():
                        log.debug("Duplicate sample name found! Overwriting: {}".format(new_sample_name))

                    try:
                        self.sargasso_data[new_sample_name] = dict(zip(items,map(int, v)))
                    except ValueError:
                        pass



        self.sargasso_keys = items

        for idx, f_name in enumerate(self.sargasso_data.keys()):



            # Reorganised parsed data for this sample
            # Collect total READ count number
            self.sargasso_data[f_name]['Total'] = 0;
            for key, value in list(self.sargasso_data[f_name].items()):   # iter on both keys and values
                if key.endswith("Reads"):
                    self.sargasso_data[f_name]['Total'] += value

            # Calculate the percent aligned if we can
            try:
                self.sargasso_data[f_name]['sargasso_percent_assigned'] = (float(self.sargasso_data[f_name]['Assigned-Reads'])/float(self.sargasso_data[f_name]['Total'])) * 100.0
            except (KeyError, ZeroDivisionError):
                pass
1036647
def blastparse(self):
        
        logging.info('Parsing BLAST results')
        # Load the NCBI 16S reference database as a dictionary
        for sample in self.runmetadata.samples:
            if sample.general.bestassemblyfile != 'NA':
                # Load the NCBI 16S reference database as a dictionary
                dbrecords = SeqIO.to_dict(SeqIO.parse(sample[self.analysistype].baitfile, 'fasta'))
                # Allow for no BLAST results
                if os.path.isfile(sample[self.analysistype].blastreport):
                    # Initialise a dictionary to store the number of times a genus is the best hit
                    sample[self.analysistype].frequency = dict()
                    # Open the sequence profile file as a dictionary
                    blastdict = DictReader(open(sample[self.analysistype].blastreport),
                                           fieldnames=self.fieldnames, dialect='excel-tab')
                    recorddict = dict()
                    for record in blastdict:
                        # Create the subject id. It will look like this: gi|1018196593|ref|NR_136472.1|
                        subject = record['subject_id']
                        # Extract the genus name. Use the subject id as a key in the dictionary of the reference db.
                        # It will return the full record e.g. gi|1018196593|ref|NR_136472.1| Escherichia marmotae
                        # strain HT073016 16S ribosomal RNA, partial sequence
                        # This full description can be manipulated to extract the genus e.g. Escherichia
                        genus = dbrecords[subject].description.split('|')[-1].split()[0]
                        # Increment the number of times this genus was found, or initialise the dictionary with this
                        # genus the first time it is seen
                        try:
                            sample[self.analysistype].frequency[genus] += 1
                        except KeyError:
                            sample[self.analysistype].frequency[genus] = 1
                        try:
                            recorddict[dbrecords[subject].description] += 1
                        except KeyError:
                            recorddict[dbrecords[subject].description] = 1
                    # Sort the dictionary based on the number of times a genus is seen
                    sample[self.analysistype].sortedgenera = sorted(sample[self.analysistype].frequency.items(),
                                                                    key=operator.itemgetter(1), reverse=True)
                    try:
                        # Extract the top result, and set it as the genus of the sample
                        sample[self.analysistype].genus = sample[self.analysistype].sortedgenera[0][0]
                        # Previous code relies on having the closest refseq genus, so set this as above
                        # sample.general.closestrefseqgenus = sample[self.analysistype].genus
                    except IndexError:
                        # Populate attributes with 'NA'
                        sample[self.analysistype].sortedgenera = 'NA'
                        sample[self.analysistype].genus = 'NA'
                        # sample.general.closestrefseqgenus = 'NA'
                else:
                    # Populate attributes with 'NA'
                    sample[self.analysistype].sortedgenera = 'NA'
                    sample[self.analysistype].genus = 'NA'
                    # sample.general.closestrefseqgenus = 'NA'
            else:
                # Populate attributes with 'NA'
                sample[self.analysistype].sortedgenera = 'NA'
                sample[self.analysistype].genus = 'NA'
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)
811437
def samples_to_records(samples, default_keys=None):
    
    from bcbio.pipeline import run_info
    RECORD_CONVERT_TO_LIST = set(["config__algorithm__tools_on", "config__algorithm__tools_off",
                                  "reference__genome_context"])
    all_keys = _get_all_cwlkeys(samples, default_keys)
    out = []
    for data in samples:
        for raw_key in sorted(list(all_keys)):
            key = raw_key.split("__")
            if tz.get_in(key, data) is None:
                data = tz.update_in(data, key, lambda x: None)
            if raw_key not in data["cwl_keys"]:
                data["cwl_keys"].append(raw_key)
            if raw_key in RECORD_CONVERT_TO_LIST:
                val = tz.get_in(key, data)
                if not val: val = []
                elif not isinstance(val, (list, tuple)): val = [val]
                data = tz.update_in(data, key, lambda x: val)
            # Booleans are problematic for CWL serialization, convert into string representation
            if isinstance(tz.get_in(key, data), bool):
                data = tz.update_in(data, key, lambda x: str(tz.get_in(key, data)))
        data["metadata"] = run_info.add_metadata_defaults(data.get("metadata", {}))
        out.append(data)
    return out