Create exon object for all exons in the database, optionally
        restrict to a particular chromosome using the `contig` argument.
492865
def exons(self, contig=None, strand=None):
        
        # DataFrame with single column called "exon_id"
        exon_ids = self.exon_ids(contig=contig, strand=strand)
        return [
            self.exon_by_id(exon_id)
            for exon_id in exon_ids
        ]
164741
def start_range(self):
      """"""
      if len(self._exons) == 0: return None
      return GenomicRange(self._exons[0].chr,
             min([x.start for x in self._exons]),# must be part of junction
             max([x.start for x in self._exons]))
1231139
def add_genes(in_file, data, max_distance=10000, work_dir=None):
    
    gene_file = regions.get_sv_bed(data, "exons", out_dir=os.path.dirname(in_file))
    if gene_file and utils.file_exists(in_file):
        out_file = "%s-annotated.bed" % utils.splitext_plus(in_file)[0]
        if work_dir:
            out_file = os.path.join(work_dir, os.path.basename(out_file))
        if not utils.file_uptodate(out_file, in_file):
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            with file_transaction(data, out_file) as tx_out_file:
                _add_genes_to_bed(in_file, gene_file, fai_file, tx_out_file, data, max_distance)
        return out_file
    else:
        return in_file
1050511
def chr(self):
    
    if len(self.exons)==0: 
      sys.stderr.write("WARNING can't return chromsome with nothing here\n")
      return None
    return self._rngs[0].chr
460551
def GetPeaksExons(bed,parsedGTF):
    

    bedtool_AB=dfTObedtool(bed)

    exonsGTF=parsedGTF[parsedGTF["feature"]=="exon"]
    exonsGTF.reset_index(inplace=True, drop=True)

    exonsBED=GTFtoBED(exonsGTF, "exon_id")
    exonsBED.columns=['chrom', 'chromStart', 'chromEnd', 'exon_id', 'score', 'strand']
    exonsBEDcols=exonsBED.columns.tolist()

    bedcols=bed.columns.tolist()
    exonsBEDcols_=[]
    for c in exonsBEDcols:
        if c in bedcols:
            exonsBEDcols_.append(c+"_exon")
        else:
            exonsBEDcols_.append(c)

    cols=[bedcols,exonsBEDcols_,["overlap"] ]
    cols=[item for sublist in cols for item in sublist]

    bedtool_exons=dfTObedtool(exonsBED)

    bedtool_target_exons=bedtool_AB.intersect(bedtool_exons, wo=True, s=True)
    dfTargetE=pd.read_table(bedtool_target_exons.fn, names=cols)
    ExonsTransGenes=parsedGTF[["exon_id","transcript_id","gene_id"]].drop_duplicates()
    dfTargets=pd.merge(dfTargetE,ExonsTransGenes,on=["exon_id"],how="left")
    dfTargets["count"]=1

    def getCounts(df,field):
        """
        For each field in a bed narrow peak returns the number or times that field is present,\
        the normalized mean of the '-log10(pValue)' and normalized mean of the signal value.

        :param df: a Pandas dataframe of a bed narrow peak
        :param field: field to analyse, ie. exons or transcripts

        :returns: a Pandas dataframe
        """

        tmp=df[[field,'name',"count"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["count"]=tmp["count"].astype(int)
        tmp.columns=[field,"%s_count" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'name',"-log10(pValue)"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["-log10(pValue)"]=tmp["-log10(pValue)"].astype(float)
        tmp=tmp.groupby(field).apply(lambda l: reduce(lambda x, y: x*y, l["-log10(pValue)"]) )
        tmp=pd.DataFrame(tmp)
        tmp.reset_index(inplace=True,drop=False)
        tmp.columns=[field,"%s norm. mean -log10(pValue)" %str(field)]
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'name',"signalValue"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["signalValue"]=tmp["signalValue"].astype(float)
        tmp=tmp.groupby(field).apply(lambda l: reduce(lambda x, y: x*y, l["signalValue"]) )
        tmp=pd.DataFrame(tmp)
        tmp.reset_index(inplace=True,drop=False)
        tmp.columns=[field,"%s signalValue" %str(field)]
        df=pd.merge(df,tmp,on=field,how="left")

        return df

    for f in ["exon_id","transcript_id"]:
        dfTargets=getCounts(dfTargets,f)

    def getCounts_GeneIDs(df):
        """
        For each gene id in a bed narrow peak returns the number or times that field is present,\
        the normalized mean of the '-log10(pValue)' and normalized mean of the signal value.

        :param df: a Pandas dataframe of a bed narrow peak

        :returns: a Pandas dataframe
        """

        field="gene_id"

        tmp=df[[field,"transcript_id","transcript_id_count"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id_count"]=tmp["transcript_id_count"].astype(int)
        tmp.columns=[field,"%s_count" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'transcript_id',"transcript_id norm. mean -log10(pValue)"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id norm. mean -log10(pValue)"]=tmp["transcript_id norm. mean -log10(pValue)"].astype(float)
        tmp.columns=[field,"%s norm. mean -log10(pValue)" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")



        tmp=df[[field,'transcript_id',"transcript_id signalValue"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id signalValue"]=tmp["transcript_id signalValue"].astype(float)
        tmp.columns=[field,"%s signalValue" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        return df

    dfTargets=getCounts_GeneIDs(dfTargets)


    dfTargets=dfTargets.drop(["count"],axis=1)
    return dfTargets
587628
def encodePolymorphicNucleotide(polySeq) :
	"""'"""
	
	if type(polySeq) is types.StringType :
		if polySeq.find("/") < 0 :
			sseq = list(polySeq)
		else :
			sseq = polySeq.split('/')
			
	else :
		sseq = polySeq
	
	seq = []
	for n in sseq :
		try :
			for n2 in polymorphicNucleotides[n] :
				seq.append(n2)
		except KeyError :
			seq.append(n)
	
	seq = set(seq)
	
	if len(seq) == 4:
		return 'N'
	elif len(seq) == 3 :
		if 'T' not in seq :
			return 'V'
		elif 'G' not in seq :
			return 'H'
		elif 'C' not in seq :
			return 'D'
		elif 'A' not in seq :
			return 'B'
	elif len(seq) == 2 :
		if 'A' in seq and 'G' in seq :
			return 'R'
		elif 'C' in seq and 'T' in seq :
			return 'Y'
		elif 'A' in seq and 'C' in seq :
			return 'M'
		elif 'T' in seq and 'G' in seq :
			return 'K'
		elif 'A' in seq and 'T' in seq :
			return 'W'
		elif 'C' in seq and 'G' in seq :
			return 'S'
	elif polySeq[0] in nucleotides :
		return polySeq[0]
	else :
		raise UnknownNucleotide(polySeq)
655228
def runblast(self, assembly, allele, sample):
        
        genome = os.path.split(assembly)[1].split('.')[0]
        # Run the BioPython BLASTn module with the genome as query, fasta(target gene) as db.
        # Do not re-perform the BLAST search each time
        make_path(sample[self.analysistype].reportdir)
        try:
            report = glob('{}{}*rawresults*'.format(sample[self.analysistype].reportdir, genome))[0]
            size = os.path.getsize(report)
            if size == 0:
                os.remove(report)
                report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                          time.strftime("%Y.%m.%d.%H.%M.%S"))
        except IndexError:
            report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                      time.strftime("%Y.%m.%d.%H.%M.%S"))
        db = allele.split('.')[0]
        # BLAST command line call. Note the mildly restrictive evalue, and the high number of alignments.
        # Due to the fact that all the targets are combined into one database, this is to ensure that all potential
        # alignments are reported. Also note the custom outfmt: the doubled quotes are necessary to get it work
        blastn = NcbiblastnCommandline(query=assembly, db=db, evalue='1E-20', num_alignments=1000000,
                                       num_threads=12,
                                       outfmt="'6 qseqid sseqid positive mismatch gaps "
                                              "evalue bitscore slen length qstart qend qseq sstart send'",
                                       out=report)
        # Save the blast command in the metadata
        sample[self.analysistype].blastcommand = str(blastn)
        sample[self.analysistype].blastreport = report
        if not os.path.isfile(report):
            # Run BLAST
            blastn()
        # Run the blast parsing module
        self.blastparser(report, sample)
958671
def _init_exons(self):
        
        exon_starts = [self.chrom_start + int(s)
                       for s in self.bed_tuple.blockStarts.strip(',').split(',')]
        exon_sizes = list(map(int, self.bed_tuple.blockSizes.strip(',').split(',')))

        # get chromosome intervals
        exons = [(exon_starts[i], exon_starts[i] + exon_sizes[i])
                 for i in range(len(exon_starts))]
        no_utr_exons = self._filter_utr(exons)
        self.exons = no_utr_exons
        self.exon_lens = [e[1] - e[0] for e in self.exons]
        self.num_exons = len(self.exons)
        self.cds_len = sum(self.exon_lens)
        self.five_ss_len = 2*(self.num_exons-1)
        self.three_ss_len = 2*(self.num_exons-1)
        self._init_splice_site_pos()
317093
private boolean isAlleleCorrect(String allele, boolean acceptAmbiguousBases) {
        if (StringUtils.isNotEmpty(allele)) {
            for (char base : allele.toCharArray()) {
                if (!isValidBase(base, acceptAmbiguousBases)) {
                    return false;
                }
            }
        }
        return true;
    }
136325
def _splitGenoGeneWindow(self,annotation_file=None,cis=1e4,funct='protein_coding',minSnps=1.,maxSnps=SP.inf):
        
         #1. load annotation
        assert annotation_file is not None, 'Splitter:: specify annotation file'
        try:
            f = h5py.File(annotation_file,'r')
            geneID = f['geneID'][:]
            gene_chrom = f['chrom'][:]
            gene_start = f['start'][:]
            gene_end = f['end'][:]
            gene_strand = f['strand'][:]
            gene_function = f['function'][:]
            f.close()
        except:
            print('Splitter:: format annotation file not valid')

        # if funct is not None, it has to be a list
        if funct is not None and funct!=list:   funct=[funct]

        windows = []
        nSnps   = []
        Igene   = []
        #2. calculates windows 
        for gene_i in range(geneID.shape[0]):
            if funct is not None:
                if gene_function[gene_i] not in funct:
                    Igene.append(False)
                    continue
            wnd = [gene_chrom[gene_i],gene_start[gene_i]-cis,gene_end[gene_i]+cis]
            Ir = (self.chrom==wnd[0])*(self.pos>=wnd[1])*(self.pos<=wnd[2])
            _nSnps = Ir.sum()
            if _nSnps>=minSnps and _nSnps<=maxSnps:
                windows.append(wnd)
                nSnps.append(_nSnps)
                Igene.append(True)
            else:
                Igene.append(False)
        Igene = SP.array(Igene)
        self.info['nSnps'] = SP.array(nSnps)
        self.info['geneID']        = geneID[Igene]
        self.info['gene_start']    = gene_start[Igene]
        self.info['gene_end']      = gene_end[Igene]
        self.info['gene_strand']   = gene_strand[Igene]
        self.info['gene_function'] = gene_function[Igene]
        return SP.array(windows)
480751
def make_gffutils_db(gtf, db):
    
    import gffutils
    out_db = gffutils.create_db(gtf,
                                db,
                                keep_order=True,
                                infer_gene_extent=False)
    return out_db
563233
def codon2aa(codon, trans_table):
    
    return Seq(''.join(codon), IUPAC.ambiguous_dna).translate(table = trans_table)[0]
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)
1046092
def _incorporate_delins(self):
        """"""
        seq, cds_start, cds_stop, start, end = self._setup_incorporate()

        ref = self._var_c.posedit.edit.ref
        alt = self._var_c.posedit.edit.alt
        ref_length = end - start if ref is not None else 0    # can't just get from ref since ref isn't always known
        alt_length = len(
            self._var_c.posedit.edit.alt) if self._var_c.posedit.edit.alt is not None else 0
        net_base_change = alt_length - ref_length
        cds_stop += net_base_change

        # incorporate the variant into the sequence (depending on the type)
        is_substitution = False
        if ref is not None and alt is not None:    # delins or SNP
            seq[start:end] = list(alt)
            if len(ref) == 1 and len(alt) == 1:
                is_substitution = True
        elif ref is not None:    # deletion
            del seq[start:end]
        else:    # insertion
            seq[start + 1:start + 1] = list(alt)    # insertion in list before python list index

        if DBG:
            print("net base change: {}".format(net_base_change))
        is_frameshift = net_base_change % 3 != 0
        # use max of mod 3 value and 1 (in event that indel starts in the 5'utr range)
        variant_start_aa = max(int(math.ceil((self._var_c.posedit.pos.start.base) / 3.0)), 1)

        alt_data = AltTranscriptData(
            seq,
            cds_start,
            cds_stop,
            is_frameshift,
            variant_start_aa,
            self._transcript_data.protein_accession,
            is_substitution=is_substitution,
            is_ambiguous=self._ref_has_multiple_stops)
        return alt_data
1152699
def copy(self):
        
        return type(self)(self.dna, suppressed=self.suppressed, name=self.name)
780980
def iTable2GFF(iTable, fa, contig = False):
    
    columns = ['#seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame', 'attribute']
    gff = {c:[] for c in columns}
    for insertion in iTable.iterrows():
        insertion = insertion[1]
        if insertion['ID'] not in fa:
            continue
        # rRNA strand
        strand = insertion['sequence'].split('strand=', 1)[1].split()[0]
        # set rRNA positions for reporting features on contig or extracted sequence
        if contig is True:
            gene = [int(i) for i in insertion['sequence'].split('pos=', 1)[1].split()[0].split('-')]
            if strand == '-':
                offset = -1 * (gene[1])
            else:
                offset = gene[0]
        else:
            strand = '+'
            gene = [1, int(insertion['sequence'].split('total-len=', 1)[1].split()[0])]
            offset = gene[0]
        insertion['strand'] = strand
        insertion['offset'] = offset
        # source for prediction
        source = insertion['sequence'].split('::model', 1)[0].rsplit(' ', 1)[-1]
        insertion['source'] = source
        # rRNA gene
        geneAnnot = '%s rRNA gene' % (source.split('from', 1)[0])
        geneNum = insertion['sequence'].split('seq=', 1)[1].split()[0]
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(source)
        gff['feature'].append('Gene')
        gff['start'].append(gene[0])
        gff['end'].append(gene[1])
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (geneNum, geneAnnot))
        # rRNA
        gff = parse_rRNA(insertion, fa[insertion['ID']], gff)
        # insertions
        gff = parse_insertion(insertion, gff)
        # orfs
        gff = parse_orf(insertion, gff)
        # catalytic RNAs
        gff = parse_catalytic(insertion, gff)
    return pd.DataFrame(gff)[columns].drop_duplicates()
1134936
def set_table(genome, table, table_name, connection_string, metadata):
    
    table = Table(table_name, genome._metadata, autoload=True,
                    autoload_with=genome.bind, extend_existing=True)

    #print "\t".join([c.name for c in table.columns])
    # need to prefix the indexes with the table name to avoid collisions
    for i, idx in enumerate(table.indexes):
        idx.name = table_name + "." + idx.name + "_ix" + str(i)

    cols = []
    for i, col in enumerate(table.columns):
        # convert mysql-specific types to varchar
        #print col.name, col.type, isinstance(col.type, ENUM)
        if isinstance(col.type, (LONGBLOB, ENUM)):

            if 'sqlite' in connection_string:
                col.type = VARCHAR()
            elif 'postgres' in connection_string:
                if isinstance(col.type, ENUM):
                    #print dir(col)
                    col.type = PG_ENUM(*col.type.enums, name=col.name,
                        create_type=True)
                else:
                    col.type = VARCHAR()
        elif str(col.type) == "VARCHAR" \
                and ("mysql" in connection_string \
                or "postgres" in connection_string):
            if col.type.length is None:
                col.type.length = 48 if col.name != "description" else None
        if not "mysql" in connection_string:
            if str(col.type).lower().startswith("set("):
                col.type = VARCHAR(15)
        cols.append(col)

    table = Table(table_name, genome._metadata, *cols,
            autoload_replace=True, extend_existing=True)

    return table