obtain the sequence for a genomic region
356759
def get_genomic_seq_for_region(self, chrom, start_pos, end_pos):
        
        
        headers = {"content-type": "text/plain"}
        
        self.attempt = 0
        ext = "/sequence/region/human/{}:{}..{}:1".format(chrom, start_pos, end_pos)
        
        return self.ensembl_request(ext, headers)
1130037
def pp_event(seq):
    """"""

    if isinstance(seq, Event):
        return str(seq)

    # Get the original sequence back if seq is a pretty name already
    rev_curses = dict((v, k) for k, v in CURSES_NAMES.items())
    rev_curtsies = dict((v, k) for k, v in CURTSIES_NAMES.items())
    if seq in rev_curses:
        seq = rev_curses[seq]
    elif seq in rev_curtsies:
        seq = rev_curtsies[seq]

    pretty = curtsies_name(seq)
    if pretty != seq:
        return pretty
    return repr(seq).lstrip('u')[1:-1]
682110
def get_variants_in_region(self, chrom, start, end):
        """"""
        region = self.get_vcf()(
            "{}:{}-{}".format(chrom, start, end)
        )
        for v in region:
            for coded_allele, g in self._make_genotypes(v.ALT, v.genotypes):
                variant = Variant(
                    v.ID, v.CHROM, v.POS, [v.REF, coded_allele]
                )
                yield Genotypes(variant, g, v.REF, coded_allele,
                                multiallelic=len(v.ALT) > 1)
218926
def get_gene_seqs(database_path, gene):
    
    gene_path = database_path + "/" + gene + ".fsa"
    gene_seq = ""
    # Open fasta file
    with open(gene_path) as gene_file:
        header = gene_file.readline()
        for line in gene_file:
            seq = line.strip()
            gene_seq += seq
    return gene_seq
296321
def get_sequences(self, chr, coords):
            
        # Check if we have an index_dir
        if not self.index_dir:
            print("Index dir is not defined!")
            sys.exit()

        # retrieve all information for this specific sequence
        fasta_file = self.fasta_file[chr]
        index_file = self.index_file[chr]
        line_size = self.line_size[chr]
        total_size = self.size[chr]
        index = open(index_file, "rb")
        fasta = open(fasta_file)
        
        seqs = []
        for coordset in coords:
            seq = ""
            for (start,end) in coordset: 
                if start > total_size:
                    raise ValueError("%s: %s, invalid start, greater than sequence length!" % (chr,start))
            
                if start < 0:
                    raise ValueError("Invalid start, < 0!")
                
                if end > total_size:
                    raise ValueError("Invalid end, greater than sequence length!")


                seq += self._read(index, fasta, start, end, line_size)
            seqs.append(seq)
        index.close()
        fasta.close()

        return seqs
1264765
def sequence(self):
        
        seq = [x.mol_letter for x in self._monomers]
        return ''.join(seq)
89495
def get_variants_in_region(self, chrom, start, end):
        """"""
        bim = self.bim.loc[
            (self.bim["chrom"] == CHROM_STR_TO_INT[chrom]) &
            (start <= self.bim["pos"]) &
            (self.bim["pos"] <= end)
        ]
        for i, g in enumerate(self.bed.iter_geno_marker(bim.index)):
            info = bim.iloc[i, :]
            name, geno = g
            yield Genotypes(
                Variant(info.name, CHROM_INT_TO_STR[info.chrom],
                        info.pos, [info.a1, info.a2]),
                self._normalize_missing(geno),
                reference=info.a2,
                coded=info.a1,
                multiallelic=info.multiallelic
            )
693187
def find_seq_homologues(self, return_raw=False):
        
        
        
        # Get sequence from structure/chain
        # We could use Bio.PDB.Polypeptide?
        
        from Bio.SCOP.Raf import to_one_letter_code
        
        s = self
        seq_iter = s.get_residues()
        seq_str = ''

        for aa in seq_iter:
            if aa.resname in to_one_letter_code:
                seq_str += to_one_letter_code[aa.resname]
        
        # Use BLAST to find homologous sequences with associated
        # structures in the PDB database.
        # Perhaps include local BLAST?
        
        from Bio.Blast.NCBIWWW import qblast

        # Adapt for short query sequences if needed
        # From http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#shortp
        
        if len(seq_str) < 15:
            word_size = 2
            expect = 20000
            matrix_name = 'PAM30'
            filter = None
        else:
            word_size = 3
            expect = 10.0
            matrix_name = 'BLOSUM62'
            filter = 'SEG'
             
        query_result = qblast(  "blastp", "pdb", seq_str, 
                                word_size, expect, matrix_name, filter)
        
        if return_raw:
           return query_result 
            
        # Parse BLAST result to yield results
        # PDBID : (E-Value, Identity, Positives, Gaps, Alignment)
        
        from Bio.Blast import NCBIXML
        
        blast_records = NCBIXML.read(query_result)
        
        results = []
        
        for alignment in blast_records.alignments:
            for hsps in alignment.hsps:
                id_perc = "%s/%s" %(hsps.identities, alignment.length)
                pos_perc = "%s/%s" %(hsps.positives, alignment.length)
                gaps_perc = "%s/%s" %(hsps.gaps, alignment.length)

                pdb_id = alignment.title.split('|')[3]
                e_value = hsps.expect
                
                results.append( ( pdb_id,
                                "%2.5e" %e_value,
                                id_perc, 
                                pos_perc, 
                                gaps_perc,
                                '\n'.join([hsps.query, hsps.match, hsps.sbjct]) ) )

        return results
1337233
def _get_fastq_in_region(region, align_bam, out_base):
    
    out_file = "{0}.fastq".format(out_base)
    if not file_exists(out_file):
        with pysam.Samfile(align_bam, "rb") as in_pysam:
            with file_transaction(out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    contig, start, end = region
                    for read in in_pysam.fetch(contig, int(start), int(end)):
                        seq = Seq.Seq(read.seq)
                        qual = list(read.qual)
                        if read.is_reverse:
                            seq = seq.reverse_complement()
                            qual.reverse()
                        out_handle.write("@{name}\n{seq}\n+\n{qual}\n".format(
                                name=read.qname, seq=str(seq), qual="".join(qual)))
    return out_file
1005073
def _join_sequence(seq, last_separator=''):
    """"""
    count = len(seq)
    return ', '.join(_format_element(element, count, i, last_separator)
                     for i, element in enumerate(seq))
175037
def make_sequence(seq, name="S"):
    
    return ["{}_{}_{}".format(name, i, x) for i, x in enumerate(seq)]
61602
def sequence_dir(self):
        """"""
        if self.root_dir:
            return op.join(self.protein_dir, 'sequences')
        else:
            log.debug('Root directory not set')
            return None
124144
def seq_str(self):
        """"""

        if not self.seq:
            return None

        return ssbio.protein.sequence.utils.cast_to_str(self.seq)
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
276526
def _sequence(self, i):
        """"""

        result = ['[']
        end_range = 0
        escape_hyphen = -1
        removed = False
        last_posix = False

        c = next(i)
        if c in ('!', '^'):
            # Handle negate char
            result.append('^')
            c = next(i)
        if c == '[':
            last_posix = self._handle_posix(i, result, 0)
            if not last_posix:
                result.append(re.escape(c))
            c = next(i)
        elif c in ('-', ']'):
            result.append(re.escape(c))
            c = next(i)

        while c != ']':
            if c == '-':
                if last_posix:
                    result.append('\\' + c)
                    last_posix = False
                elif i.index - 1 > escape_hyphen:
                    # Found a range delimiter.
                    # Mark the next two characters as needing to be escaped if hyphens.
                    # The next character would be the end char range (s-e),
                    # and the one after that would be the potential start char range
                    # of a new range (s-es-e), so neither can be legitimate range delimiters.
                    result.append(c)
                    escape_hyphen = i.index + 1
                    end_range = i.index
                elif end_range and i.index - 1 >= end_range:
                    if self._sequence_range_check(result, '\\' + c):
                        removed = True
                    end_range = 0
                else:
                    result.append('\\' + c)
                c = next(i)
                continue
            last_posix = False

            if c == '[':
                last_posix = self._handle_posix(i, result, end_range)
                if last_posix:
                    c = next(i)
                    continue

            if c == '\\':
                # Handle escapes
                subindex = i.index
                try:
                    value = self._references(i, True)
                except PathNameException:
                    raise StopIteration
                except StopIteration:
                    i.rewind(i.index - subindex)
                    value = r'\\'
            elif c == '/':
                if self.pathname:
                    raise StopIteration
                value = c
            elif c in SET_OPERATORS:
                # Escape &, |, and ~ to avoid &&, ||, and ~~
                value = '\\' + c
            else:
                # Anything else
                value = c

            if end_range and i.index - 1 >= end_range:
                if self._sequence_range_check(result, value):
                    removed = True
                end_range = 0
            else:
                result.append(value)

            c = next(i)

        result.append(']')
        # Bad range removed.
        if removed:
            value = "".join(result)
            if value == '[]':
                # We specified some ranges, but they are all
                # out of reach.  Create an impossible sequence to match.
                result = ['[^%s]' % ('\x00-\xff' if self.is_bytes else uniprops.UNICODE_RANGE)]
            elif value == '[^]':
                # We specified some range, but hey are all
                # out of reach. Since this is exclusive
                # that means we can match *anything*.
                result = ['[%s]' % ('\x00-\xff' if self.is_bytes else uniprops.UNICODE_RANGE)]
            else:
                result = [value]

        if self.pathname or self.after_start:
            return self._restrict_sequence() + ''.join(result)

        return ''.join(result)
688879
def _extract_seq(extractor, seq):
    """"""
    subtype = Any
    if seq.__args__ and seq.__args__[0] is not Any:
        subtype = seq.__args__[0]
    return _array_type(extractor.extract(subtype))
1140752
def get_orfs(_orf_length = 100, prediction = @prediction)
      '-' if prediction.type != 'nucleotide'

      seq = Bio::Sequence::NA.new(prediction.raw_sequence)

      result = {}
      key = 0
      (1..6).each do |f|
        s = seq.translate(f)
        f = -1 if f == 4
        f = -2 if f == 5
        f = -3 if f == 6
        s.scan(/(\w{30,})/) do |_orf|
          orf_start = $LAST_MATCH_INFO.offset(0)[0] + 1
          orf_end   = $LAST_MATCH_INFO.offset(0)[1] + 1
          coverage = (((orf_end - orf_start) / s.length.to_f) * 100).ceil
          # reduce the orf_end and the translated length by 2% to increase the
          #   width between ORFs on the plot
          chopping = s.length * 0.02
          orf_end = (orf_end.to_f - chopping).floor
          translated_length = (s.length - chopping).ceil
          key += 1
          result[key] = { frame: f, orf_start: orf_start, orf_end: orf_end,
                          coverage: coverage,
                          translated_length: translated_length }
        end
      end
      result
    end