Given a FASTA database, proteins are trypsinized and resulting peptides
    stored in a database or dict for lookups
189322
def create_searchspace(lookup, fastafn, proline_cut=False,
                       reverse_seqs=True, do_trypsinize=True):
    """"""
    allpeps = []
    for record in SeqIO.parse(fastafn, 'fasta'):
        if do_trypsinize:
            pepseqs = trypsinize(record.seq, proline_cut)
        else:
            pepseqs = [record.seq]
        # Exchange all leucines to isoleucines because MS can't differ
        pepseqs = [(str(pep).replace('L', 'I'),) for pep in pepseqs]
        allpeps.extend(pepseqs)
        if len(allpeps) > 1000000:  # more than x peps, write to SQLite
            lookup.write_peps(allpeps, reverse_seqs)
            allpeps = []
    # write remaining peps to sqlite
    lookup.write_peps(allpeps, reverse_seqs)
    lookup.index_peps(reverse_seqs)
    lookup.close_connection()
1320219
def fetch(self, is_dl_forced=True):
        

        username = config.get_config()['dbauth']['udp']['user']
        password = config.get_config()['dbauth']['udp']['password']

        credentials = (username, password)

        # Get patient map file:
        patient_id_map = self.open_and_parse_yaml(self.map_files['patient_ids'])
        udp_internal_ids = patient_id_map.keys()

        phenotype_fields = ['Patient', 'HPID', 'Present']

        # Get phenotype ids for each patient
        phenotype_params = {
            'method': 'search_subjects',
            'subject_type': 'Phenotype',
            'search_mode': 'DEEP',
            'fields': 'Patient',
            'conditions': 'equals',
            'values': ','.join(udp_internal_ids),
            'user_fields': ','.join(phenotype_fields)
        }

        prioritized_variants = [
            'Patient', 'Gene', 'Chromosome Position', 'Variant Allele', 'Transcript']

        prioritized_params = {
            'method': 'search_subjects',
            'subject_type': 'Variant Prioritization',
            'search_mode': 'DEEP',
            'fields': 'Patient',
            'conditions': 'equals',
            'values': ','.join(udp_internal_ids),
            'user_fields': ','.join(prioritized_variants),
            'format': 'json'}

        variant_fields = [
            'Patient', 'Family', 'Chr', 'Build', 'Chromosome Position',
            'Reference Allele', 'Variant Allele', 'Parent of origin',
            'Allele Type', 'Mutation Type', 'Gene', 'Transcript', 'Original Amino Acid',
            'Variant Amino Acid', 'Amino Acid Change', 'Segregates with',
            'Position', 'Exon', 'Inheritance model', 'Zygosity', 'dbSNP ID',
            '1K Frequency', 'Number of Alleles']

        variant_params = {
            'method': 'search_subjects',
            'subject_type': 'Exome Analysis Results',
            'search_mode': 'DEEP',
            'fields': 'Patient',
            'conditions': 'equals',
            'user_fields': ','.join(variant_fields),
            'format': 'json'}

        pheno_file = open(
            '/'.join((self.rawdir, self.files['patient_phenotypes']['file'])), 'w')
        variant_file = open(
            '/'.join((self.rawdir, self.files['patient_variants']['file'])), 'w')

        pheno_file.write('{0}\n'.format('\t'.join(phenotype_fields)))
        variant_file.write('{0}\n'.format('\t'.join(variant_fields)))

        variant_gene = self._fetch_data_from_udp(
            udp_internal_ids, prioritized_params, prioritized_variants, credentials)

        variant_gene_map = dict()
        for line in variant_gene:
            variant_gene_map.setdefault(line[0], []).append(
                # Try to make a unique value based on gene-pos-variantAlele-transcript
                # TODO make this a dict for readability purposes
                "{0}-{1}-{2}-{3}".format(line[1], line[2], line[3], line[4]))

        variant_info = self._fetch_data_from_udp(
            udp_internal_ids, variant_params, variant_fields, credentials)

        for line in variant_info:
            variant = "{0}-{1}-{2}-{3}".format(line[10], line[4], line[6], line[11])
            if variant in variant_gene_map[line[0]]:
                line[0] = patient_id_map[line[0]]
                line[4] = re.sub(r'\.0$', '', line[4])
                variant_file.write('{0}\n'.format('\t'.join(line)))

        phenotype_info = self._fetch_data_from_udp(
            udp_internal_ids, phenotype_params, phenotype_fields, credentials)

        for line in phenotype_info:
            line[0] = patient_id_map[line[0]]
            pheno_file.write('{0}\n'.format('\t'.join(line)))

        variant_file.close()
        pheno_file.close()

        return
1008397
def generate_combined_fasta(self, genome_list, genome_dir):
        '''
        '''

        fasta = []
            
        for genome in genome_list:
            
            full_path = genome_dir + genome
            handle = open(full_path, "rU")
            print 'making combined fasta for', genome
            try:
                seq_record = SeqIO.read(handle, 'genbank')
                org_accession = seq_record.name
                
            except AssertionError,e:
                print str(e), genome
                        
            for feature in seq_record.features:
                if feature.type == 'CDS':
                    try:                                  
                        prot_accession = feature.qualifiers['protein_id'][0]
                        prot_translation = feature.qualifiers['translation'][0]
                        
                        newfast = '>' + org_accession + ',' + prot_accession + \
                         '\n' + prot_translation + '\n'
                        
                        #if newfast not in fasta:
                        fasta.append(newfast)
                            
                    except AttributeError,e:
                        print "organism %s, protein %s did not have \
                        the right attributes" % (org_accession, prot_accession)
                        print str(e)
                    
                    except KeyError,e:
                        print "organism %s, protein %s did not have \
                        the right key" % (org_accession, prot_accession)
                        print str(e)
                          
            handle.close()
            
            print "%s proteins were added" % len(fasta)

            set_fasta = set(fasta)

            print "%s unique proteins were added -- dropping redundant ones" % len(set_fasta)
                                                
        faastring = "".join(set_fasta)

        write_fasta = open('combined_fasta', 'w')
        write_fasta.write(faastring)
        write_fasta.close()

        return set_fasta
128310
def sort_protein_sequences(protein_sequences):
    
    return list(
        sorted(
            protein_sequences,
            key=ProteinSequence.ascending_sort_key,
            reverse=True))
313893
def extract_db_info(self, obj, keys):
        """"""
        result = {}
        if isinstance(obj, dict):
            try:
                qc = obj['quality_control']
            except KeyError:
                qc = QC.UNKNOWN
        elif isinstance(obj, DataFrame):
            with obj.open() as hdulist:
                qc = self.datamodel.get_quality_control(hdulist)
        else:
            qc = QC.UNKNOWN

        result['quality_control'] = qc
        other = super(DataProductMixin, self).extract_db_info(obj, keys)
        result.update(other)
        return result
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
1109120
def protein_sequences_generator_to_dataframe(variant_and_protein_sequences_generator):
    
    return dataframe_from_generator(
        element_class=ProteinSequence,
        variant_and_elements_generator=variant_and_protein_sequences_generator,
        converters=dict(
            gene=lambda x: ";".join(x)))
759328
def sb_filter(fastq, bc, cores, nedit):
    ''' 
    '''
    barcodes = set(sb.strip() for sb in bc)
    if nedit == 0:
        filter_sb = partial(exact_sample_filter2, barcodes=barcodes)
    else:
        barcodehash = MutationHash(barcodes, nedit)
        filter_sb = partial(correcting_sample_filter2, barcodehash=barcodehash)
    p = multiprocessing.Pool(cores)

    chunks = tz.partition_all(10000, read_fastq(fastq))
    bigchunks = tz.partition_all(cores, chunks)
    for bigchunk in bigchunks:
        for chunk in p.map(filter_sb, list(bigchunk)):
            for read in chunk:
                sys.stdout.write(read)
396978
def protein_sequence(self, protein_id):
        
        if self.protein_sequences is None:
            raise ValueError(
                "No protein FASTA supplied to this Genome: %s" % self)
        return self.protein_sequences.get(protein_id)
1040430
def demultiplex_samples(fastq, out_dir, nedit, barcodes):
    ''' 
    '''
    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    if barcodes:
        barcodes = set(barcode.strip() for barcode in barcodes)
    else:
        barcodes = set()

    if nedit == 0:
        filter_bc = partial(exact_sample_filter, barcodes=barcodes)
    else:
        barcodehash = MutationHash(barcodes, nedit)
        filter_bc = partial(correcting_sample_filter, barcodehash=barcodehash)

    sample_set = set()
    batch = collections.defaultdict(list)
    parsed = 0
    safe_makedir(out_dir)
    for read in read_fastq(fastq):
        parsed += 1
        read = filter_bc(read)
        if not read:
            continue
        match = parser_re.search(read).groupdict()
        sample = match['SB']
        sample_set.add(sample)
        batch[sample].append(read)
        # write in batches to avoid opening up file handles repeatedly
        if not parsed % 10000000:
            for sample, reads in batch.items():
                out_file = os.path.join(out_dir, sample + ".fq")
                with open(out_file, "a") as out_handle:
                    for read in reads:
                        fixed = filter_bc(read)
                        if fixed:
                            out_handle.write(fixed)
            batch = collections.defaultdict(list)

    for sample, reads in batch.items():
        out_file = os.path.join(out_dir, sample + ".fq")
        with open(out_file, "a") as out_handle:
            for read in reads:
                fixed = filter_bc(read)
                if fixed:
                    out_handle.write(read)
523769
def create_merge_psm_map(peptides, ns):
    """"""
    psmmap = {}
    for peptide in peptides:
        seq = reader.get_peptide_seq(peptide, ns)
        psm_ids = reader.get_psm_ids_from_peptide(peptide, ns)
        for psm_id in psm_ids:
            try:
                psmmap[seq][psm_id.text] = 1
            except KeyError:
                psmmap[seq] = {psm_id.text: 2}
    for seq, psm_id_dict in psmmap.items():
        psmmap[seq] = [x for x in psm_id_dict]
    return psmmap
1101348
def reference_fasta(self):
        """"""
        if self._db_location:
            ref_files = glob.glob(os.path.join(self._db_location, "*", self._REF_FASTA))
            if ref_files:
                return ref_files[0]
712932
def create_isobaric_quant_lookup(fns, tablefn_map, featmap, pqdb,
                                 featcolnr, allquantcols, psmcolmap):
    """"""
    pqdb.store_quant_channels(
            map_psmnrcol_to_quantcol(allquantcols, psmcolmap, tablefn_map),
            psmcolmap)
    quantmap = pqdb.get_quantchannel_map()
    to_store = []
    for fn, header, pquant in tsvreader.generate_tsv_pep_protein_quants(fns):
        pqdata = get_isob_quant_data(pquant, header, tablefn_map[fn],
                                     featmap, featcolnr, quantmap)
        to_store.extend(pqdata)
        if len(to_store) > 10000:
            pqdb.store_isobaric_quants(to_store, psmcolmap)
            to_store = []
    pqdb.store_isobaric_quants(to_store, psmcolmap)
291260
def _guess_seqtype_from_file(handle):
    ""
    if isinstance(handle, basestring):
        handle = StringIO(handle)

    for line in handle:
        if not line.strip():
            continue
        if line.lstrip().split()[0] in ('LOCUS', 'FEATURES', 'source', 'CDS',
                                        'gene'):
            return 'genbank'
        if len(line) > 2 and line[:3] in ('ID ', 'FT '):
            return 'embl'
        if line.startswith('>'):
            return 'fasta'
    handle.seek(0)
    import string
    from Bio.Data import IUPACData as iupac
    all_input_letters = set(handle.read().lower())
    all_valid = set(string.digits)
    all_valid.update(set(iupac.protein_letters.lower()))
    all_valid.update(set(iupac.unambiguous_dna_letters.lower()))
    all_valid.update(set('- \n'))
    if all_valid.issuperset(all_input_letters):
        return 'fasta'

    raise ValueError("Failed to guess format for input")
983376
def entrez(args):
    
    p = OptionParser(entrez.__doc__)

    allowed_databases = {"fasta": ["genome", "nuccore", "nucgss", "protein", "nucest"],
                         "asn.1": ["genome", "nuccore", "nucgss", "protein", "gene"],
                         "xml": ["genome", "nuccore", "nucgss", "nucest", "gene"],
                         "gb": ["genome", "nuccore", "nucgss"],
                         "est": ["nucest"],
                         "gss": ["nucgss"],
                         "acc": ["nuccore"],
                         }

    valid_formats = tuple(allowed_databases.keys())
    valid_databases = ("genome", "nuccore", "nucest",
                       "nucgss", "protein", "gene")

    p.add_option("--noversion", dest="noversion",
                 default=False, action="store_true",
                 help="Remove trailing accession versions")
    p.add_option("--format", default="fasta", choices=valid_formats,
                 help="download format [default: %default]")
    p.add_option("--database", default="nuccore", choices=valid_databases,
                 help="search database [default: %default]")
    p.add_option("--retmax", default=1000000, type="int",
                 help="how many results to return [default: %default]")
    p.add_option("--skipcheck", default=False, action="store_true",
                 help="turn off prompt to check file existence [default: %default]")
    p.add_option("--batchsize", default=500, type="int",
                 help="download the results in batch for speed-up [default: %default]")
    p.set_outdir(outdir=None)
    p.add_option("--outprefix", default="out",
                 help="output file name prefix [default: %default]")
    p.set_email()
    opts, args = p.parse_args(args)

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

    filename, = args
    if op.exists(filename):
        pf = filename.rsplit(".", 1)[0]
        list_of_terms = [row.strip() for row in open(filename)]
        if opts.noversion:
            list_of_terms = [x.rsplit(".", 1)[0] for x in list_of_terms]
    else:
        pf = filename
        # the filename is the search term
        list_of_terms = [filename.strip()]

    fmt = opts.format
    database = opts.database
    batchsize = opts.batchsize

    assert database in allowed_databases[fmt], \
        "For output format '{0}', allowed databases are: {1}".\
        format(fmt, allowed_databases[fmt])
    assert batchsize >= 1, "batchsize must >= 1"

    if " " in pf:
        pf = opts.outprefix

    outfile = "{0}.{1}".format(pf, fmt)

    outdir = opts.outdir
    if outdir:
        mkdir(outdir)

    # If noprompt, will not check file existence
    if not outdir:
        fw = must_open(outfile, "w", checkexists=True,
                       skipcheck=opts.skipcheck)
        if fw is None:
            return

    seen = set()
    totalsize = 0
    for id, size, term, handle in batch_entrez(list_of_terms, retmax=opts.retmax,
                                               rettype=fmt, db=database, batchsize=batchsize,
                                               email=opts.email):
        if outdir:
            outfile = urljoin(outdir, "{0}.{1}".format(term, fmt))
            fw = must_open(outfile, "w", checkexists=True,
                           skipcheck=opts.skipcheck)
            if fw is None:
                continue

        rec = handle.read()
        if id in seen:
            logging.error("Duplicate key ({0}) found".format(rec))
            continue

        totalsize += size
        print(rec, file=fw)
        print(file=fw)

        seen.add(id)

    if seen:
        print("A total of {0} {1} records downloaded.".
              format(totalsize, fmt.upper()), file=sys.stderr)

    return outfile
1024169
def setup_jobs(outpath, options, input_files):
    ''' 
    '''

    job_inputs = None
    reverse_mapping = None
    fasta_file_contents = {}

    # Generate FASTA files for PDB inputs
    # fasta_file_contents is a mapping from a file path to a pair (FASTA contents, file type). We remember the file type
    # since we offset residue IDs depending on file type i.e. for FASTA files, we treat each sequence separately and do
    # not renumber the fragments in postprocessing. For PDB files, however, we need to respect the order and length of
    # sequences so that we renumber the fragments appropriately in postprocessing - we assume that if a PDB file is passed in
    # then all chains (protein, RNA, or DNA) will be used in a Rosetta run.
    for input_file in input_files:
        assert(not(fasta_file_contents.get(input_file)))
        if any(fnmatch(input_file, x) for x in pdb_file_wildcards):
            pdb = PDB.from_filepath(input_file, strict=True)
            pdb.pdb_id = os.path.basename(input_file).split('.')[0]            
            if pdb.pdb_id.startswith('pdb') and len(pdb.pdb_id) >= 7:
                # Hack to rename FASTA identifiers for pdb*.ent files which are present in mirrors of the PDB
                pdb.pdb_id = pdb.pdb_id.replace('pdb', '')    
            fasta_file_contents[input_file] = (pdb.create_fasta(prefer_seqres_order = False), 'PDB')
        else:
            fasta_file_contents[input_file] = (read_file(input_file), 'FASTA')

    # Extract sequences from the input FASTA files.
    found_sequences, reverse_mapping, errors = get_sequences(options, fasta_file_contents)
    if found_sequences:
        reformat(found_sequences)
    if errors:
        return None, False, errors

    # Discard sequences that are the wrong chain.
    desired_sequences = {}
    for key, sequence in found_sequences.iteritems():
        pdb_id, chain, file_name = key
        if options.chain is None or chain == options.chain:
            desired_sequences[key] = sequence

    # Create the input FASTA and script files.
    job_inputs, errors = create_inputs(options, outpath, desired_sequences)

    # Create the reverse mapping file
    if reverse_mapping:
        segment_mapping_file = os.path.join(outpath, "segment_map.json")
        colorprinter.message("Creating a reverse mapping file %s." % segment_mapping_file)
        write_file(segment_mapping_file, json.dumps(reverse_mapping))

    # Create the post-processing script file
    post_processing_script = read_file(os.path.join(os.path.split(os.path.realpath(__file__))[0], 'post_processing.py'))
    write_file(os.path.join(outpath, 'post_processing.py'), post_processing_script, 'w')

    # Create the secondary structure filter file
    if options.secondary_structure_file:
        write_file(os.path.join(outpath, 'ss_filter.json'), json.dumps({'secondary_structure_filter' : SecondaryStructureDefinition.from_filepath(options.secondary_structure_file).data}), 'w')

    return job_inputs, reverse_mapping != None, errors
1052698
def df_proteins(self):
        """"""
        pre_df = []
        df_cols = ['gene', 'id', 'sequences', 'num_sequences', 'representative_sequence',
                   'repseq_gene_name', 'repseq_uniprot', 'repseq_description',
                   'num_structures', 'experimental_structures', 'num_experimental_structures',
                   'homology_models', 'num_homology_models',
                   'representative_structure', 'representative_chain', 'representative_chain_seq_coverage',
                   'repstruct_description', 'repstruct_resolution',
                   'num_sequence_alignments', 'num_structure_alignments']

        for g in self.genes:
            # Get per protein DataFrame
            protein_dict = g.protein.protein_statistics
            protein_dict['gene'] = g.id
            pre_df.append(protein_dict)

        df = pd.DataFrame.from_records(pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)