Loops peptides from multiple files, fetches PSMs from
    sequence:PSM map, outputs correctly PSM mapped peptides
770085
def merge_peptides(fns, ns):
    """"""
    peptides_to_map = reader.generate_peptides_multiple_fractions(fns, ns)
    psmmap = create_merge_psm_map(peptides_to_map, ns)
    peptides = reader.generate_peptides_multiple_fractions(fns, ns)
    for peptide in peptides:
        seq = reader.get_peptide_seq(peptide, ns)
        psm_ids = reader.get_psm_ids_from_peptide(peptide, ns)
        # remove current psm ids, repopulate with stored ones
        psm_ids.clear()
        for new_psm_id in psmmap[seq]:
            etree.SubElement(psm_ids, 'psm_id').text = new_psm_id
        yield formatting.string_and_clear(peptide, ns)
188179
def create_psm_lookup(fn, fastafn, mapfn, header, pgdb, unroll=False,
                      specfncol=None, decoy=False,
                      fastadelim=None, genefield=None):
    
    proteins = store_proteins_descriptions(pgdb, fastafn, fn, mapfn, header,
                                           decoy, fastadelim, genefield)
    mzmlmap = pgdb.get_mzmlfile_map()
    sequences = {}
    for psm in tsvreader.generate_tsv_psms(fn, header):
        seq = tsvreader.get_psm_sequence(psm, unroll)
        sequences[seq] = 1
    pgdb.store_pepseqs(((seq,) for seq in sequences))
    pepseqmap = pgdb.get_peptide_seq_map()
    psms = []
    for row, psm in enumerate(tsvreader.generate_tsv_psms(fn, header)):
        specfn, psm_id, scan, seq, score = tsvreader.get_psm(psm, unroll,
                                                             specfncol)
        if len(psms) % DB_STORE_CHUNK == 0:
            pgdb.store_psms(psms)
            psms = []
        psms.append({'rownr': row,
                     'psm_id': psm_id,
                     'seq': pepseqmap[seq],
                     'score': score,
                     'specfn': mzmlmap[specfn],
                     'scannr': scan,
                     'spec_id': '{}_{}'.format(mzmlmap[specfn], scan),
                     })
    pgdb.store_psms(psms)
    pgdb.index_psms()
    store_psm_protein_relations(fn, header, pgdb, proteins)
919986
def get_proteins_psms_for_map(self):
        """"""
        fields = ['p.gene_acc', 'sets.set_name',
                  'pep.sequence', 'psm.psm_id']
        firstjoin = ('protein_psm', 'pp', 'protein_acc')
        genetable = 'genes'
        return self.get_unique_gene_psms(genetable, fields, firstjoin)
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
664633
def get_unique_gene_psms(self, genetable, fields, firstjoin):
        """"""
        lastgene = None
        gpsms_out, gp_ids = [], []
        for gpsm in self.get_proteins_psms(genetable, fields, firstjoin):
            if gpsm[0] != lastgene:
                for outpsm in gpsms_out:
                    yield outpsm
                lastgene = gpsm[0]
                gpsms_out, gp_ids = [], []
            gp_id = gpsm[0] + gpsm[1] + gpsm[3]
            if gp_id not in gp_ids:
                gp_ids.append(gp_id)
                gpsms_out.append(gpsm)
        for outpsm in gpsms_out:
            yield outpsm
759627
def get_proteins_from_psm(line):
    """"""
    proteins = line[mzidtsvdata.HEADER_PROTEIN].split(';')
    outproteins = []
    for protein in proteins:
        prepost_protein = re.sub('\(pre=.*post=.*\)', '', protein).strip()
        outproteins.append(prepost_protein)
    return outproteins
458327
def map_all(self, prot_alignment, nucl_sequences):
        
        zipped = itertools.zip_longest(prot_alignment, nucl_sequences)
        for p, n in zipped:
            if p is None:
                raise ValueError("Exhausted protein sequences")
            elif n is None:
                raise ValueError("Exhausted nucleotide sequences")
            yield self.map_alignment(p, n)
122079
def epcr_parse(self):
        
        logging.info('Parsing ePCR outputs')
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                # Create a set to store all the unique results
                toxin_set = set()
                if os.path.isfile(sample[self.analysistype].resultsfile):
                    with open(sample[self.analysistype].resultsfile) as epcrresults:
                        for result in epcrresults:
                            # Only the lines without a # contain results
                            if "#" not in result:
                                # Split on \t
                                data = result.split('\t')
                                # The subtyping primer pair is the first entry on lines with results
                                vttype = data[0].split('_')[0]
                                # Add the verotoxin subtype to the set of detected subtypes
                                toxin_set.add(vttype)
                # Create a string of the entries in the sorted list of toxins joined with ";"
                sample[self.analysistype].toxinprofile = ";".join(sorted(list(toxin_set))) if toxin_set else 'ND'
            else:
                setattr(sample, self.analysistype, GenObject())
                sample[self.analysistype].toxinprofile = 'NA'
1331152
def prune_structure_according_to_loops_file(self, loops_file_content, search_radius, expected_loop_length = None, generate_pymol_session = True):
        ''''''
        lf = LoopsFile(loops_file_content)
        assert(len(lf.data) == 1) # todo: remove

        # Extract the list of Residues defined in the loops file
        loop_residues = []
        for loop in lf.data:
            # Identify the Residues corresponding to the loops file definition
            start_id = str(loop['start'])
            end_id = str(loop['end'])
            start = []
            end = []
            for chain, chain_residues in self.residues.iteritems():
                for res, residue_object in chain_residues.iteritems():
                    if res.strip() == start_id:
                        start.append((chain, res, residue_object))
                    if res.strip() == end_id:
                        end.append((chain, res, residue_object))
            if len(start) != len(end) != 1:
                raise Exception('The PDB is ambiguous with respect to the loops file i.e. more than one PDB residue corresponds to the start or end residue defined in the loops file.')
            start, end = start[0], end[0]

            # We assume that the loops file identifies a range in document order
            loop_start = False
            for l in self.indexed_lines:
                if l[0] == 'ATOM' or l[0] == 'HETATM':
                    atom_residue = l[3].residue
                    if not loop_start:
                        if atom_residue == start[2]:
                            loop_start = True
                            if atom_residue not in loop_residues:
                                loop_residues.append(atom_residue)
                    else:
                        if atom_residue not in loop_residues:
                            loop_residues.append(atom_residue)
                        if atom_residue == end[2]:
                            break
        if expected_loop_length != None and expected_loop_length != len(loop_residues):
            raise Exception('Expected to identify {0} residues but {1} were identified.'.format(expected_loop_length, len(loop_residues)))

        # Determine the sidechain atoms to be removed
        sidechain_atom_serial_numbers = self.find_sidechain_atoms_within_radius_of_residue_objects(loop_residues, search_radius)

        # Determine the loop residue atoms to be removed
        loop_residues_ids = set()
        loop_atom_serial_numbers = set()
        for loop_residue in loop_residues:
            for aatom in loop_residue.get('ATOM'):
                loop_atom_serial_numbers.add(aatom.serial_number)
        assert(len(sidechain_atom_serial_numbers.intersection(loop_atom_serial_numbers)) == 0)

        return self.prune(loop_atom_serial_numbers, sidechain_atom_serial_numbers, generate_pymol_session = generate_pymol_session)
420755
def generate_psms_split(fn, oldheader, bioset, splitcol):
    """"""
    try:
        splitcolnr = get_splitcolnr(oldheader, bioset, splitcol)
    except IndexError:
        raise RuntimeError('Cannot find bioset header column in '
                           'input file {}, though --bioset has '
                           'been passed'.format(fn))
    for psm in tsvreader.generate_tsv_psms(fn, oldheader):
        yield {'psm': psm,
               'split_pool': psm[oldheader[splitcolnr]]
               }
7395
def create_pep_protein_quant_lookup(fns, pqdb, poolnames, featcolnr, patterns,
                                    storefuns, isobqcolpattern=None,
                                    psmnrpattern=None):
    """"""
    tablefn_map = create_tablefn_map(fns, pqdb, poolnames)
    feat_map = pqdb.get_feature_map()
    for pattern, storefun in zip(patterns, storefuns):
        if pattern is None:
            continue
        colmap = get_colmap(fns, pattern, single_col=True)
        if colmap:
            store_single_col_data(fns, tablefn_map, feat_map,
                                  storefun, featcolnr, colmap)
    if isobqcolpattern is not None:
        isocolmap = get_colmap(fns, isobqcolpattern, antipattern=psmnrpattern)
    else:
        return
    if psmnrpattern is not None:
        psmcolmap = get_colmap(fns, psmnrpattern)
    else:
        psmcolmap = False
    create_isobaric_quant_lookup(fns, tablefn_map,
                                 feat_map, pqdb,
                                 featcolnr,
                                 isocolmap, psmcolmap)
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
889831
def get_peptable_headerfields(headertypes, lookup=False, poolnames=False):
    """"""
    field_defs = {'isoquant': get_isoquant_fields,
                  'precursorquant': get_precursorquant_fields,
                  'peptidefdr': get_peptidefdr_fields,
                  'peptidepep': get_peptidepep_fields,
                  'proteindata': get_proteininfo_fields,
                  }
    return generate_headerfields(headertypes, field_defs, poolnames, lookup)
720076
def read_fastas(input_files):
    
    tumor_file = [y for x, y in input_files.items() if x.startswith('T')][0]
    normal_file = [y for x, y in input_files.items() if x.startswith('N')][0]
    output_files = defaultdict(list)
    output_files = _read_fasta(tumor_file, output_files)
    num_entries = len(output_files)
    output_files = _read_fasta(normal_file, output_files)
    assert len(output_files) == num_entries
    return output_files
1055630
def create_featuredata_map(pgdb, psm_fill_fun, pgene_fill_fun, genecentric=False, 
        count_fun=None, pool_to_output=False, get_uniques=False, is_peptides=False):
    
    if genecentric:
        pgcontentmap = None
    else:
        pgcontentmap = get_proteingroup_content(pgdb)
    proteindata = {}
    for pgenedata in pgdb.get_protein_gene_symbol_for_map():
        if is_peptides:
            # include sequence (index 1 of record) to fill func
            pgene_fill_fun(proteindata, pgenedata[1], pgenedata[0], 
                    pgenedata, genecentric, pgcontentmap)
        else:
            pgene_fill_fun(proteindata, pgenedata[0], 
                    pgenedata, genecentric, pgcontentmap)
    protein_psms_data = pgdb.get_proteins_psms_for_map()
    psmdata = next(protein_psms_data)
    last_prot, last_pool = psmdata[0], psmdata[1]
    psm_fill_fun(proteindata, last_prot, last_pool, psmdata)
    for psmdata in protein_psms_data:
        p_acc, samplepool = psmdata[0], psmdata[1]
        if pool_to_output and samplepool != pool_to_output:
            continue
        if samplepool != last_pool or p_acc != last_prot:
            if count_fun is not None:
                count_fun(proteindata, last_prot, last_pool)
            last_pool, last_prot = samplepool, p_acc
        psm_fill_fun(proteindata, p_acc, samplepool, psmdata)
    if count_fun is not None:
        count_fun(proteindata, last_prot, last_pool)
    
    if get_uniques:
        for no_uni, pool, pacc in pgdb.get_unique_peptide_nrs():
            proteindata[pacc]['pools'][pool]['unipeps'] = no_uni
    return proteindata
916640
def add_ms1_quant_from_top3_mzidtsv(proteins, psms, headerfields, protcol):
    """"""
    if not protcol:
        protcol = mzidtsvdata.HEADER_MASTER_PROT
    top_ms1_psms = generate_top_psms(psms, protcol)
    for protein in proteins:
        prot_acc = protein[prottabledata.HEADER_PROTEIN]
        prec_area = calculate_protein_precursor_quant(top_ms1_psms, prot_acc)
        outprotein = {k: v for k, v in protein.items()}
        outprotein[headerfields['precursorquant'][
            prottabledata.HEADER_AREA][None]] = str(prec_area)
        yield outprotein
14829
private Map<String, Set<SequenceEntityReference>> organismToProteinRefsMap(Set<SequenceEntityReference> seqErs)
	{
		Map<String,Set<SequenceEntityReference>> map = new HashMap<String, Set<SequenceEntityReference>>();

		if(seqErs.isEmpty())
			throw new IllegalArgumentException("Empty set");
		
		if (crossSpeciesCheckEnabled) {
			for (SequenceEntityReference r : seqErs) {
				String key = getOrganismKey(r.getOrganism()); // null org. is ok (key == "")
				//collect PRs only from allowed organisms
				if(allowedOrganisms==null || allowedOrganisms.isEmpty() || allowedOrganisms.contains(key)) {
					Set<SequenceEntityReference> sers = map.get(key);
					if (sers == null) {
						sers = new HashSet<SequenceEntityReference>();
						map.put(key, sers);
					}
					sers.add(r);
				}
			}
		} else {
			final Set<SequenceEntityReference> sers = new HashSet<SequenceEntityReference>();
			for (SequenceEntityReference r : seqErs) {
				String key = getOrganismKey(r.getOrganism());
				//collect PRs only from allowed organisms
				if(allowedOrganisms==null || allowedOrganisms.isEmpty() || allowedOrganisms.contains(key)) {
					sers.add(r);
				}
			}
			map.put("", sers);
		}
				
		return map;
	}