find consensus base based on nucleotide
    frequencies
774988
def find_consensus(bases):
    
    nucs = ['A', 'T', 'G', 'C', 'N']
    total = sum([bases[nuc] for nuc in nucs if nuc in bases])
    # save most common base as consensus (random nuc if there is a tie)
    try:
        top = max([bases[nuc] for nuc in nucs if nuc in bases])
    except:
        bases['consensus'] = ('N', 'n/a')
        bases['consensus frequency'] = 'n/a'
        bases['reference frequency'] = 'n/a'
        return bases
    top = [(nuc, bases[nuc]) for nuc in bases if bases[nuc] == top]
    if top[0][1] == 0:
        bases['consensus'] = ('n/a', 0)
    else:
        bases['consensus'] = random.choice(top)
    if total == 0:
        c_freq = 'n/a'
        ref_freq = 'n/a'
    else:
        c_freq = float(bases['consensus'][1]) / float(total)
        if bases['ref'] not in bases:
            ref_freq = 0
        else:
            ref_freq = float(bases[bases['ref']]) / float(total)
    bases['consensus frequency'] = c_freq
    bases['reference frequency'] = ref_freq
    return bases
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
716984
def _rngs(self):
      """"""
      outputs = []
      if len(self._exon_groups)==1:
         return [self._exon_groups.consensus('single')]
      z = 0 #output count
      begin = 0
      meeting_criteria = [i for i in range(0,len(self._exon_groups)) if  len(self._exon_groups) >= self._evidence]
      if len(meeting_criteria) == 0: return []
      finish = len(meeting_criteria)
      if len(meeting_criteria) > 0:
         begin = meeting_criteria[0]
         finish = meeting_criteria[-1]
      for i in range(0,len(self._exon_groups)):
         if z == begin:
            outputs.append(self._exon_groups[i].consensus('leftmost'))
         elif z == finish:
            #len(self._exon_groups)-1:
            outputs.append(self._exon_groups[i].consensus('rightmost'))
         else:
            outputs.append(self._exon_groups[i].consensus('internal'))
         z += 1
      v = [seqtools.structure.transcript.Exon(x) for x in outputs]
      v[0].set_leftmost()
      v[-1].set_rightmost()
      return v
1203448
private static SubstitutionMatrix<NucleotideCompound> getNucleotideMatrix(String file) {
		if (!nucleotideMatrices.containsKey(file)) {
			nucleotideMatrices.put(file, new SimpleSubstitutionMatrix<NucleotideCompound>(
					AmbiguityDNACompoundSet.getDNACompoundSet(), getReader(file), file));
		}
		return nucleotideMatrices.get(file);
	}
452393
def novo2(args):
    
    p = OptionParser(novo2.__doc__)
    p.set_fastq_names()
    p.set_align(pctid=95)
    opts, args = p.parse_args(args)

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

    trimmed, pf = args
    pctid = opts.pctid
    reads, samples = scan_read_files(trimmed, opts.names)

    # Set up directory structure
    clustdir = "uclust"
    acdir = "allele_counts"
    for d in (clustdir, acdir):
        mkdir(d)

    mm = MakeManager()
    clustfiles = []
    # Step 0 - clustering within sample
    for s in samples:
        flist = [x for x in reads if op.basename(x).split(".")[0] == s]
        outfile = s + ".P{0}.clustS".format(pctid)
        outfile = op.join(clustdir, outfile)
        cmd = "python -m jcvi.apps.uclust cluster --cpus=8"
        cmd += " {0} {1}".format(s, " ".join(flist))
        cmd += " --outdir={0}".format(clustdir)
        cmd += " --pctid={0}".format(pctid)
        mm.add(flist, outfile, cmd)
        clustfiles.append(outfile)

    # Step 1 - make consensus within sample
    allcons = []
    for s, clustfile in zip(samples, clustfiles):
        outfile = s + ".P{0}.consensus".format(pctid)
        outfile = op.join(clustdir, outfile)
        cmd = "python -m jcvi.apps.uclust consensus"
        cmd += " {0}".format(clustfile)
        mm.add(clustfile, outfile, cmd)
        allcons.append(outfile)

    # Step 2 - clustering across samples
    clustSfile = pf + ".P{0}.clustS".format(pctid)
    cmd = "python -m jcvi.apps.uclust mcluster {0}".format(" ".join(allcons))
    cmd += " --prefix={0}".format(pf)
    mm.add(allcons, clustSfile, cmd)

    # Step 3 - make consensus across samples
    locifile = pf + ".P{0}.loci".format(pctid)
    cmd = "python -m jcvi.apps.uclust mconsensus {0}".format(" ".join(allcons))
    cmd += " --prefix={0}".format(pf)
    mm.add(allcons + [clustSfile], locifile, cmd)

    mm.write()
1171181
def _normalize_seqs(s, t):
    """"""
    for ids in s:
        obj = s[ids]
        [obj.norm_freq.update({sample: 1.0 * obj.freq[sample] / (t[sample]+1) * 1000000}) for sample in obj.norm_freq]
        s[ids] = obj
    return s
735422
def _validate_residue_types(self):
        ''''''

        for chain_id, sequence_map in self.rosetta_to_atom_sequence_maps.iteritems():
            rosetta_sequence = self.rosetta_sequences[chain_id]
            atom_sequence = self.atom_sequences[chain_id]
            for rosetta_id, atom_id, _ in sequence_map:
                assert(rosetta_sequence[rosetta_id].ResidueAA == atom_sequence[atom_id].ResidueAA)

        for chain_id, sequence_map in self.atom_to_seqres_sequence_maps.iteritems():
            atom_sequence = self.atom_sequences[chain_id]
            seqres_sequence = self.seqres_sequences[chain_id]
            for atom_id, seqres_id, _ in sorted(sequence_map):
                assert(atom_sequence[atom_id].ResidueAA == seqres_sequence[seqres_id].ResidueAA)

        for chain_id, sequence_map in self.seqres_to_uniparc_sequence_maps.iteritems():
            if self.pdb_chain_to_uniparc_chain_mapping.get(chain_id):
                seqres_sequence = self.seqres_sequences[chain_id]
                uniparc_sequence = self.uniparc_sequences[self.pdb_chain_to_uniparc_chain_mapping[chain_id]]
                for seqres_id, uniparc_id_resid_pair, substitution_match in sequence_map:
                    uniparc_id = uniparc_id_resid_pair[1]
                    # Some of the matches may not be identical but all the '*' Clustal Omega matches should be identical
                    if substitution_match and substitution_match.clustal == 1:
                        assert(seqres_sequence[seqres_id].ResidueAA == uniparc_sequence[uniparc_id].ResidueAA)
1154157
def periodicity(freq_or_frame):
    
    if hasattr(freq_or_frame, 'rule_code'):
        rc = freq_or_frame.rule_code
        rc = rc.split('-')[0]
        factor = PER_YEAR_MAP.get(rc, None)
        if factor is not None:
            return factor / abs(freq_or_frame.n)
        else:
            raise Exception('Failed to determine periodicity. No factor mapping for %s' % freq_or_frame)
    elif isinstance(freq_or_frame, basestring):
        factor = PER_YEAR_MAP.get(freq_or_frame, None)
        if factor is not None:
            return factor
        else:
            raise Exception('Failed to determine periodicity. No factor mapping for %s' % freq_or_frame)
    elif isinstance(freq_or_frame, (pd.Series, pd.DataFrame, pd.TimeSeries)):
        freq = freq_or_frame.index.freq
        if not freq:
            freq = pd.infer_freq(freq_or_frame.index)
            if freq:
                return periodicity(freq)
            else:
                # Attempt to resolve it
                import warnings

                freq = guess_freq(freq_or_frame.index)
                warnings.warn('frequency not set. guessed it to be %s' % freq)
                return periodicity(freq)
        else:
            return periodicity(freq)
    else:
        raise ValueError("periodicity expects DataFrame, Series, or rule_code property")
12374
def consensus(args):
    
    p = OptionParser(consensus.__doc__)
    p.add_option("--fasta", default=False, action="store_true",
            help="Generate consensus FASTA sequences [default: %default]")
    p.add_option("--mask", default=0, type="int",
            help="Mask bases with quality lower than")
    opts, args = p.parse_args(args)

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

    fastafile, bamfile = args
    fasta = opts.fasta
    suffix = "fasta" if fasta else "fastq"
    pf = bamfile.rsplit(".", 1)[0]
    cnsfile = pf + ".cns.{0}".format(suffix)
    vcfgzfile = pf + ".vcf.gz"
    vcf([fastafile, bamfile, "-o", vcfgzfile])
    cmd += "zcat {0} | vcfutils.pl vcf2fq".format(vcfgzfile)
    if fasta:
        cmd += " | seqtk seq -q {0} -A -".format(opts.mask)

    sh(cmd, outfile=cnsfile)
515491
def find_coordinates(hmms, bit_thresh):
    
    # get coordinates from cmsearch output
    seq2hmm = parse_hmm(hmms, bit_thresh)
    seq2hmm = best_model(seq2hmm)
    group2hmm = {} # group2hmm[seq][group] = [model, strand, coordinates, matches, gaps]
    for seq, info in list(seq2hmm.items()):
        group2hmm[seq] = {}
        # info = [model, [[hit1], [hit2], ...]]
        for group_num, group in enumerate(hit_groups(info[1])):
            # group is a group of hits to a single 16S gene
            # determine matching strand based on best hit
            best = sorted(group, reverse = True, key = itemgetter(-1))[0]
            strand = best[5]
            coordinates = [i[0] for i in group] + [i[1] for i in group]
            coordinates = [min(coordinates), max(coordinates), strand]
            # make sure all hits are to the same strand
            matches = [i for i in group if i[5] == strand]
            # gaps = [[gstart, gend], [gstart2, gend2]]
            gaps = check_gaps(matches)
            group2hmm[seq][group_num] = [info[0], strand, coordinates, matches, gaps]
    return group2hmm
216779
def filter_clades(self):
        ""
        passed = []
        carrs = np.array([list(i[0]) for i in self.clade_counts], dtype=int)
        freqs = np.array([i[1] for i in self.clade_counts])

        for idx in range(carrs.shape[0]):
            conflict = False
            if freqs[idx] < self.cutoff:
                continue
            
            for pidx in passed:
                intersect = np.max(carrs[idx] + carrs[pidx]) > 1
                
                # is either one a subset of the other?
                subset_test0 = np.all(carrs[idx] - carrs[pidx] >= 0)
                subset_test1 = np.all(carrs[pidx] - carrs[idx] >= 0)
                if intersect:
                    if (not subset_test0) and (not subset_test1):
                        conflict = True

            if not conflict:
                passed.append(idx)

        rclades = []
        for idx in passed:
            rclades.append((carrs[idx], freqs[idx]))
        self.fclade_counts = rclades
222479
def generate_complementary_strand(strand1):
        """"""
        rise_adjust = (
            strand1.rise_per_nucleotide * strand1.axis.unit_tangent) * 2
        strand2 = NucleicAcidStrand.from_start_and_end(
            strand1.helix_end - rise_adjust, strand1.helix_start - rise_adjust,
            generate_antisense_sequence(strand1.base_sequence),
            phos_3_prime=strand1.phos_3_prime)
        ad_ang = dihedral(strand1[0]["C1'"]._vector, strand1.axis.start,
                          strand2.axis.start + rise_adjust,
                          strand2[-1]["C1'"]._vector)
        strand2.rotate(
            225.0 + ad_ang, strand2.axis.unit_tangent,
            point=strand2.helix_start)  # 225 is the base adjust
        return strand2
883344
def build(self):
        """"""
        ang_per_res = (2 * numpy.pi) / self.nucleotides_per_turn
        atom_offset_coords = _backbone_properties[self.helix_type]['atoms']
        if self.handedness == 'l':
            handedness = -1
        else:
            handedness = 1
        base_atom_labels = _backbone_properties[self.helix_type]['labels']
        monomers = []
        mol_code_format = _backbone_properties[self.helix_type]['mol_code_format']
        for i, b in enumerate(self.base_sequence):
            nucleotide = Nucleotide(
                mol_code=mol_code_format.format(b), ampal_parent=self)
            atoms_dict = OrderedDict()
            if (i == (len(self.base_sequence) - 1)) and not self.phos_3_prime:
                # Do not include phosphate on last nucleotide
                atom_labels = base_atom_labels[3:] + [_bases[b]['labels'][2]]
                atom_offsets = {k: v for k, v in zip(
                    atom_labels, atom_offset_coords[3:])}
            else:
                atom_labels = base_atom_labels + [_bases[b]['labels'][2]]
                atom_offsets = {k: v for k, v in zip(
                    atom_labels, atom_offset_coords)}
            for atom_label in atom_labels:
                r, zeta, z_shift = atom_offsets[atom_label]
                rot_ang = ((i * ang_per_res) + zeta) * handedness
                z = (self.rise_per_nucleotide * i) + z_shift
                coords = cylindrical_to_cartesian(
                    radius=r, azimuth=rot_ang, z=z, radians=True)
                atom = Atom(
                    coordinates=coords, element=atom_label[0],
                    ampal_parent=nucleotide, res_label=atom_label)
                atoms_dict[atom_label] = atom
            base_ref = _bases[b]['ref_atom']
            rot_adj = _bases[b]['rot_adj']
            base_dict = OrderedDict(
                zip(_bases[b]['labels'], _bases[b]['atoms']))
            translation, angle, axis, point = find_transformations(
                base_dict[base_ref], base_dict["C1'"],
                atoms_dict[base_ref]._vector, atoms_dict["C1'"]._vector)
            q1 = Quaternion.angle_and_axis(angle, axis)
            # Align N9 C1'
            for k, v in base_dict.items():
                base_dict[k] = q1.rotate_vector(v, point) + translation
            # Rotate to align O4'
            axis = numpy.array(base_dict["C1'"]) - base_dict[base_ref]
            angle = dihedral(base_dict["O4'"], base_dict[base_ref],
                             base_dict["C1'"], atoms_dict["O4'"]) - rot_adj
            q2 = Quaternion.angle_and_axis(angle, axis)
            for k, v in list(base_dict.items()):
                if k not in atoms_dict:
                    atom = Atom(q2.rotate_vector(v, base_dict[base_ref]),
                                element=k[0], ampal_parent=nucleotide,
                                res_label=k)
                    atoms_dict[k] = atom
            nucleotide.atoms = atoms_dict
            monomers.append(nucleotide)
        self._monomers = monomers
        self.relabel_monomers()
        self.relabel_atoms()
        return
1221766
def collapse_to_consensus(seqrecords, strict=False, do_iron=True):
    
    level = 0
    name = seqrecords[0].id
    # If this is a CMA alignment, extract additional info:
    if hasattr(seqrecords, '_records'):
        if hasattr(seqrecords, 'level'):
            level = seqrecords.level
        if hasattr(seqrecords, 'name'):
            name = seqrecords.name
        seqrecords = seqrecords._records

    consensus = seqrecords.pop(0)
    cons_length = len(consensus)
    for i, s in enumerate(seqrecords):
        if len(s) != cons_length:
            raise ValueError(
            "Sequence #%d has length %d, consensus is %d"
            % (i+2, len(s), cons_length))

    if '.' in str(consensus.seq):
        # Strict -- error if there's a '-'
        if '-' in str(consensus.seq):
            if strict:
                raise ValueError("Consensus contains '-' gap characters")
            logging.warn("Consensus sequence contains both '.' and '-' gap "
                         "characters -- is it really the consensus?")
            aligned_cols = [(c not in '.-') for c in str(consensus.seq)]
        else:
            aligned_cols = [c != '.' for c in str(consensus.seq)]
    else:
        # A little more ambiguous...
        aligned_cols = [c != '-' for c in str(consensus.seq)]

    consensus.seq = replace_asterisks(consensus.seq, 'consensus')
    # Start a block with the consensus sequence
    block = consensus2block(consensus, level=level, name=name)
    qlen = block['query_length']

    # Collapse & add remaining sequences to the block
    for index, rec in zip(xrange(2, len(seqrecords)+2), seqrecords):
        # Collapse rec.seq down to aligned size
        new_mol_seq = []
        is_beginning = True
        for aligned_col, char in zip(aligned_cols,
                                     replace_asterisks(rec.seq, index)):
            if aligned_col:
                is_beginning = False
                if char in '-.':
                    # deletion
                    new_mol_seq.append('-')
                else:
                    # aligned character
                    new_mol_seq.append(char.upper())
            else:
                # it's an insert or nothing
                # (also, skip any left-side inserts)
                if char not in '-.' and not is_beginning:
                    new_mol_seq.append(char.lower())
        rec.seq = ''.join(new_mol_seq)
        if do_iron:
            rec.seq = iron(rec.seq)
        block['sequences'].append(seqrecord2sequence(rec, qlen, index))

    return block
387768
def get_frequency_shift(
            self,
            grid_points,
            temperatures=np.arange(0, 1001, 10, dtype='double'),
            epsilons=None,
            output_filename=None):
        

        if self._interaction is None:
            self.set_phph_interaction()
        if epsilons is None:
            _epsilons = [0.1]
        else:
            _epsilons = epsilons
        self._grid_points = grid_points
        get_frequency_shift(self._interaction,
                            self._grid_points,
                            self._band_indices,
                            _epsilons,
                            temperatures,
                            output_filename=output_filename,
                            log_level=self._log_level)
144171
def col_frequencies(col, weights=None, gap_chars='-.'):
    """"""
    counts = col_counts(col, weights, gap_chars)
    # Reduce to frequencies
    scale = 1.0 / sum(counts.values())
    return dict((aa, cnt * scale) for aa, cnt in counts.iteritems())
1149178
def build_bam_tags():
    ''''''
    #pylint: disable=unused-argument
    def _combine_filters(fam, paired_align, align):
        filters = [x.filter_value for x in [fam, align] if x and x.filter_value]
        if filters:
            return ";".join(filters).replace('; ', ';')
        return None

    boolean_tag_value = {True:1}
    tags = [
        BamTag("X0", "Z",
               ("filter (why the alignment was excluded)"),
               _combine_filters),
        BamTag("X1", "Z",
               ("leftmost~rightmost matched pair positions"),
               lambda fam, pair, align: pair.positions('{left}~{right}')),
        BamTag("X2", "Z",
               ("L~R CIGARs"),
               lambda fam, pair, align: pair.cigars('{left}~{right}')),
        BamTag("X3", "i",
               "unique identifier for this alignment family",
               lambda fam, pair, align: fam.umi_sequence),
        BamTag("X4", "Z",
               ("L~R UMT barcodes for this alignment family; because "
                "of fuzzy matching the family UMT may be distinct "
                "from the UMT of the original alignment"),
               lambda fam, pair, align: fam.umt('{left}~{right}')),
        BamTag("X5", "i",
               "family size (number of align pairs in this family)",
               lambda fam, pair, align: fam.included_pair_count),
        BamTag("X6", "i",
               ("presence of this tag signals that this alignment "
                "would be the template for the consensus alignment"),
               lambda fam, pair, align: boolean_tag_value.get(fam.is_consensus_template(align), None))]
    return tags