Creates the self.seqpos2genome dictionary that converts positions
        relative to the sequence to genome coordinates.
321299
def init_genome_coordinates(self) :
        """"""
        self.seqpos2genome = {}

        # record genome positions for each sequence position
        seq_pos = 0
        for estart, eend in self.exons:
            for genome_pos in range(estart, eend):
                if self.strand == '+':
                    self.seqpos2genome[seq_pos] = genome_pos
                elif self.strand == '-':
                    tmp = self.cds_len - seq_pos - 1
                    self.seqpos2genome[tmp] = genome_pos
                seq_pos += 1

        # recode 5' splice site locations
        for i in range(0, self.five_ss_len):
            seq_pos = self.cds_len + i
            ss_ix = i // 2  # the ss_ix'th 5'ss starting from upstream tx
            pos_in_ss = i % 2  # whether first/second nuc in splice site

            # determine genome coordinates for 5' splice site
            if self.strand == '+':
                self.seqpos2genome[seq_pos] = self.exons[ss_ix][1] + pos_in_ss
            else:
                exon_pos = -1 - ss_ix
                self.seqpos2genome[seq_pos] = self.exons[exon_pos][0] - pos_in_ss - 1

        # recode 3' splice site locations
        for i in range(0, self.three_ss_len):
            seq_pos = self.cds_len + self.five_ss_len + i
            ss_ix = i // 2  # the ss_ix'th 3'ss starting from upstream tx
            pos_in_ss = i % 2  # whether first/second nuc in splice site

            # determine genome coordinates for 3' splice site
            if self.strand == '+':
                self.seqpos2genome[seq_pos] = self.exons[ss_ix+1][0] - 2 + pos_in_ss
            else:
                exon_pos = -1 - ss_ix
                self.seqpos2genome[seq_pos] = self.exons[exon_pos-1][1] + 1 - pos_in_ss
242529
def split_coords_2d(seq):
    
    lons, lats = [], []
    for i, el in enumerate(seq):
        if i % 2 == 0:
            lons.append(valid.longitude(el))
        elif i % 2 == 1:
            lats.append(valid.latitude(el))
    return list(zip(lons, lats))
1176321
def generate_antisense_sequence(sequence):
    """"""
    dna_antisense = {
        'A': 'T',
        'T': 'A',
        'C': 'G',
        'G': 'C'
    }
    antisense = [dna_antisense[x] for x in sequence[::-1]]
    return ''.join(antisense)
92958
function sequence ( parserSyntax ) {
				var i, res,
					originalPos = pos,
					result = [];

				for ( i = 0; i < parserSyntax.length; i++ ) {
					res = parserSyntax[i]();

					if ( res === null ) {
						pos = originalPos;

						return null;
					}

					result.push( res );
				}

				return result;
			}
604776
def create_features(bam_in, loci_file, reference, out_dir):
    
    lenvec_plus = op.join(out_dir, 'genomic_lenvec.plus')
    lenvec_minus = op.join(out_dir, 'genomic_lenvec.minus')
    compute_genomic_cmd = ("compute_genomic_lenvectors "
                           "{bam_in} {lenvec_plus} "
                           "{lenvec_minus} "
                           "{min_len} "
                           "{max_len} ")
    index_genomic_cmd = ("index_genomic_lenvectors "
                         "{lenvec} ")
    genomic_lenvec = op.join(out_dir, 'genomic_lenvec')
    feat_len_file = op.join(out_dir, 'feat_lengths.txt')
    compute_locus_cmd = ("compute_locus_lenvectors "
                         "{loci_file} "
                         "{genomic_lenvec} "
                         "{min_len} "
                         "{max_len} "
                         "> {feat_len_file}")
    cov_S_file = op.join(out_dir, 'loci.cov_anti')
    coverage_anti_cmd = ("coverageBed -S -counts -b "
                         "{bam_in} -a {loci_file} "
                         "> {cov_S_file}")
    feat_posentropy = op.join(out_dir, 'feat_posentropy.txt')
    entropy_cmd = ("compute_locus_entropy.rb "
                   "{counts_reads} "
                   "> {feat_posentropy}")
    with utils.chdir(out_dir):
        run(compute_genomic_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute_genomic")
        run(index_genomic_cmd.format(lenvec=lenvec_plus), "Run index in plus")
        run(index_genomic_cmd.format(lenvec=lenvec_minus), "Run index in minus")
        run(compute_locus_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute locus")
        run(coverage_anti_cmd.format(**locals()), "Run coverage antisense")
        feat_antisense = _order_antisense_column(cov_S_file, min_trimmed_read_len)

        counts_reads = _reads_per_position(bam_in, loci_file, out_dir)
        run(entropy_cmd.format(**locals()), "Run entropy")

        rnafold = calculate_structure(loci_file, reference)
863295
def calculate_oobatake_dS(seq, temp):
    

    seq = ssbio.protein.sequence.utils.cast_to_str(seq)

    dS = 0
    temp += 273.15
    T0 = 298.15
    dCp_sum = _sum_of_dCp(seq)
    for aa in seq:
        S0 = oobatake_dictionary[aa]['dS']
        dS += S0
    return dS + dCp_sum * math.log(temp / T0)
1318665
def save_refs(sam, fastas, genomes, s2b):
    
    if s2b is False:
        s2b = {}
    for fasta in fastas:
        # if no reference sequence supplied, get length from sam file
        if fasta is False:
            for line in sam:
                if line.startswith('@PG'):
                    break
                if line.startswith('@SQ') is False:
                    continue
                line = line.strip().split()
                contig, length = line[1].split(':', 1)[1], int(line[2].split(':', 1)[1])
                if contig not in s2b:
                    genome = 'n/a'
                    s2b[contig] = genome
                else:
                    genome = s2b[contig]
                if genome not in genomes:
                    genomes[genome] = {}
                if contig not in genomes[genome]:
                    genomes[genome][contig] = {}
                bp_stats = []
                for i in range(0, length):
                    bp_stats.append({'A':0, 'T':0, 'G':0, 'C':0, 'N':0, \
                                        'In':[], 'Del':[], 'ref':'N'})
                genomes[genome][contig][sam.name] = {'bp_stats':bp_stats}
                return genomes, s2b
        # save reference sequences, if available
        genome = fasta.rsplit('.', 1)[0]
        if genome not in genomes:
            genomes[genome] = {}
        bp_stats = []
        for seq in parse_fasta(fasta):
            contig = seq[0].split('>')[1].split()[0]
            s2b[contig] = genome
            if contig not in genomes[genome]:
                genomes[genome][contig] = {}
            bp_stats = []
            for base in seq[1]:
                bp_stats.append({'A':0, 'T':0, 'G':0, 'C':0, 'N':0, \
                                    'In':[], 'Del':[], 'ref':base.upper()})
            genomes[genome][contig][sam.name] = {'bp_stats':bp_stats}
        return genomes, s2b
1237691
def build_vcf_inversion(x1, x2, genome_2bit):
    
    id1 = "hydra{0}".format(x1.name)
    start_coords = sorted([x1.start1, x1.end1, x2.start1, x2.end1])
    end_coords = sorted([x1.start2, x1.end2, x2.start2, x2.start2])
    start_pos = (start_coords[1] + start_coords[2]) // 2
    end_pos = (end_coords[1] + end_coords[2]) // 2
    base1 = genome_2bit[x1.chrom1].get(start_pos, start_pos + 1).upper()
    info = "SVTYPE=INV;IMPRECISE;CIPOS={cip1},{cip2};CIEND={cie1},{cie2};" \
           "END={end};SVLEN={length}".format(cip1=start_pos - start_coords[0],
                                             cip2=start_coords[-1] - start_pos,
                                             cie1=end_pos - end_coords[0],
                                             cie2=end_coords[-1] - end_pos,
                                             end=end_pos,
                                             length=end_pos-start_pos)
    return VcfLine(x1.chrom1, start_pos, id1, base1, "<INV>", info)
1012037
def get_genome_refs(loc_file, loc_type):
    
    if not file_exists(loc_file):
        return None
    refs = {}
    with open(loc_file) as in_handle:
        for line in in_handle:
            if not line.startswith("#"):
                parts = line.strip().split()
                if loc_type in ["bowtie2", "samtools", "alignseq"]:
                    refs[parts[1]] = parts[-1]
                else:
                    refs[parts[0]] = parts[-1]
    return refs
247026
def _codon_positions(self, feature):
        
        results = self._transcript_feature_positions(feature)
        if len(results) != 3:
            raise ValueError(
                "Expected 3 positions for %s of %s but got %d" % (
                    feature,
                    self.id,
                    len(results)))
        return results
645005
def map_to_precursors(seqs, names, loci, out_file, args):
    """"""
    with make_temp_directory() as temp:
        pre_fasta = os.path.join(temp, "pre.fa")
        seqs_fasta = os.path.join(temp, "seqs.fa")
        out_sam = os.path.join(temp, "out.sam")
        pre_fasta = get_loci_fasta(loci, pre_fasta, args.ref)
        out_precursor_file = out_file.replace("tsv", "fa")
        seqs_fasta = get_seqs_fasta(seqs, names, seqs_fasta)

        # print(open(pre_fasta).read().split("\n")[1])
        if find_cmd("razers3"):
            cmd = "razers3 -dr 2 -i 80 -rr 90 -f -o {out_sam} {temp}/pre.fa  {seqs_fasta}"
            run(cmd.format(**locals()))
            out_file = read_alignment(out_sam, loci, seqs, out_file)
            shutil.copy(pre_fasta, out_precursor_file)
    return out_file
1174249
def six_frame(genome, table, minimum = 10):
    
    for seq in parse_fasta(genome):
        dna = Seq(seq[1].upper().replace('U', 'T'), IUPAC.ambiguous_dna)
        counter = 0
        for sequence in ['f', dna], ['rc', dna.reverse_complement()]:
            direction, sequence = sequence
            for frame in range(0, 3):
                for prot in \
                            sequence[frame:].\
                            translate(table = table, to_stop = False).split('*'):
                    if len(prot) < minimum:
                        continue
                    counter += 1
                    header = '%s_%s table=%s frame=%s-%s %s' % \
                                (seq[0].split()[0], counter, table, frame+1, \
                                direction, ' '.join(seq[0].split()[1:]))
                    yield [header, prot]
706195
def seq_info(names, id2names, insertions, sequences):
    
    seqs = {} # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns]], ...]]
    for name in names:
        id = id2names[name]
        gene = name.split('fromHMM::', 1)[0].rsplit(' ', 1)[1]
        model = name.split('fromHMM::', 1)[1].split('=', 1)[1].split()[0]
        i_gene_pos = insertions[id] # coordinates of each insertion wrt gene
        i_model_pos = name.split('fromHMM::', 1)[1].split('model-pos(ins-len)=')[1].split()[0].split(';') # model overlap
        i_info = []
        for i, ins in enumerate(i_gene_pos):
            model_pos = i_model_pos[i].split('-')[1].split('(')[0]
            length = i_model_pos[i].split('(')[1].split(')')[0]
            iheader = '>%s_%s insertion::seq=%s type=insertion strand=n/a gene-pos=%s-%s model-pos=%s'\
                    % (id, (i + 1), (i + 1), ins[0], ins[1], model_pos)
            iseq = sequences[id][1][ins[0]:(ins[1] + 1)]
            iseq = [iheader, iseq]
            info = [ins, model_pos, length, iseq, [], []]
            i_info.append(info)
        seqs[id] = [gene, model, i_info]
    return seqs
831514
private String _insdc_feature_position_string(Point location, int increment) {
			StringBuilder sb = new StringBuilder();
			Formatter formatter = new Formatter(sb,Locale.US);
			formatter.format("%s", location.getPosition() + increment);
			String output = formatter.toString();
			formatter.close();
			return output;

	/*
	if isinstance(pos, SeqFeature.ExactPosition):
		return "%i" % (pos.position+offset)
	elif isinstance(pos, SeqFeature.WithinPosition):
		return "(%i.%i)" % (pos.position + offset,
							pos.position + pos.extension + offset)
	elif isinstance(pos, SeqFeature.BetweenPosition):
		return "(%i^%i)" % (pos.position + offset,
							pos.position + pos.extension + offset)
	elif isinstance(pos, SeqFeature.BeforePosition):
		return "<%i" % (pos.position + offset)
	elif isinstance(pos, SeqFeature.AfterPosition):
		return ">%i" % (pos.position + offset)
	elif isinstance(pos, SeqFeature.OneOfPosition):
		return "one-of(%s)" \
			   % ",".join([_insdc_feature_position_string(p,offset) \
						   for p in pos.position_choices])
	elif isinstance(pos, SeqFeature.AbstractPosition):
		raise NotImplementedError("Please report this as a bug in Biopython.")
	else:
		raise ValueError("Expected a SeqFeature position object.")
		 */
	}
620417
def add_flanking_seqs(self, ref_seq, new_start, new_end):
        ''''''
        if new_start > self.POS or new_end < self.ref_end_pos():
            raise Error('new start and end positions must not try to shrink VCF record. new_start=' + str(new_start) + ', new_end=' + str(new_end) + '. VCF=' + str(self))

        new_start_nucleotides = ref_seq[new_start:self.POS]
        new_end_nucleotodes = ref_seq[self.ref_end_pos() + 1:new_end + 1]
        self.POS = new_start
        self.REF = new_start_nucleotides + self.REF + new_end_nucleotodes
        self.ALT = [new_start_nucleotides + x + new_end_nucleotodes for x in self.ALT]
1081527
def pair_hmm_align_unaligned_seqs(seqs, moltype=DNA_cogent, params={}):
    

    seqs = LoadSeqs(data=seqs, moltype=moltype, aligned=False)
    try:
        s1, s2 = seqs.values()
    except ValueError:
        raise ValueError(
            "Pairwise aligning of seqs requires exactly two seqs.")

    try:
        gap_open = params['gap_open']
    except KeyError:
        gap_open = 5
    try:
        gap_extend = params['gap_extend']
    except KeyError:
        gap_extend = 2
    try:
        score_matrix = params['score_matrix']
    except KeyError:
        score_matrix = make_dna_scoring_dict(
            match=1, transition=-1, transversion=-1)

    return local_pairwise(s1, s2, score_matrix, gap_open, gap_extend)
1126892
def positional_encoding(inputs,
                        num_units=None,
                        zero_pad=True,
                        scale=True,
                        scope="positional_encoding",
                        reuse=None):
    '''
    
    '''
    Shape = tf.shape(inputs)
    N = Shape[0]
    T = Shape[1]
    num_units = Shape[2]
    with tf.variable_scope(scope, reuse=reuse):
        position_ind = tf.tile(tf.expand_dims(tf.range(T), 0), [N, 1])

        # First part of the PE function: sin and cos argument
        #  Second part, apply the cosine to even columns and sin to odds.
        X = tf.expand_dims(tf.cast(tf.range(T), tf.float32), axis=1)
        Y = tf.expand_dims(
            tf.cast(10000 ** -(2 * tf.range(num_units) / num_units), tf.float32), axis=0)
        h1 = tf.cast((tf.range(num_units) + 1) % 2, tf.float32)
        h2 = tf.cast((tf.range(num_units) % 2), tf.float32)
        position_enc = tf.multiply(X, Y)
        position_enc = tf.sin(position_enc) * tf.multiply(tf.ones_like(X), h1) + \
            tf.cos(position_enc) * tf.multiply(tf.ones_like(X), h2)

        # Convert to a tensor
        lookup_table = position_enc

        if zero_pad:
            lookup_table = tf.concat((tf.zeros(shape=[1, num_units]),
                                      lookup_table[1:, :]), 0)
        outputs = tf.nn.embedding_lookup(lookup_table, position_ind)

        if scale:
            outputs = outputs * tf.sqrt(tf.cast(num_units, tf.float32))

        return outputs