The residue numbers of the protein can be reformated in case of misallignment with the convention.
            Takes:
                 * offset * - a number that represents by how many residues the numbering has to be shifted.
973190
def renumber_system(self, offset):
        

        self.universe.protein = self.universe.select_atoms("protein")
        self.universe.protein.residues.resids = self.universe.protein.residues.resids+int(offset)
114591
def check_padding_around_mutation(given_padding, epitope_lengths):
    
    min_required_padding = max(epitope_lengths) - 1
    if not given_padding:
        return min_required_padding
    else:
        require_integer(given_padding, "Padding around mutation")
        if given_padding < min_required_padding:
            raise ValueError(
                "Padding around mutation %d cannot be less than %d "
                "for epitope lengths %s" % (
                    given_padding,
                    min_required_padding,
                    epitope_lengths))
        return given_padding
382566
def read_and_calculate(workspace, pdb_paths):
    

    # Parse the given restraints file.  The restraints definitions are used to
    # calculate the "restraint_dist" metric, which reflects how well each
    # structure achieves the desired geometry. Note that this is calculated
    # whether or not restraints were used to create the structures in question.
    # For example, the validation runs don't use restraints but the restraint
    # distance is a very important metric for deciding which designs worked.

    restraints = parse_restraints(workspace.restraints_path)

    # Calculate score and distance metrics for each structure.

    from klab.bio.basics import residue_type_3to1_map

    records = []
    metadata = {}
    num_restraints = len(restraints) + 1
    atom_xyzs = {}
    fragment_size = 0

    # It's kinda hard to tell which lines are part of the score table.  The 
    # first column has some pretty heterogeneous strings (examples below) and 
    # all the other columns are just numbers.  My strategy here is to try to 
    # make a regular expression that matches all of these examples, with the 
    # exception of the ligand.  I think the ligand will simply be too 
    # heterogeneous to match, and the purpose of this is to get dunbrack 
    # scores, which the ligand doesn't have.
    #
    #   MET:NtermProteinFull_1
    #   ASN_2
    #   LYS:protein_cutpoint_lower_39
    #   ASP:protein_cutpoint_upper_40
    #   ALA:CtermProteinFull_124
    #   HIS_D_224
    #   pdb_EQU_250

    score_table_pattern = re.compile(
            r'^[A-Z]{3}(?:_[A-Z])?'  # Residue name with optional tautomer.
            r'(?::[A-Za-z_]+)?'      # Optional patch type.
            r'_([0-9]+) '            # Residue number preceded by underscore.
    )                                # The terminal space is important to match
                                     # the full residue number.

    for i, path in enumerate(sorted(pdb_paths)):
        record = {'path': os.path.basename(path)}
        sequence = ""
        sequence_map = {}
        last_residue_id = None
        dunbrack_index = None
        dunbrack_scores = {}

        # Update the user on our progress, because this is often slow.

        sys.stdout.write("\rReading '{}' [{}/{}]".format(
            os.path.relpath(os.path.dirname(path)), i+1, len(pdb_paths)))
        sys.stdout.flush()

        # Read the PDB file, which we are assuming is gzipped.

        try:
            with gzip.open(path) as file:
                lines = file.readlines()
        except IOError:
            print "\nFailed to read '{}'".format(path)
            continue

        if not lines:
            print "\n{} is empty".format(path)
            continue

        # Get different information from different lines in the PDB file.  Some
        # of these lines are specific to different simulations.

        for line in lines:
            line = line.decode('utf8')

            score_table_match = \
                    dunbrack_index and score_table_pattern.match(line)

            if line.startswith('pose'):
                meta = ScoreMetadata(
                        name='total_score',
                        title='Total Score',
                        unit='REU',
                        order=1,
                )
                record['total_score'] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('label'):
                fields = line.split()
                dunbrack_index = fields.index('fa_dun')

            elif score_table_match:
                residue_id = int(score_table_match.group(1))
                for restraint in restraints:
                    if residue_id in restraint.residue_ids:
                        dunbrack_score = float(line.split()[dunbrack_index])
                        dunbrack_scores[residue_id] = dunbrack_score
                        break

            elif line.startswith('rmsd'):
                meta = ScoreMetadata(
                        name='loop_rmsd',
                        title='Loop RMSD (Backbone Heavy-Atom)',
                        unit='Å',
                        guide=1.0, lower=0.0, upper='95%', order=4,
                )
                record[meta.name] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('  all_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats',
                        title='Buried Unsatsified H-Bonds',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('  sc_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats_sidechain',
                        title='Buried Unsatisfied H-Bonds (Sidechain)',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('  bb_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats_backbone',
                        title='Buried Unsatisfied H-Bonds (Backbone)',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('time'):
                meta = ScoreMetadata(
                        name='simulation_time',
                        title='Simulation Time',
                        unit='sec',
                        order=5,
                )
                record[meta.name] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('FragmentScoreFilter '):
                fragment_size = line.split()[2].split('-')[0]

            elif line.startswith('FSF') or line.startswith('FragmentScoreFilter_metric'):
                splitline = line.split()
                if splitline[1] == 'Max':
                    max_res = 0
                    max_crmsd = 0
                    if splitline[3] == 'res:':
                        max_res = splitline[3]
                        meta = ScoreMetadata(
                                name='max_fragment_crmsd_position',
                                title = 'Max {}-Residue Fragment RMSD \
(C-Alpha) Position'.format(fragment_size),
                                order=7)
                    elif splitline[3] == 'score:':
                        max_crmsd = splitline[3]
                        meta = ScoreMetadata(
                                name='max_fragment_crmsd_score',
                                title = 'Max {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                                order=7)

                elif splitline[1] == 'Min':
                    min_res = 0
                    min_crmsd = 0
                    if splitline[3] == 'res:':
                        min_res = splitline[3]
                        meta = ScoreMetadata(
                                name='min_fragment_crmsd_position',
                                title = 'Min {}-Residue Fragment RMSD \
(C-Alpha) Position'.format(fragment_size),
                                order=8)
                    elif splitline[3] == 'score:':
                        min_crmsd = splitline[3]
                        meta = ScoreMetadata(
                                name='min_fragment_crmsd_score',
                                title = 'Min {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                                order=8)

                elif splitline[1] == 'Avg':
                    meta = ScoreMetadata(
                            name='avg_fragment_crmsd',
                            title='Avg {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                            order=9)
                else:
                    position = splitline[2]
                    crmsd = splitline[4]
                    meta = ScoreMetadata(
                            name='fragment_crmsd_pos_{}'.format(position),
                            title='{}-Residue Fragment RMSD at Res {} \
(C-Alpha)'.format(fragment_size,position),
                            order=6)

                record[meta.name] = float(splitline[4])
                metadata[meta.name] = meta

            elif line.startswith('EXTRA_SCORE'):
                tokens = line[len('EXTRA_SCORE_'):].rsplit(None, 1)
                meta = parse_extra_metric(tokens[0], 5)
                record[meta.name] = float(tokens[1])
                metadata[meta.name] = meta

            elif line.startswith('EXTRA_METRIC'):
                tokens = line[len('EXTRA_METRIC '):].rsplit(None, 1)

                # Ignore the BuriedUnsat filter.  It just reports 911 every 
                # time, and we extract the actual buried unsat information from 
                # some other lines it adds to the PDB.
                if tokens[0] == 'IGNORE':
                    continue
                if tokens[0] == 'Buried Unsatisfied H-Bonds [-|#]':
                    continue

                meta = parse_extra_metric(tokens[0], 5)
                record[meta.name] = float(tokens[1])
                metadata[meta.name] = meta

            elif (line.startswith('ATOM') or line.startswith('HETATM')):
                atom_name = line[12:16].strip()
                residue_id = int(line[22:26].strip())
                residue_name = line[17:20].strip()

                # Keep track of this model's sequence.

                if residue_id != last_residue_id:
                    if line.startswith('ATOM'): 
                        one_letter_code = residue_type_3to1_map.get(residue_name, 'X')
                        sequence += one_letter_code
                        sequence_map[residue_id] = one_letter_code
                        last_residue_id = residue_id
                    elif line.startswith('HETATM'):
                        sequence_map[residue_id] = 'X'
                        last_residue_id = residue_id

                # Save the coordinate for this atom.  This will be used later 
                # to calculate restraint distances.

                atom_xyzs[atom_name, residue_id] = xyz_to_array((
                        line[30:38], line[38:46], line[46:54]))

        # Calculate how well each restraint was satisfied.

        restraint_values = {}
        restraint_values_by_residue = {}
        is_sidechain_restraint = {}
        restraint_units = {
                'dist': 'Å',
                'angle': '°',
        }

        for restraint in restraints:
            d = restraint.distance_from_ideal(atom_xyzs)
            metric = restraint.metric
            backbone_atoms = set(['N', 'C', 'CA', 'O'])
            backbone_restraint = backbone_atoms.issuperset(restraint.atom_names)

            restraint_values.setdefault(metric, []).append(d)
            restraint_values_by_residue.setdefault(metric, {})
            for i in restraint.residue_ids:
                restraint_values_by_residue[metric].setdefault(i, []).append(d)
                is_sidechain_restraint[i] = (not backbone_restraint) \
                        or is_sidechain_restraint.get(i, False)

        for metric, values in restraint_values.items():
            meta = ScoreMetadata(
                    name='restraint_{0}'.format(metric),
                    title='Restraint Satisfaction',
                    unit=restraint_units[metric],
                    guide=1.0, lower=0.0, upper='95%', order=2,
            )
            record[meta.name] = np.max(values)
            metadata[meta.name] = meta

        for metric, values_by_residue in restraint_values_by_residue.items():
            if len(values_by_residue) <= 1:
                continue

            for i in values_by_residue:
                # I want to put the amino acid in these names, because I think 
                # it looks nice, but it causes problems for positions that can 
                # mutate.  So I assume that if a position has a sidechain 
                # restraint, it must not be allowed to mutate.
                aa = sequence_map[i] if is_sidechain_restraint[i] else 'X'
                res = '{0}{1}'.format(aa, i)
                meta = ScoreMetadata(
                        name='restraint_{0}_{1}'.format(metric, res.lower()),
                        title='Restraint Satisfaction for {0}'.format(res),
                        unit=restraint_units[metric],
                        guide=1.0, lower=0.0, upper='95%', order=3,
                )
                record[meta.name] = np.max(values_by_residue[i])
                metadata[meta.name] = meta

        # Finish calculating some records that depend on the whole structure.

        record['sequence'] = sequence
        for i, score in dunbrack_scores.items():
            aa = sequence_map[i] if is_sidechain_restraint[i] else 'X'
            res = '{0}{1}'.format(aa, i)
            meta = ScoreMetadata(
                    name='dunbrack_score_{0}'.format(res.lower()),
                    title='Dunbrack Score for {0}'.format(res),
                    unit='REU',
                    order=5,
            )
            record[meta.name] = score
            metadata[meta.name] = meta

        records.append(record)

    if pdb_paths:
        sys.stdout.write('\n')

    return records, metadata
945140
public ResidueRangeAndLength trimToValidResidues(ResidueRange rr) {
		ResidueNumber start = rr.getStart();
		ResidueNumber end = rr.getEnd();
		String chain = rr.getChainName();
		// Add chainName
		if(start.getChainName() == null) {
			start = new ResidueNumber(chain,start.getSeqNum(),start.getInsCode());
		}
		if(end.getChainName() == null) {
			end = new ResidueNumber(chain,end.getSeqNum(),end.getInsCode());
		}
		// Check that start and end are present in the map.
		// If not, try to find the next valid residue
		// (terminal residues sometimes lack CA atoms, so they don't appear)
		Integer startIndex = getPosition(start);
		if( startIndex == null) {
			// Assume that the residue numbers are sequential
			// Find startIndex such that the SeqNum is bigger than start's seqNum
			for(ResidueNumber key :  treeMap.keySet()) {
				if( !key.getChainName().equals(chain) )
					continue;
				if( start.getSeqNum() <= key.getSeqNum() ) {
					start = key;
					startIndex = getPosition(key);
					break;
				}
			}
			if( startIndex == null ) {
				logger.error("Unable to find Residue {} in AtomPositionMap, and no plausible substitute.",start);
				return null;
			} else {
				logger.warn("Unable to find Residue {}, so substituting {}.",rr.getStart(),start);
			}
		}
		Integer endIndex = getPosition(end);
		if( endIndex == null) {
			// Assume that the residue numbers are sequential
			// Find startIndex such that the SeqNum is bigger than start's seqNum
			for(ResidueNumber key :  treeMap.descendingKeySet()) {
				if( !key.getChainName().equals(chain) )
					continue;
				Integer value = getPosition(key);
				if( value < startIndex ) {
					// start is before the end!
					break;
				}
				if( end.getSeqNum() >= key.getSeqNum() ) {
					end = key;
					endIndex = value;
					break;
				}
			}
			if( endIndex == null ) {
				logger.error("Unable to find Residue {} in AtomPositionMap, and no plausible substitute.",end);
				return null;
			} else {
				logger.warn("Unable to find Residue {}, so substituting {}.",rr.getEnd(),end);
			}
		}

		// now use those to calculate the length
		// if start or end is null, will throw NPE
		int length = getLength(startIndex, endIndex,chain);

		return new ResidueRangeAndLength(chain, start, end, length);
	}
682772
def remove_isolated_residues(seq, len = 2)
      gap_starts = seq.to_enum(:scan, /(-\w{1,#{len}}-)/i).map { |_m| $`.size + 1 }
      # remove isolated residues
      gap_starts.each do |i|
        (i..i + len - 1).each { |j| seq[j] = '-' if isalpha(seq[j]) }
      end
      # remove isolated gaps
      res_starts = seq.to_enum(:scan, /([?\w]-{1,2}[?\w])/i).map { |_m| $`.size + 1 }
      res_starts.each do |i|
        (i..i + len - 1).each { |j| seq[j] = '?' if seq[j] == '-' }
      end
      seq
    end
620515
def _get_offset(self):
        
        sx, sxy, syx, sy, ox, oy = self.transformation
        return (ox, oy)
1064105
def coding_sequence(self):
        
        if self.sequence is None:
            return None

        start = self.first_start_codon_spliced_offset
        end = self.last_stop_codon_spliced_offset

        # If start codon is the at nucleotide offsets [3,4,5] and
        # stop codon is at nucleotide offsets  [20,21,22]
        # then start = 3 and end = 22.
        #
        # Adding 1 to end since Python uses non-inclusive ends in slices/ranges.

        # pylint: disable=invalid-slice-index
        # TODO(tavi) Figure out pylint is not happy with this slice
        return self.sequence[start:end + 1]
1038882
def check_residues(self):
        ''''''
        for chain, residue_ids in self.residues.iteritems():
            for residue_id, residue in residue_ids.iteritems():
                for record_type, atoms in residue.records.iteritems():
                    freq = {}
                    for atom in atoms:
                        rec_id = atom.name + atom.conformation
                        freq[rec_id] = freq.get(rec_id, 0)
                        freq[rec_id] += 1
                    for atom_type, count in freq.items():
                        if count > 1:
                            raise Exception('{0} occurrences of atom type {1} for record type {2} occur in residue {3} in chain {4}.'.format(count, atom_type, record_type, residue_id.strip(), chain))
1233213
def _reductionForOffset(self, offset):
        
        reduction = 0
        for (thisOffset, thisReduction) in self._adjustments:
            if offset >= thisOffset:
                reduction += thisReduction
            else:
                break
        return reduction
454444
def bytealign(self):
        
        skipped = (8 - (self._pos % 8)) % 8
        self.pos += self._offset + skipped
        assert self._assertsanity()
        return skipped
994924
def refine_hydrophobic(self, all_h, pistacks):
        """"""
        sel = {}
        #  1. Rings interacting via stacking can't have additional hydrophobic contacts between each other.
        for pistack, h in itertools.product(pistacks, all_h):
            h1, h2 = h.bsatom.idx, h.ligatom.idx
            brs, lrs = [p1.idx for p1 in pistack.proteinring.atoms], [p2.idx for p2 in pistack.ligandring.atoms]
            if h1 in brs and h2 in lrs:
                sel[(h1, h2)] = "EXCLUDE"
        hydroph = [h for h in all_h if not (h.bsatom.idx, h.ligatom.idx) in sel]
        sel2 = {}
        #  2. If a ligand atom interacts with several binding site atoms in the same residue,
        #  keep only the one with the closest distance
        for h in hydroph:
            if not (h.ligatom.idx, h.resnr) in sel2:
                sel2[(h.ligatom.idx, h.resnr)] = h
            else:
                if sel2[(h.ligatom.idx, h.resnr)].distance > h.distance:
                    sel2[(h.ligatom.idx, h.resnr)] = h
        hydroph = [h for h in sel2.values()]
        hydroph_final = []
        bsclust = {}
        #  3. If a protein atom interacts with several neighboring ligand atoms, just keep the one with the closest dist
        for h in hydroph:
            if h.bsatom.idx not in bsclust:
                bsclust[h.bsatom.idx] = [h, ]
            else:
                bsclust[h.bsatom.idx].append(h)

        idx_to_h = {}
        for bs in [a for a in bsclust if len(bsclust[a]) == 1]:
            hydroph_final.append(bsclust[bs][0])

        # A list of tuples with the idx of an atom and one of its neighbours is created
        for bs in [a for a in bsclust if not len(bsclust[a]) == 1]:
            tuples = []
            all_idx = [i.ligatom.idx for i in bsclust[bs]]
            for b in bsclust[bs]:
                idx = b.ligatom.idx
                neigh = [na for na in pybel.ob.OBAtomAtomIter(b.ligatom.OBAtom)]
                for n in neigh:
                    n_idx = n.GetIdx()
                    if n_idx in all_idx:
                        if n_idx < idx:
                            tuples.append((n_idx, idx))
                        else:
                            tuples.append((idx, n_idx))
                        idx_to_h[idx] = b

            tuples = list(set(tuples))
            tuples = sorted(tuples, key=itemgetter(1))
            clusters = cluster_doubles(tuples)  # Cluster connected atoms (i.e. find hydrophobic patches)

            for cluster in clusters:
                min_dist = float('inf')
                min_h = None
                for atm_idx in cluster:
                    h = idx_to_h[atm_idx]
                    if h.distance < min_dist:
                        min_dist = h.distance
                        min_h = h
                hydroph_final.append(min_h)
        before, reduced = len(all_h), len(hydroph_final)
        if not before == 0 and not before == reduced:
            write_message('Reduced number of hydrophobic contacts from %i to %i.\n' % (before, reduced), indent=True)
        return hydroph_final
647909
def pistacking(rings_bs, rings_lig):
    """"""
    data = namedtuple(
        'pistack', 'proteinring ligandring distance angle offset type restype resnr reschain restype_l resnr_l reschain_l')
    pairings = []
    for r, l in itertools.product(rings_bs, rings_lig):
        # DISTANCE AND RING ANGLE CALCULATION
        d = euclidean3d(r.center, l.center)
        b = vecangle(r.normal, l.normal)
        a = min(b, 180 - b if not 180 - b < 0 else b)  # Smallest of two angles, depending on direction of normal

        # RING CENTER OFFSET CALCULATION (project each ring center into the other ring)
        proj1 = projection(l.normal, l.center, r.center)
        proj2 = projection(r.normal, r.center, l.center)
        offset = min(euclidean3d(proj1, l.center), euclidean3d(proj2, r.center))

        # RECEPTOR DATA
        resnr, restype, reschain = whichresnumber(r.atoms[0]), whichrestype(r.atoms[0]), whichchain(r.atoms[0])
        resnr_l, restype_l, reschain_l = whichresnumber(l.orig_atoms[0]), whichrestype(
            l.orig_atoms[0]), whichchain(l.orig_atoms[0])

        # SELECTION BY DISTANCE, ANGLE AND OFFSET
        passed = False
        if not config.MIN_DIST < d < config.PISTACK_DIST_MAX:
            continue
        if 0 < a < config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
            ptype = 'P'
            passed = True
        if 90 - config.PISTACK_ANG_DEV < a < 90 + config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
            ptype = 'T'
            passed = True
        if passed:
            contact = data(proteinring=r, ligandring=l, distance=d, angle=a, offset=offset,
                           type=ptype, resnr=resnr, restype=restype, reschain=reschain,
                           resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l)
            pairings.append(contact)
    return filter_contacts(pairings)
1220400
def p_reference_variable_string_offset(p):
    ''
    p[0] = ast.StringOffset(p[1], p[3], lineno=p.lineno(2))
851633
def line_to_offset(self, line, column):
    
    line -= 1
    if line >= len(self._line_offsets):
      return self._text_len
    elif line < 0:
      return 0
    else:
      return min(self._line_offsets[line] + max(0, column), self._text_len)
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();
	}
224894
def split_pdb_residue(s):
    ''''''
    if s.isdigit():
        return (int(s), ' ')
    else:
        assert(s[:-1].isdigit())
        return ((s[:-1], s[-1]))
532580
private String getShiftDecimal(long number, int shift) {
        StringBuilder builder = new StringBuilder();
        builder.append(number);

        int index = builder.length() - shift;
        builder.insert(index, ".");

        return builder.toString();
    }