We assume a Rosetta-generated structure where residues are uniquely identified by number.
765816
def getAtomLinesForResidueInRosettaStructure(self, resid):
        ''''''
        lines = [line for line in self.lines if line[0:4] == "ATOM" and resid == int(line[22:27])]
        if not lines:
            #print('Failed searching for residue %d.' % resid)
            #print("".join([line for line in self.lines if line[0:4] == "ATOM"]))
            raise Exception("Could not find the ATOM/HETATM line corresponding to residue '%(resid)s'." % vars())
        return lines
364155
def cylinder_and_strand(self):
        '''

        '''
        top = self.topology
        # We build a  mini-state machine to find the
        # start end of helices and such
        in_helix = False
        helices_starts = []
        helices_ends = []
        coils = []

        coil = []
        for i, typ in enumerate(top['secondary_structure']):
            if typ == 'H':
                if in_helix == False:
                    # We become helices
                    helices_starts.append(top['residue_indices'][i][0])
                    in_helix = True

                    # We end the previous coil
                    coil.append(top['residue_indices'][i][0])
            else:
                if in_helix == True:
                    # We stop being helices
                    helices_ends.append(top['residue_indices'][i][0])

                    # We start a new coil
                    coil = []
                    coils.append(coil)
                    in_helix = False

                # We add control points
                coil.append(top['residue_indices'][i][0])
                [coil.append(j) for j in top['residue_indices'][i] if top['atom_names'][j] == 'CA']

        # We add the coils
        coil_representations = []
        for control_points in coils:
            rid = self.add_representation('smoothtube', {'coordinates': self.coordinates[control_points],
                                                   'radius': 0.05,
                                                   'resolution': 4,
                                                   'color': 0xffffff})
            coil_representations.append(rid)


        start_idx, end_idx = helices_starts, helices_ends
        cylinders = self.add_representation('cylinders', {'startCoords': self.coordinates[list(start_idx)],
                                              'endCoords': self.coordinates[list(end_idx)],
                                              'colors': [0xffff00] * len(self.coordinates),
                                              'radii': [0.15] * len(self.coordinates)})
        def update(self=self, cylinders=cylinders, coils=coils,
                   coil_representations=coil_representations,
                   start_idx=start_idx, end_idx=end_idx, control_points=control_points):
            for i, control_points in enumerate(coils):
                rid = self.update_representation(coil_representations[i],
                                                 {'coordinates': self.coordinates[control_points]})

            self.update_representation(cylinders, {'startCoords': self.coordinates[list(start_idx)],
                                                  'endCoords': self.coordinates[list(end_idx)]})

        self.update_callbacks.append(update)
        self.autozoom(self.coordinates)
150800
def _get_SEQRES_sequences(self):
        ''''''

        pdb_id = self.get_pdb_id()
        SEQRES_lines = self.parsed_lines["SEQRES"]

        modified_residue_mapping_3 = self.modified_residue_mapping_3
        # I commented this out since we do not need it for my current test cases
        #for k, v in self.modified_residues.iteritems():
        #    assert(v['modified_residue'] not in modified_residues)
        #    modified_residues[v['modified_residue']] = v['original_residue_3']

        for x in range(0, len(SEQRES_lines)):
            assert(SEQRES_lines[x][7:10].strip().isdigit())

        if not SEQRES_lines:
            #colortext.warning("WARNING: No SEQRES records were found. Kyle is trying to handle this gracefully, but Shane may need to fix it")
            return

        seqres_chain_order = []
        SEQRES_lines = [line[11:].rstrip() for line in SEQRES_lines] # we cannot strip the left characters as some cases e.g. 2MBP are missing chain identifiers

        # Collect all residues for all chains, remembering the chain order
        chain_tokens = {}
        for line in SEQRES_lines:
            chainID = line[0]
            if missing_chain_ids.get(self.pdb_id):
                chainID = missing_chain_ids[self.pdb_id]
            if chainID not in seqres_chain_order:
                seqres_chain_order.append(chainID)
            chain_tokens[chainID] = chain_tokens.get(chainID, [])
            chain_tokens[chainID].extend(line[6:].strip().split())

        sequences = {}
        self.chain_types = {}

        canonical_acid_types = set(residue_type_3to1_map.keys())
        canonical_acid_types.remove('UNK')
        assert(len(canonical_acid_types) == 20)

        for chain_id, tokens in chain_tokens.iteritems():

            # Determine whether the chain is DNA, RNA, or a protein chain
            # 1H38 is a good test for this - it contains DNA (chains E and G and repeated by H, K, N, J, M, P), RNA (chain F, repeated by I, L, O) and protein (chain D, repeated by A,B,C) sequences
            # 1ZC8 is similar but also has examples of DU
            # 4IHY has examples of DI (I is inosine)
            # 2GRB has RNA examples of I and U
            # 1LRP has protein chains with only CA atoms
            # This will throw an exception when a non-canonical is found which is not listed in basics.py. In that case, the list in basics.py should be updated.

            chain_type = None
            set_of_tokens = set(tokens)
            if (set_of_tokens.union(all_recognized_dna) == all_recognized_dna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'DNA'
            elif (set_of_tokens.union(all_recognized_rna) == all_recognized_rna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'RNA'
            elif len(set_of_tokens) == 1 and 'UNK' in set_of_tokens:
                chain_type = 'Unknown'
            elif not(set_of_tokens.intersection(canonical_acid_types)):
                # Zero canonical residues may imply a ligand or a heterogen chain
                chain_type = PDB._determine_heterogen_chain_type(set_of_tokens)
            else:
                assert(len(set(tokens).intersection(dna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                chain_type = 'Protein'
                if not self.chain_atoms.get(chain_id):
                    # possible for biological unit files
                    continue
                if self.chain_atoms[chain_id] == set(['CA']):
                    chain_type = 'Protein skeleton'

            # Get the sequence, mapping non-canonicals to the appropriate letter
            self.chain_types[chain_id] = chain_type

            sequence = []
            if chain_type == 'DNA':
                for r in tokens:
                    if dna_nucleotides_2to1_map.get(r):
                        sequence.append(dna_nucleotides_2to1_map[r])
                    else:
                        if non_canonical_dna.get(r):
                            sequence.append(non_canonical_dna[r])
                        else:
                            raise Exception("Unknown DNA residue %s." % r)
            elif chain_type == 'RNA':
                for r in tokens:
                    if r in rna_nucleotides:
                        sequence.append(r)
                    else:
                        if non_canonical_rna.get(r):
                            sequence.append(non_canonical_rna[r])
                        else:
                            raise Exception("Unknown RNA residue %s." % r)
            else:
                token_counter = 0
                for r in tokens:
                    token_counter += 1
                    if residue_type_3to1_map.get(r):
                        sequence.append(residue_type_3to1_map[r])
                    else:

                        if self.modified_residue_mapping_3.get(r):
                            sequence.append(residue_type_3to1_map[self.modified_residue_mapping_3.get(r)])
                        elif non_canonical_amino_acids.get(r):
                            #print('Mapping non-canonical residue %s to %s.' % (r, non_canonical_amino_acids[r]))
                            #print(SEQRES_lines)
                            #print(line)
                            sequence.append(non_canonical_amino_acids[r])
                        elif r == 'UNK':
                            continue
                        # Skip these residues
                        elif r == 'ACE' and token_counter == 1:
                            # Always allow ACE as the first residue of a chain
                            sequence.append('X')
                        elif r == 'ACE' and pdb_id in cases_with_ACE_residues_we_can_ignore:
                            sequence.append('X')
                            #continue
                        # End of skipped residues
                        else:
                            #print(modified_residue_mapping_3)
                            if modified_residue_mapping_3.get(r):
                                if modified_residue_mapping_3[r] == 'UNK':
                                    sequence.append('X')
                                else:
                                    assert(modified_residue_mapping_3[r] in residue_types_3)
                                    sequence.append(residue_type_3to1_map[modified_residue_mapping_3[r]])
                            else:
                                raise Exception("Unknown protein residue %s in chain %s." % (r, chain_id))
            sequences[chain_id] = "".join(sequence)

        self.seqres_chain_order = seqres_chain_order

        # Create Sequence objects for the SEQRES sequences
        for chain_id, sequence in sequences.iteritems():
            self.seqres_sequences[chain_id] = Sequence.from_sequence(chain_id, sequence, self.chain_types[chain_id])
930542
public String getProteinSequenceString() {

		if (sequence != null)
			return sequence.toString();

		StringBuilder builder = new StringBuilder();
		for (Atom a : reprAtoms)
			// This method preferred over getChemComp.getOneLetterCode because
			// it returns always X for Unknown residues
			builder.append(StructureTools.get1LetterCode(a.getGroup()
					.getPDBName()));

		return builder.toString();
	}
1096227
def _create_sequence_maps(self):
        ''''''

        if self.sifts:
            self.sifts_atom_to_seqres_sequence_maps = self.sifts.atom_to_seqres_sequence_maps
            self.sifts_seqres_to_uniparc_sequence_maps = self.sifts.seqres_to_uniparc_sequence_maps
            self.sifts_atom_to_uniparc_sequence_maps = self.sifts.atom_to_uniparc_sequence_maps
            if self.pdb_id in pdbs_with_do_not_use_SIFTS_for_these_chains:
                for chain_id in self.sifts_atom_to_seqres_sequence_maps.keys() + self.sifts_seqres_to_uniparc_sequence_maps.keys() + self.sifts_atom_to_uniparc_sequence_maps.keys():
                    if (self.pdb_id, chain_id) in do_not_use_SIFTS_for_these_chains:
                        self.sifts_atom_to_seqres_sequence_maps[chain_id] = SequenceMap()
                        self.sifts_seqres_to_uniparc_sequence_maps = SequenceMap()
                        self.sifts_atom_to_uniparc_sequence_maps = SequenceMap()

        if self.pdb_to_rosetta_residue_map_error:
            self.rosetta_to_atom_sequence_maps = {}
            for c in self.atom_sequences.keys():
                self.rosetta_to_atom_sequence_maps[c] = SequenceMap()
        else:
            self.rosetta_to_atom_sequence_maps = self.pdb.rosetta_to_atom_sequence_maps

        # If we removed atoms from the PDB file, we need to remove them from the maps so that our validations hold later on
        self.pdbml_atom_to_seqres_sequence_maps = self.pdbml.atom_to_seqres_sequence_maps
        if self.pdb_id in ROSETTA_HACKS_residues_to_remove:
            for residue_to_remove in ROSETTA_HACKS_residues_to_remove[self.pdb_id]:
                chain_id = residue_to_remove[0]
                self.pdbml_atom_to_seqres_sequence_maps[chain_id].remove(residue_to_remove)
                #if self.sifts:
                #    self.sifts_atom_to_seqres_sequence_maps[chain_id].remove(residue_to_remove)

        if self.pdb_id not in do_not_use_the_sequence_aligner:
            self.clustal_seqres_to_uniparc_sequence_maps = self.PDB_UniParc_SA.seqres_to_uniparc_sequence_maps
70614
def invariants(mol):
    """"""
    atom_ids = {}
    for a in mol.atoms:
        components = []
        components.append(a.number)
        components.append(len(a.oatoms))
        components.append(a.hcount)
        components.append(a.charge)
        components.append(a.mass)
        if len(a.rings) > 0:
            components.append(1)

        atom_ids[a.index] = gen_hash(components)

    return atom_ids
749747
def num_mutations(self):
        '''
        
        '''
        self._initialize()
        res = super(Container, self).num_mutations()
        return res
437854
def unique_id(self):
        
        chain = self.ampal_parent.ampal_parent.id
        residue = self.ampal_parent.id
        return chain, residue, self.id
673884
def build(self):
        """"""
        ang_per_res = (2 * numpy.pi) / self.residues_per_turn
        atom_offsets = _atom_offsets[self.helix_type]
        if self.handedness == 'l':
            handedness = -1
        else:
            handedness = 1

        atom_labels = ['N', 'CA', 'C', 'O']
        if all([x in atom_offsets.keys() for x in atom_labels]):
            res_label = 'GLY'
        else:
            res_label = 'UNK'

        monomers = []
        for i in range(self.num_monomers):
            residue = Residue(mol_code=res_label, ampal_parent=self)
            atoms_dict = OrderedDict()
            for atom_label in atom_labels:
                r, zeta, z_shift = atom_offsets[atom_label]
                rot_ang = ((i * ang_per_res) + zeta) * handedness
                z = (self.rise_per_residue * i) + z_shift
                coords = cylindrical_to_cartesian(
                    radius=r, azimuth=rot_ang, z=z, radians=True)
                atom = Atom(
                    coordinates=coords, element=atom_label[0],
                    ampal_parent=residue, res_label=atom_label)
                atoms_dict[atom_label] = atom
            residue.atoms = atoms_dict
            monomers.append(residue)
        self._monomers = monomers
        self.relabel_monomers()
        self.relabel_atoms()
        return
20343
def _topology_from_parmed(structure, non_element_types):
    """"""
    topology = app.Topology()
    residues = dict()
    for pmd_residue in structure.residues:
        chain = topology.addChain()
        omm_residue = topology.addResidue(pmd_residue.name, chain)
        residues[pmd_residue] = omm_residue
    atoms = dict()  # pmd.Atom: omm.Atom

    for pmd_atom in structure.atoms:
        name = pmd_atom.name
        if pmd_atom.name in non_element_types:
            element = non_element_types[pmd_atom.name]
        else:
            if (isinstance(pmd_atom.atomic_number, int) and
                    pmd_atom.atomic_number != 0):
                element = elem.Element.getByAtomicNumber(pmd_atom.atomic_number)
            else:
                element = elem.Element.getBySymbol(pmd_atom.name)

        omm_atom = topology.addAtom(name, element, residues[pmd_atom.residue])
        atoms[pmd_atom] = omm_atom
        omm_atom.bond_partners = []

    for bond in structure.bonds:
        atom1 = atoms[bond.atom1]
        atom2 = atoms[bond.atom2]
        topology.addBond(atom1, atom2)
        atom1.bond_partners.append(atom2)
        atom2.bond_partners.append(atom1)
    if structure.box_vectors and np.any([x._value for x in structure.box_vectors]):
        topology.setPeriodicBoxVectors(structure.box_vectors)

    positions = structure.positions
    return topology, positions
1322607
def get_differing_atom_residue_ids(self, pdb_name, pdb_list = []):
        ''''''

        # partition_by_sequence is a map pdb_name -> Int where two pdb names in the same equivalence class map to the same integer (i.e. it is a partition)
        # representative_pdbs is a map Int -> pdb_object mapping the equivalence classes (represented by an integer) to a representative PDB file
        # self.pdb_name_to_structure_mapping : pdb_name -> pdb_object

        # Sanity checks
        assert(pdb_name in self.pdb_names)
        assert(set(pdb_list).intersection(set(self.pdb_names)) == set(pdb_list)) # the names in pdb_list must be in pdb_names

        # 1. Get the representative structure for pdb_name
        representative_pdb_id = self.partition_by_sequence[pdb_name]
        representative_pdb = self.representative_pdbs[representative_pdb_id]

        # 2. Get the other representative structures as dictated by pdb_list
        other_representative_pdbs = set()
        other_representative_pdb_ids = set()
        if not pdb_list:
            pdb_list = self.pdb_names
        for opdb_name in pdb_list:
            orepresentative_pdb_id = self.partition_by_sequence[opdb_name]
            other_representative_pdb_ids.add(orepresentative_pdb_id)
            other_representative_pdbs.add(self.representative_pdbs[orepresentative_pdb_id])
        other_representative_pdbs.discard(representative_pdb)
        other_representative_pdb_ids.discard(representative_pdb_id)

        # Early out if pdb_list was empty (or all pdbs were in the same equivalence class)
        if not other_representative_pdbs:
            return []

        # 3. Return all residues of pdb_name's representative which differ from all the other representatives
        differing_atom_residue_ids = set()
        for other_representative_pdb_id in other_representative_pdb_ids:
            differing_atom_residue_ids = differing_atom_residue_ids.union(set(self.differing_atom_residue_ids[(representative_pdb_id, other_representative_pdb_id)]))

        return sorted(differing_atom_residue_ids)
392289
private int getGroupTypIndicator(PolymerType polymerType) {
		if(PolymerType.PROTEIN_ONLY.contains(polymerType)){
			return 1;
		}
		if(PolymerType.POLYNUCLEOTIDE_ONLY.contains(polymerType)){
			return 2;
		}
		return 0;
	}
1094214
def protein_subsequences_around_mutations(effects, padding_around_mutation):
    
    protein_subsequences = {}
    protein_subsequence_start_offsets = {}
    for effect in effects:
        protein_sequence = effect.mutant_protein_sequence
        # some effects will lack a mutant protein sequence since
        # they are either silent or unpredictable
        if protein_sequence:
            mutation_start = effect.aa_mutation_start_offset
            mutation_end = effect.aa_mutation_end_offset
            seq_start_offset = max(
                0,
                mutation_start - padding_around_mutation)
            # some pseudogenes have stop codons in the reference sequence,
            # if we try to use them for epitope prediction we should trim
            # the sequence to not include the stop character '*'
            first_stop_codon_index = protein_sequence.find("*")
            if first_stop_codon_index < 0:
                first_stop_codon_index = len(protein_sequence)

            seq_end_offset = min(
                first_stop_codon_index,
                mutation_end + padding_around_mutation)
            subsequence = protein_sequence[seq_start_offset:seq_end_offset]
            protein_subsequences[effect] = subsequence
            protein_subsequence_start_offsets[effect] = seq_start_offset
    return protein_subsequences, protein_subsequence_start_offsets
371009
def get_generic_subseq_3D(protein, cutoff, prop, condition):
        """"""
        if not protein.representative_structure:
            log.error('{}: no representative structure, cannot search for subseq'.format(protein.id))
            return {'subseq_len': 0, 'subseq': None, 'subseq_resnums': []}

        subseq, subseq_resnums = protein.get_seqprop_subsequence_from_structchain_property(property_key=prop,
                                                                                           property_value=cutoff,
                                                                                           condition=condition,
                                                                                           use_representatives=True,
                                                                                           return_resnums=True) or (
                                 None, [])

        return {'subseq_len': len(subseq_resnums), 'subseq': subseq, 'subseq_resnums': subseq_resnums}
622301
private SubstructureIdentifier toCanonical() {
		StructureIdentifier real = getStructureIdentifier();
		if(real != null) {
			try {
				return real.toCanonical();
			} catch (StructureException e) {
				// generate fake one if needed
			}
		}

		// No identifier set, so generate based on residues present in the structure
		List<ResidueRange> range = new ArrayList<>();
		for (Chain chain : getChains()) {
			List<Group> groups = chain.getAtomGroups();
			ListIterator<Group> groupsIt = groups.listIterator();
			if(!groupsIt.hasNext()) {
				continue; // no groups in chain
			}
			Group g = groupsIt.next();
			ResidueNumber first = g.getResidueNumber();

			//TODO Detect missing intermediate residues -sbliven, 2015-01-28
			//Already better than previous whole-chain representation

			// get last residue
			while(groupsIt.hasNext()) {
				g = groupsIt.next();
			}
			ResidueNumber last = g.getResidueNumber();

			range.add(new ResidueRange(chain.getName(),first,last));
		}
		return new SubstructureIdentifier(getPDBCode(),range);
	}
1064025
private void assembleLinkages(List<List<Atom[]>> matchedAtomsOfLinkages,
			ProteinModification mod, List<ModifiedCompound> ret) {
		ModificationCondition condition = mod.getCondition();
		List<ModificationLinkage> modLinks = condition.getLinkages();

		int nLink = matchedAtomsOfLinkages.size();
		int[] indices = new int[nLink];
		Set<ModifiedCompound> identifiedCompounds = new HashSet<ModifiedCompound>();
		while (indices[0]<matchedAtomsOfLinkages.get(0).size()) {
			List<Atom[]> atomLinkages = new ArrayList<Atom[]>(nLink);
			for (int iLink=0; iLink<nLink; iLink++) {
				Atom[] atoms = matchedAtomsOfLinkages.get(iLink).get(indices[iLink]);
				atomLinkages.add(atoms);
			}
			if (matchLinkages(modLinks, atomLinkages)) {
				// matched

				int n = atomLinkages.size();
				List<StructureAtomLinkage> linkages = new ArrayList<StructureAtomLinkage>(n);
				for (int i=0; i<n; i++) {
					Atom[] linkage = atomLinkages.get(i);
					StructureAtomLinkage link = StructureUtil.getStructureAtomLinkage(
							linkage[0], residues.contains(linkage[0].getGroup()),
							linkage[1], residues.contains(linkage[1].getGroup()));
					linkages.add(link);
				}

				ModifiedCompound mc = new ModifiedCompoundImpl(mod, linkages);
				if (!identifiedCompounds.contains(mc)) {
					ret.add(mc);
					identifiedCompounds.add(mc);
				}
			}

			// indices++ (e.g. [0,0,1]=>[0,0,2]=>[1,2,0])
			int i = nLink-1;
			while (i>=0) {
				if (i==0 || indices[i]<matchedAtomsOfLinkages.get(i).size()-1) {
					indices[i]++;
					break;
				} else {
					indices[i] = 0;
					i--;
				}
			}
		}
	}
524291
def get_protein_modification_language(identifier_qualified: ParserElement) -> ParserElement:
    """"""
    pmod_identifier = MatchFirst([
        identifier_qualified,
        pmod_default_ns,
        pmod_legacy_ns
    ])

    return pmod_tag + nest(
        Group(pmod_identifier)(IDENTIFIER) +
        Optional(
            WCW +
            amino_acid(PMOD_CODE) +
            Optional(
                WCW +
                ppc.integer(PMOD_POSITION)
            )
        )
    )