Read read length and genome size from the config and calculate
        the approximate number of counts (or base pairs) in 1x of depth
511218
def _get_counts_in_1x(data_is_basepairs):
    
    read_length = float(getattr(config, 'preseq', {}).get('read_length', 0))
    genome_size = getattr(config, 'preseq', {}).get('genome_size')
    if genome_size:
        try:
            genome_size = float(genome_size)
        except ValueError:
            presets = {'hg19_genome': 2897310462,
                       'hg38_genome': 3049315783,
                       'mm10_genome': 2652783500}
            if genome_size in presets:
                genome_size = presets[genome_size]
            else:
                log.warn('The size for genome ' + genome_size + ' is unknown to MultiQC, ' +
                         'please specify it explicitly or choose one of the following: ' +
                         ', '.join(presets.keys()) + '. Falling back to molecule counts.')
                genome_size = None

    if genome_size:
        if data_is_basepairs:
            return genome_size
        elif read_length:
            return genome_size / read_length
    else:
        return None
314743
def reads(args):
    
    p = OptionParser(reads.__doc__)
    p.add_option("-p", dest="prefix_length", default=4, type="int",
            help="group the reads based on the first N chars [default: %default]")
    opts, args = p.parse_args(args)

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

    frgscffile, = args
    prefix_length = opts.prefix_length

    fp = open(frgscffile)
    keyfn = lambda: defaultdict(int)
    counts = defaultdict(keyfn)
    for row in fp:
        f = FrgScfLine(row)
        fi = f.fragmentID[:prefix_length]
        counts[f.scaffoldID][fi] += 1

    for scf, count in sorted(counts.items()):
        print("{0}\t{1}".format(scf,
                ", ".join("{0}:{1}".format(*x) for x in sorted(count.items()))))
186677
def normalise_reads(self):
        
        logging.info('Normalising reads to a kmer depth of 100')
        for sample in self.metadata:
            # Set the name of the normalised read files
            sample.general.normalisedreads = [fastq.split('.fastq.gz')[0] + '_normalised.fastq.gz'
                                              for fastq in sorted(sample.general.fastqfiles)]
            try:
                # Run the normalisation command
                out, err, cmd = bbtools.bbnorm(forward_in=sorted(sample.general.trimmedcorrectedfastqfiles)[0],
                                               forward_out=sample.general.normalisedreads[0],
                                               returncmd=True,
                                               threads=self.cpus)
                sample[self.analysistype].normalisecmd = cmd
                write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
            except CalledProcessError:
                sample.general.normalisedreads = sample.general.trimmedfastqfiles
            except IndexError:
                sample.general.normalisedreads = list()
446800
def overlap_size(self,in_genomic_range):
    
    if self.chr != in_genomic_range.chr:
      return 0
    if self.end < in_genomic_range.start:
      return 0
    if in_genomic_range.end < self.start:
      return 0
    if self.start > in_genomic_range.end:
      return 0
    if self.start >= in_genomic_range.start and self.end <= in_genomic_range.end:
      return self.end-self.start+1
    if self.start <= in_genomic_range.start and self.end >= in_genomic_range.end:
      return in_genomic_range.end-in_genomic_range.start+1
    if self.start <= in_genomic_range.start and self.end >= in_genomic_range.start:
      return self.end-in_genomic_range.start+1
    if self.start <= in_genomic_range.end and self.end >= in_genomic_range.end:
      return in_genomic_range.end-self.start+1
    if in_genomic_range.start <= self.start and in_genomic_range.end >= self.start:
      return in_genomic_range.end-self.start+1
    if in_genomic_range.start <= self.end and in_genomic_range.end >= self.end:
      return self.end-in_genomic_range.start+1
    sys.stderr.write("overlap_size: unprogrammed error\n")
    return 0
242093
def len_on_depth(d, depth):
        
        counter = 0
        for node in DictTree.v_depth(d, depth-1):
            counter += DictTree.length(node)
        return counter
466433
def load(self, read_tuple_name):
        
        self.prefix_width = 0
        self.read_tuple_id_width = 0
        self.genome_id_width = 0
        self.chr_id_width = 0
        self.coor_width = 0

        parts = read_tuple_name.split("__")
        self.prefix_width = len(parts[0])
        self.read_tuple_id_width = len(parts[1])

        segments = parts[2][1:-1].split("),(")
        for segment in segments:
            int_widths = list(map(len, segment.split(",")))
            self.genome_id_width = max(self.genome_id_width, int_widths[0])
            self.chr_id_width = max(self.chr_id_width, int_widths[1])
            self.coor_width = max(self.coor_width, int_widths[2], int_widths[3])
41894
def depth(args):
    
    from jcvi.utils.range import range_depth

    p = OptionParser(depth.__doc__)
    p.add_option("--depthfile",
                 help="Generate file with gene and depth [default: %default]")
    p.add_option("--histogram", default=False, action="store_true",
                 help="Plot histograms in PDF")
    p.add_option("--xmax", type="int", help="x-axis maximum to display in plot")
    p.add_option("--title", default=None, help="Title to display in plot")
    p.add_option("--quota", help="Force to use this quota, e.g. 1:1, 1:2 ...")
    p.set_beds()

    opts, args = p.parse_args(args)

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

    anchorfile, = args
    qbed, sbed, qorder, sorder, is_self = check_beds(anchorfile, p, opts)
    depthfile = opts.depthfile
    ac = AnchorFile(anchorfile)
    qranges = []
    sranges = []
    blocks = ac.blocks
    for ib in blocks:
        q, s, t = zip(*ib)
        q = [qorder[x] for x in q]
        s = [sorder[x] for x in s]
        qrange = (min(q)[0], max(q)[0])
        srange = (min(s)[0], max(s)[0])
        qranges.append(qrange)
        sranges.append(srange)
        if is_self:
            qranges.append(srange)

    qgenome = op.basename(qbed.filename).split(".")[0]
    sgenome = op.basename(sbed.filename).split(".")[0]
    qtag = "Genome {0} depths".format(qgenome)
    print("{}:".format(qtag), file=sys.stderr)
    dsq, details = range_depth(qranges, len(qbed))
    if depthfile:
        fw = open(depthfile, "w")
        write_details(fw, details, qbed)

    if is_self:
        return

    stag = "Genome {0} depths".format(sgenome)
    print("{}:".format(stag), file=sys.stderr)
    dss, details = range_depth(sranges, len(sbed))
    if depthfile:
        write_details(fw, details, sbed)
        fw.close()
        logging.debug("Depth written to `{0}`.".format(depthfile))

    if not opts.histogram:
        return

    from jcvi.graphics.base import plt, quickplot_ax, savefig, normalize_axes

    # Plot two histograms one for query genome, one for subject genome
    plt.figure(1, (6, 3))
    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)

    xmax = opts.xmax or max(4, max(dsq.keys() + dss.keys()))
    if opts.quota:
        speak, qpeak = opts.quota.split(":")
        qpeak, speak = int(qpeak), int(speak)
    else:
        qpeak = find_peak(dsq)
        speak = find_peak(dss)

    qtag = "# of {} blocks per {} gene".format(sgenome, qgenome)
    stag = "# of {} blocks per {} gene".format(qgenome, sgenome)
    quickplot_ax(ax1, dss, 0, xmax, stag, ylabel="Percentage of genome",
                 highlight=range(1, speak + 1))
    quickplot_ax(ax2, dsq, 0, xmax, qtag, ylabel=None,
                 highlight=range(1, qpeak + 1))

    title = opts.title or "{} vs {} syntenic depths\n{}:{} pattern"\
                    .format(qgenome, sgenome, speak, qpeak)
    root = f.add_axes([0, 0, 1, 1])
    vs, pattern = title.split('\n')
    root.text(.5, .97, vs, ha="center", va="center", color="darkslategray")
    root.text(.5, .925, pattern, ha="center", va="center",
                                 color="tomato", size=16)
    print(title, file=sys.stderr)

    normalize_axes(root)

    pf = anchorfile.rsplit(".", 1)[0] + ".depth"
    image_name = pf + ".pdf"
    savefig(image_name)
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
179872
def find_n50(self):
        
        for sample in self.metadata:
            # Initialise the N50 attribute in case there is no assembly, and the attribute is not created in the loop
            sample[self.analysistype].n50 = '-'
            # Initialise a variable to store a running total of contig lengths
            currentlength = 0
            for contig_length in sample[self.analysistype].contig_lengths:
                # Increment the current length with the length of the current contig
                currentlength += contig_length
                # If the current length is now greater than the total genome / 2, the current contig length is the N50
                if currentlength >= sample[self.analysistype].genome_length * 0.5:
                    # Populate the dictionary, and break the loop
                    sample[self.analysistype].n50 = contig_length
                    break
1230699
def g_to_n(self, g_interval):
        """"""

        # This code is extremely convoluted. To begin with, it
        # confuses interbase intervals for a *single* hgvs position
        # with intervals of *two* hgvs positions.  This is the origin
        # of bug #205 (now closed). Although the bug is itself fixed,
        # the mappers really need an overhaul. See #208.

        def _hgvs_offset(g_position, grs, gre, strand):
            """ Calculates the HGVS coordinate offset from a given genomic position"""
            if g_position == grs or g_position == gre:
                return 0
            mid = (grs + gre) / 2
            if g_position < mid or (g_position == mid and strand == 1):
                offset = g_position - grs
            else:
                offset = g_position - gre
            return offset

        def map_g_to_n_pos(pos):
            g_ci = _hgvs_coord_to_ci(pos, pos)
            # frs, fre = (f)orward (r)na (s)tart & (e)nd; forward w.r.t. genome
            frs, fre = self.im.map_ref_to_tgt(
                g_ci[0] - self.gc_offset, g_ci[1] - self.gc_offset, max_extent=False)
            if self.strand == -1:
                frs, fre = self.tgt_len - fre, self.tgt_len - frs
            # get BaseOffsetPosition information for n.
            if self.strand == 1:
                grs, gre = self.im.map_tgt_to_ref(frs, fre, max_extent=False)
            elif self.strand == -1:
                grs, gre = self.im.map_tgt_to_ref((self.tgt_len - fre), (self.tgt_len - frs),
                                                  max_extent=False)
            grs, gre = grs + self.gc_offset, gre + self.gc_offset
            # 0-width interval indicates an intron.  Need to calculate offsets but we"re are in ci coordinates
            # requires adding 1 strategically to get the HGVS position
            # (shift coordinates to 3' end of the ref nucleotide)

            if frs == fre:
                start_offset = _hgvs_offset(g_ci[0] + 1, grs, gre + 1, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre + 1, self.strand)
            else:
                start_offset = _hgvs_offset(g_ci[0], grs, gre, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre, self.strand)
            if self.strand == -1:
                start_offset, end_offset = self.strand * end_offset, self.strand * start_offset
            if start_offset > 0:
                frs -= 1
            if end_offset < 0:
                fre += 1
            start_base = _ci_to_hgvs_coord(frs, fre)[0]
            end_base = _ci_to_hgvs_coord(frs, fre)[1]
            return start_base, start_offset, end_base, end_offset

        start_bo = map_g_to_n_pos(g_interval.start.base)[0:2]
        end_bo = start_bo if g_interval.start.base == g_interval.end.base else map_g_to_n_pos(
            g_interval.end.base)[2:4]
        if self.strand == -1:
            start_bo, end_bo = end_bo, start_bo

        return hgvs.location.BaseOffsetInterval(
            start=hgvs.location.BaseOffsetPosition(
                base=start_bo[0], offset=start_bo[1], datum=Datum.SEQ_START),
            end=hgvs.location.BaseOffsetPosition(
                base=end_bo[0], offset=end_bo[1], datum=Datum.SEQ_START),
            uncertain=g_interval.uncertain)
1197641
def calc_frequencies(genomes, bp_table, min_cov, min_per):
    
    nucs = ['A', 'T', 'G', 'C', 'N']
    if bp_table is not False:
        bp_table = open(bp_table, 'w')
        header = ['#genome', 'contig', 'sample', 'position', \
                    'reference', 'ref. frequency', \
                    'consensus', 'con. frequency', \
                    'A', 'T', 'G', 'C', 'N', '# insertions', '# deletions']
        print('\t'.join(header), file=bp_table)
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample, stats in list(samples.items()):
                for pos, ps in enumerate(stats['bp_stats'], 1):
                    coverage = sum([ps[nuc] for nuc in nucs])
                    for nuc in nucs:
                        # make sure support for base passes thresholds
                        nuc_cov = ps[nuc]
                        if coverage == 0:
                            nuc_per = 0
                        else:
                            nuc_per = (float(nuc_cov)/coverage)*100
                        if nuc_cov < min_cov or nuc_per < min_per:
                            del ps[nuc]
                    ps = find_consensus(ps)
                    genomes[genome][contig][sample][pos] = ps
                    if bp_table is not False:
                        out = [genome, contig, sample, pos]
                        for i in ['ref', 'reference frequency', \
                                  'consensus', 'consensus frequency', \
                                  'A', 'T', 'G', 'C', 'N', \
                                  'In', 'Del']:
                            try:
                                if i == 'consensus':
                                    out.append(ps[i][0])
                                elif i in ['In', 'Del']:
                                    out.append(len(ps[i]))
                                else:
                                    out.append(ps[i])
                            except:
                                out.append('n/a')
                        print('\t'.join([str(i) for i in out]), file=bp_table)
    return genomes
478777
def calculateProbableRootOfGeneTree(speciesTree, geneTree, processID=lambda x : x):
    
    #get all rooted trees
    #run dup calc on each tree
    #return tree with fewest number of dups
    if geneTree.traversalID.midEnd <= 3:
        return (0, 0, geneTree)
    checkGeneTreeMatchesSpeciesTree(speciesTree, geneTree, processID)
    l = []
    def fn(tree):
        if tree.traversalID.mid != geneTree.left.traversalID.mid and tree.traversalID.mid != geneTree.right.traversalID.mid:
            newGeneTree = moveRoot(geneTree, tree.traversalID.mid)
            binaryTree_depthFirstNumbers(newGeneTree)
            dupCount, lossCount = calculateDupsAndLossesByReconcilingTrees(speciesTree, newGeneTree, processID)
            l.append((dupCount, lossCount, newGeneTree))
        if tree.internal:
            fn(tree.left)
            fn(tree.right)
    fn(geneTree)
    l.sort()
    return l[0][2], l[0][0], l[0][1]
677102
def scan(self, seqs, nreport=100, scan_rc=True, normalize=False):
        

        if not self.threshold:
            sys.stderr.write(
                "Using default threshold of 0.95. "
                "This is likely not optimal!\n"
                )
            self.set_threshold(threshold=0.95)

        seqs = as_fasta(seqs, genome=self.genome)
           
        it = self._scan_sequences(seqs.seqs, 
                    nreport, scan_rc)
       
        if normalize:
            if len(self.meanstd) == 0:
                self.set_meanstd()
            mean_std = [self.meanstd.get(m_id) for m_id in self.motif_ids]
            means = [x[0] for x in  mean_std]
            stds = [x[1] for x in  mean_std]


        for result in it:
            if normalize:
                zresult = [] 
                for i,mrow in enumerate(result):
                    mrow = [((x[0] - means[i]) / stds[i], x[1], x[2]) for x in mrow]
                    zresult.append(mrow)
                yield zresult
            else:
                yield result
637636
def _check_reads_hit(self, alignment_io, min_aligned_fraction):
        ''''''
        to_return = []
        alignment_length = None
        for s in SeqIO.parse(alignment_io, "fasta"):
            if not alignment_length:
                alignment_length = len(s.seq)
                min_length = int(min_aligned_fraction * alignment_length)
                logging.debug("Determined min number of aligned bases to be %s" % min_length)
            elif len(s.seq) != alignment_length:
                raise Exception("Alignment file appears to not be of uniform length")

            num_unaligned = s.seq.count('-')
            num_aligned = alignment_length-num_unaligned
            logging.debug("Sequence %s has %d aligned positions" % (s.name, alignment_length-num_unaligned))
            if num_aligned <= min_length:
                to_return.append(s.name)
        return to_return
1273208
def create_genome_size_dict(genome):
    # type: (str) -> Dict[str,int]
    """"""

    size_file = get_genome_size_file(genome)
    size_lines = open(size_file).readlines()

    size_dict = {}
    for line in size_lines:
        genome, length = line.split()
        size_dict[genome] = int(length)

    return size_dict
1169418
def _generate_reads(seq, name):
    """"""
    reads = dict()
    if len(seq) < 130 and len(seq) > 70:
        reads.update(_mature(seq[:40], 0, name))
        reads.update(_mature(seq[-40:], len(seq) - 40, name))
        reads.update(_noise(seq, name))
        reads.update(_noise(seq, name, 25))
    return reads
1141982
def count(args):
    
    from jcvi.graphics.histogram import stem_leaf_plot
    from jcvi.utils.cbook import SummaryStats

    p = OptionParser(count.__doc__)
    p.add_option("--csv", help="Write depth per contig to file")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    csv = open(opts.csv, "w") if opts.csv else None

    f = Fasta(fastafile, lazy=True)
    sizes = []
    for desc, rec in f.iterdescriptions_ordered():
        if desc.startswith("singleton"):
            sizes.append(1)
            continue

        # consensus_for_cluster_0 with 63 sequences
        if "with" in desc:
            name, w, size, seqs = desc.split()
            if csv:
                print("\t".join(str(x)
                                for x in (name, size, len(rec))), file=csv)
            assert w == "with"
            sizes.append(int(size))
        # MRD85:00603:02472;size=167;
        else:
            name, size, tail = desc.split(";")
            sizes.append(int(size.replace("size=", "")))

    if csv:
        csv.close()
        logging.debug("File written to `{0}`".format(opts.csv))

    s = SummaryStats(sizes)
    print(s, file=sys.stderr)
    stem_leaf_plot(s.data, 0, 100, 20, title="Cluster size")