Extract the sequences from the PDB file, perform the alignment,
        and load the coordinates of the CA of the common residues.
965669
def _load_data_alignment(self, chain1, chain2):
        
        parser = PDB.PDBParser(QUIET=True)
        ppb = PDB.PPBuilder()

        structure1 = parser.get_structure(chain1, self.pdb1)
        structure2 = parser.get_structure(chain2, self.pdb2)

        seq1 = str(ppb.build_peptides(structure1)[0].get_sequence())
        seq2 = str(ppb.build_peptides(structure2)[0].get_sequence())

        # Alignment parameters taken from PconsFold renumbering script.
        align = pairwise2.align.globalms(seq1, seq2, 2, -1, -0.5, -0.1)[0]
        indexes = set(i for i, (s1, s2) in enumerate(zip(align[0], align[1]))
                      if s1 != '-' and s2 != '-')
        coord1 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure1.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)
        coord2 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure2.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)

        self.coord1 = coord1
        self.coord2 = coord2
        self.N = len(seq1)
790806
def retrieve_data_from_rcsb(cls, ligand_code, pdb_id = None, silent = True, cached_dir = None):
        ''''''
        if not silent:
            colortext.printf("Retrieving data from RCSB")
        if cached_dir:
            assert(os.path.exists(cached_dir))

        ligand_info_path, ligand_info, pdb_ligand_info, pdb_ligand_info_path = None, None, None, None
        if cached_dir:
            ligand_info_path = os.path.join(cached_dir, '{0}.cif'.format(ligand_code))
            if os.path.exists(ligand_info_path):
                ligand_info = read_file(ligand_info_path)
        if not ligand_info:
            ligand_info = retrieve_ligand_cif(ligand_code)
            if cached_dir:
                write_file(ligand_info_path, ligand_info)

        # Parse .cif
        l = cls(ligand_code)
        l.parse_cif(ligand_info)
        l.pdb_id = pdb_id or l.pdb_id
        has_pdb_id = l.pdb_id and (len(l.pdb_id) == 4) and (l.pdb_id != '?')  # the last case is unnecessary and will be short-cut but I included it to show possible values

        # Parse PDB XML
        if has_pdb_id:
            if cached_dir:
                pdb_ligand_info_path = os.path.join(cached_dir, '{0}.pdb.ligandinfo'.format(l.pdb_id.lower()))
                if os.path.exists(pdb_ligand_info_path):
                    pdb_ligand_info = read_file(pdb_ligand_info_path)
                else:
                    pdb_ligand_info = retrieve_pdb_ligand_info(l.pdb_id)
                    write_file(pdb_ligand_info_path, pdb_ligand_info)
            else:
                pdb_ligand_info = retrieve_pdb_ligand_info(l.pdb_id)
        if pdb_ligand_info:
            l.parse_pdb_ligand_info(pdb_ligand_info)

        # Retrive the diagram image
        l.get_diagram()

        return l
43882
def _create_sequences(self):
        ''''''

        # Create the Rosetta sequences and the maps from the Rosetta sequences to the ATOM sequences
        try:
            self.pdb.construct_pdb_to_rosetta_residue_map(self.rosetta_scripts_path, rosetta_database_path = self.rosetta_database_path, cache_dir = self.cache_dir)
        except PDBMissingMainchainAtomsException:
            self.pdb_to_rosetta_residue_map_error = True

        # Get all the Sequences
        if self.pdb_id not in do_not_use_the_sequence_aligner:
            self.uniparc_sequences = self.PDB_UniParc_SA.uniparc_sequences
        else:
            self.uniparc_sequences = self.sifts.get_uniparc_sequences()

        self.fasta_sequences = self.FASTA.get_sequences(self.pdb_id)
        self.seqres_sequences = self.pdb.seqres_sequences
        self.atom_sequences = self.pdb.atom_sequences

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

        # Update the chain types for the UniParc sequences
        uniparc_pdb_chain_mapping = {}
        if self.pdb_id not in do_not_use_the_sequence_aligner:
            for pdb_chain_id, matches in self.PDB_UniParc_SA.clustal_matches.iteritems():
                if matches:
                    # we are not guaranteed to have a match e.g. the short chain J in 1A2C, chimeras, etc.
                    uniparc_chain_id = matches.keys()[0]
                    assert(len(matches) == 1)
                    uniparc_pdb_chain_mapping[uniparc_chain_id] = uniparc_pdb_chain_mapping.get(uniparc_chain_id, [])
                    uniparc_pdb_chain_mapping[uniparc_chain_id].append(pdb_chain_id)
        else:
            for pdb_chain_id, uniparc_chain_ids in self.sifts.get_pdb_chain_to_uniparc_id_map().iteritems():
                for uniparc_chain_id in uniparc_chain_ids:
                    uniparc_pdb_chain_mapping[uniparc_chain_id] = uniparc_pdb_chain_mapping.get(uniparc_chain_id, [])
                    uniparc_pdb_chain_mapping[uniparc_chain_id].append(pdb_chain_id)

        for uniparc_chain_id, pdb_chain_ids in uniparc_pdb_chain_mapping.iteritems():
            sequence_type = set([self.seqres_sequences[p].sequence_type for p in pdb_chain_ids])
            assert(len(sequence_type) == 1)
            sequence_type = sequence_type.pop()
            assert(self.uniparc_sequences[uniparc_chain_id].sequence_type == None)
            self.uniparc_sequences[uniparc_chain_id].set_type(sequence_type)
            for p in pdb_chain_ids:
                self.pdb_chain_to_uniparc_chain_mapping[p] = uniparc_chain_id

        # Update the chain types for the FASTA sequences
        for chain_id, sequence in self.seqres_sequences.iteritems():
            self.fasta_sequences[chain_id].set_type(sequence.sequence_type)
749177
def fetch_pdb(self, pdb_code):
        
        self.pdb_path, self.pdb_text = self._fetch_pdb(pdb_code)
        self._df = self._construct_df(pdb_lines=self.pdb_text.splitlines(True))
        return self
1110300
def _extract_residues(self):
                
        
        r_list = self.handle.getElementsByTagName("response")
        
        r_data = {}
        
        for r in r_list:
            data = self._parse_residue(r)
            res_id = (data['model'], data['chain'], data['number']) # (A, 1, 1), (A, 1, 2), ...
            if not r_data.has_key(res_id):
                r_data[res_id] = data
            else: # Some atoms get repeated at the end with TER Hydrogens/oxygens
                r_data[res_id]['atoms'] += data['atoms'] # Append Atoms
        
        for key in sorted(r_data.keys()):
            yield r_data[key]
1156988
def reporter(self):
        
        # Initialise variables
        header = ''
        row = ''
        databasedict = dict()
        # Load the database sequence type into a dictionary
        strainprofile = os.path.join(self.profilelocation, 'strainprofiles.txt')
        databaseprofile = DictReader(open(strainprofile))
        # Put the strain profile dictionary into a more easily searchable format
        for data in databaseprofile:
            databasedict[data['Strain']] = data['SequenceType']
        for sample in self.metadata.samples:
            closestmatches = list()
            if sample[self.analysistype].reportdir != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Populate the header with the appropriate data, including all the genes in the list of targets
                    header = 'Strain,SequenceType,Matches,Mismatches,NA,TotalGenes,ClosestDatabaseMatch,{},\n' \
                        .format(','.join(sorted(sample[self.analysistype].allelenames)))
                sortedmatches = sorted(sample[self.analysistype].profilematches.items(),
                                       key=operator.itemgetter(1), reverse=True)[0]
                closestseqtype = sortedmatches[0]
                # Pull out the closest database match
                for strain, seqtype in databasedict.items():
                    if seqtype == closestseqtype:
                        closestmatches.append(strain)
                sample[self.analysistype].closestseqtype = closestseqtype
                nummatches = int(sortedmatches[1])
                numna = 0
                queryallele = list()
                # Get all the alleles into a list
                for gene, allele in sorted(sample[self.analysistype].profiledata[closestseqtype].items()):
                    try:
                        # Extract the allele (anything after the -) from the allele matches
                        query = sample[self.analysistype].allelematches[gene].split('-')[1]
                        if allele == query:
                            queryallele.append(query)
                        else:
                            queryallele.append('{} ({})'.format(query, allele))
                    except KeyError:
                        queryallele.append('NA')
                        numna += 1
                mismatches = len(sample[self.analysistype].alleles) - nummatches - numna
                row += '{},{},{},{},{},{},{},{}'\
                    .format(sample.name, closestseqtype, nummatches, mismatches, numna,
                            len(sample[self.analysistype].alleles), ','.join(closestmatches), ','.join(queryallele))
                row += '\n'
        # Create the report folder
        make_path(self.reportpath)
        # Create the report containing all the data from all samples
        with open(os.path.join(self.reportpath, '{}.csv'.format(self.analysistype)), 'w') as combinedreport:
            # Write the results to this report
            combinedreport.write(header)
            combinedreport.write(row)
382755
def main(inputstructs, inputpdbids):
    """"""
    pdbid, pdbpath = None, None
    # #@todo For multiprocessing, implement better stacktracing for errors
    # Print title and version
    title = "* Protein-Ligand Interaction Profiler v%s *" % __version__
    write_message('\n' + '*' * len(title) + '\n')
    write_message(title)
    write_message('\n' + '*' * len(title) + '\n\n')
    outputprefix = config.OUTPUTFILENAME

    if inputstructs is not None:  # Process PDB file(s)
        num_structures = len(inputstructs)
        inputstructs = remove_duplicates(inputstructs)
        read_from_stdin = False
        for inputstruct in inputstructs:
            if inputstruct == '-':
                inputstruct = sys.stdin.read()
                read_from_stdin = True
                if config.RAWSTRING:
                    if sys.version_info < (3,):
                        inputstruct = bytes(inputstruct).decode('unicode_escape')
                    else:
                        inputstruct = bytes(inputstruct, 'utf8').decode('unicode_escape')
            else:
                if os.path.getsize(inputstruct) == 0:
                    sysexit(2, 'Empty PDB file\n')  # Exit if input file is empty
                if num_structures > 1:
                    basename = inputstruct.split('.')[-2].split('/')[-1]
                    config.OUTPATH = '/'.join([config.BASEPATH, basename])
                    outputprefix = 'report'
            process_pdb(inputstruct, config.OUTPATH, as_string=read_from_stdin, outputprefix=outputprefix)
    else:  # Try to fetch the current PDB structure(s) directly from the RCBS server
        num_pdbids = len(inputpdbids)
        inputpdbids = remove_duplicates(inputpdbids)
        for inputpdbid in inputpdbids:
            pdbpath, pdbid = download_structure(inputpdbid)
            if num_pdbids > 1:
                config.OUTPATH = '/'.join([config.BASEPATH, pdbid[1:3].upper(), pdbid.upper()])
                outputprefix = 'report'
            process_pdb(pdbpath, config.OUTPATH, outputprefix=outputprefix)

    if (pdbid is not None or inputstructs is not None) and config.BASEPATH is not None:
        if config.BASEPATH in ['.', './']:
            write_message('\nFinished analysis. Find the result files in the working directory.\n\n')
        else:
            write_message('\nFinished analysis. Find the result files in %s\n\n' % config.BASEPATH)
5742
def get_msms_annotations(self, outdir, force_rerun=False):
        
        # Now can run on Biopython Model objects exclusively thanks to Biopython updates
        # if self.file_type != 'pdb':
        #     raise ValueError('{}: unable to run MSMS with "{}" file type. Please change file type to "pdb"'.format(self.id,
        #                                                                                                     self.file_type))

        if self.structure:
            parsed = self.structure
        else:
            parsed = self.parse_structure()

        if not parsed:
            log.error('{}: unable to open structure to run MSMS'.format(self.id))
            return

        log.debug('{}: running MSMS'.format(self.id))
        # PDB ID is currently set to the structure file so the output name is the same with _msms.df appended to it
        msms_results = ssbio.protein.structure.properties.msms.get_msms_df(model=parsed.first_model, pdb_id=self.structure_path,
                                                                           outdir=outdir, force_rerun=force_rerun)
        if msms_results.empty:
            log.error('{}: unable to run MSMS'.format(self.id))
            return

        chains = msms_results.chain.unique()

        for chain in chains:
            res_depths = msms_results[msms_results.chain == chain].res_depth.tolist()
            ca_depths = msms_results[msms_results.chain == chain].ca_depth.tolist()

            chain_prop = self.chains.get_by_id(chain)
            chain_seq = chain_prop.seq_record

            # Making sure the X's are filled in
            res_depths = ssbio.protein.structure.properties.residues.match_structure_sequence(orig_seq=chain_seq,
                                                                                              new_seq=res_depths,
                                                                                              fill_with=float('Inf'))

            ca_depths = ssbio.protein.structure.properties.residues.match_structure_sequence(orig_seq=chain_seq,
                                                                                             new_seq=ca_depths,
                                                                                             fill_with=float('Inf'))

            chain_prop.seq_record.letter_annotations['RES_DEPTH-msms'] = res_depths
            chain_prop.seq_record.letter_annotations['CA_DEPTH-msms'] = ca_depths
            log.debug('{}: stored residue depths in chain seq_record letter_annotations'.format(chain))
661615
def get_sequences(self, pdb_id = None):
        ''''''
        sequences = {}
        if pdb_id:
            for chain_id, sequence in self.get(pdb_id, {}).iteritems():
                sequences[chain_id] = Sequence.from_sequence(chain_id, sequence)
        else:
            for pdb_id, v in self.iteritems():
                sequences[pdb_id] = {}
                for chain_id, sequence in v.iteritems():
                    sequences[pdb_id][chain_id] = Sequence.from_sequence(chain_id, sequence)
        return sequences
218593
def get_column(self, chrom, position,
                 missing_seqs=MissingSequenceHandler.TREAT_AS_ALL_GAPS,
                 species=None):
    """"""
    blocks = self.get_blocks(chrom, position, position + 1)
    if len(blocks) == 0:
      raise NoSuchAlignmentColumnError("Request for column on chrom " +
                                       chrom + " at position " +
                                       str(position) + " not possible; " +
                                       "genome alignment not defined at " +
                                       "that locus.")
    if len(blocks) > 1:
      raise NoUniqueColumnError("Request for column on chrom " + chrom +
                                " at position " + str(position) + "not " +
                                "possible; ambiguous alignment of that locus.")

    return blocks[0].get_column_absolute(position, missing_seqs, species)
599145
def get_sequences_from_cluster(c1, c2, data):
    """"""
    seqs1 = data[c1]['seqs']
    seqs2 = data[c2]['seqs']
    seqs = list(set(seqs1 + seqs2))
    names = []
    for s in seqs:
        if s in seqs1 and s in seqs2:
            names.append("both")
        elif s in seqs1:
            names.append(c1)
        else:
            names.append(c2)
    return seqs, names
1060006
def parse_pdb_file(self):
        """"""
        self.pdb_parse_tree = {'info': {},
                               'data': {
                                   self.state: {}}
                               }
        try:
            for line in self.pdb_lines:
                self.current_line = line
                record_name = line[:6].strip()
                if record_name in self.proc_functions:
                    self.proc_functions[record_name]()
                else:
                    if record_name not in self.pdb_parse_tree['info']:
                        self.pdb_parse_tree['info'][record_name] = []
                    self.pdb_parse_tree['info'][record_name].append(line)
        except EOFError:
            # Raised by END record
            pass
        if self.new_labels:
            ampal_data_session.commit()
        return
1141671
def df_pdb_metadata(self):
        """"""
        if self.num_structures == 0:
            log.error('No experimental PDB structures have been mapped to protein')
            return pd.DataFrame()

        pdb_pre_df = []

        not_showing_info = []
        for p in self.get_experimental_structures():
            if not p.structure_file:
                not_showing_info.append(p.id)
                log.debug('{}: PDB file has not been downloaded, run "pdb_downloader_and_metadata" \
                    to download all structures'.format(p.id))
                continue

            infodict = p.get_dict(df_format=True)
            infodict['pdb_id'] = p.id
            pdb_pre_df.append(infodict)

        log.warning('Not showing info for mapped PDB entries, as they have not been downloaded yet. Run pdb_downloader_and_metadata')

        cols = ['pdb_id', 'pdb_title', 'description', 'experimental_method', 'mapped_chains',
                'resolution', 'chemicals', 'date', 'taxonomy_name',
                'structure_file']
        df = pd.DataFrame.from_records(pdb_pre_df, columns=cols).set_index('pdb_id')
        return ssbio.utils.clean_df(df)
164922
def static_get_pdb_object(pdb_id, bio_cache = None, cache_dir = None):
        ''''''
        pdb_id = pdb_id.upper()

        if bio_cache:
            return bio_cache.get_pdb_object(pdb_id)

        if cache_dir:
            # Check to see whether we have a cached copy of the PDB file
            filepath = os.path.join(cache_dir, '{0}.pdb'.format(pdb_id))
            if os.path.exists(filepath):
                return PDB.from_filepath(filepath)

        # Get any missing files from the RCSB and create cached copies if appropriate
        pdb_contents = retrieve_pdb(pdb_id)
        if cache_dir:
            write_file(os.path.join(cache_dir, "%s.pdb" % pdb_id), pdb_contents)
        return PDB(pdb_contents)
413847
def parse_atoms(self):
        ''''''

        atom_site_header_tag = self.main_tag.getElementsByTagName("PDBx:atom_siteCategory")
        assert(len(atom_site_header_tag) == 1)
        atom_site_header_tag = atom_site_header_tag[0]

        atom_site_tags = atom_site_header_tag.getElementsByTagName("PDBx:atom_site")

        residue_map = {}
        residues_read = {}
        int_type = types.IntType
        for t in atom_site_tags:
            r, seqres, ResidueAA, Residue3AA = PDBML_slow.parse_atom_site(t, self.modified_residues)
            if r:
                # skip certain ACE residues
                if not(self.pdb_id in cases_with_ACE_residues_we_can_ignore and Residue3AA == 'ACE'):
                    full_residue_id = str(r)
                    if residues_read.get(full_residue_id):
                        assert(residues_read[full_residue_id] == (r.ResidueAA, seqres))
                    else:
                        residues_read[full_residue_id] = (r.ResidueAA, seqres)
                        residue_map[r.Chain] = residue_map.get(r.Chain, {})
                        assert(type(seqres) == int_type)
                        residue_map[r.Chain][str(r)] = seqres

        ## Create SequenceMap objects to map the ATOM Sequences to the SEQRES Sequences
        atom_to_seqres_sequence_maps = {}
        for chain_id, atom_seqres_mapping in residue_map.iteritems():
            atom_to_seqres_sequence_maps[chain_id] = SequenceMap.from_dict(atom_seqres_mapping)

        self.atom_to_seqres_sequence_maps = atom_to_seqres_sequence_maps
1371184
def _guess_sequence_type_from_string(self, seq):
        ''''''
        # Define expected residues for each sequence type
        aa_chars = ['P','V','L','I','M','F','Y','W','H','K','R','Q','N','E','D','S','X','*']
        aas = set(itertools.chain(aa_chars, [lower(a) for a in aa_chars]))
        na_chars = ['A', 'T', 'G', 'C', 'N', 'U']
        nas = set(itertools.chain(na_chars, [lower(a) for a in na_chars]))

        num_nucleotide = 0
        num_protein = 0
        count = 0
        for residue in seq:
            if residue in nas:
                num_nucleotide += 1
            elif residue in aas:
                num_protein += 1
            else:
                raise Exception('Encountered unexpected character when attempting to guess sequence type: %s' % (residue))
            count += 1
            if count >300: break
        if float(num_protein) / (num_protein+num_nucleotide) > 0.1:
            return self.PROTEIN_SEQUENCE_TYPE
        else:
            return self.NUCLEOTIDE_SEQUENCE_TYPE
146192
def sequences(self):
        """"""
        lines = Lines()
        lines.add(0, '@cython.final')
        lines.add(0, 'cdef class Sequences(object):')
        for subseqs in self.model.sequences:
            lines.add(1, 'cdef public %s %s'
                         % (objecttools.classname(subseqs), subseqs.name))
        if getattr(self.model.sequences, 'states', None) is not None:
            lines.add(1, 'cdef public StateSequences old_states')
            lines.add(1, 'cdef public StateSequences new_states')
        for subseqs in self.model.sequences:
            print('        - %s' % subseqs.name)
            lines.add(0, '@cython.final')
            lines.add(0, 'cdef class %s(object):'
                         % objecttools.classname(subseqs))
            for seq in subseqs:
                ctype = 'double' + NDIM2STR[seq.NDIM]
                if isinstance(subseqs, sequencetools.LinkSequences):
                    if seq.NDIM == 0:
                        lines.add(1, 'cdef double *%s' % seq.name)
                    elif seq.NDIM == 1:
                        lines.add(1, 'cdef double **%s' % seq.name)
                        lines.add(1, 'cdef public int len_%s' % seq.name)
                else:
                    lines.add(1, 'cdef public %s %s' % (ctype, seq.name))
                lines.add(1, 'cdef public int _%s_ndim' % seq.name)
                lines.add(1, 'cdef public int _%s_length' % seq.name)
                for idx in range(seq.NDIM):
                    lines.add(1, 'cdef public int _%s_length_%d'
                                 % (seq.name, idx))
                if seq.NUMERIC:
                    ctype_numeric = 'double' + NDIM2STR[seq.NDIM+1]
                    lines.add(1, 'cdef public %s _%s_points'
                                 % (ctype_numeric, seq.name))
                    lines.add(1, 'cdef public %s _%s_results'
                                 % (ctype_numeric, seq.name))
                    if isinstance(subseqs, sequencetools.FluxSequences):
                        lines.add(1, 'cdef public %s _%s_integrals'
                                     % (ctype_numeric, seq.name))
                        lines.add(1, 'cdef public %s _%s_sum'
                                     % (ctype, seq.name))
                if isinstance(subseqs, sequencetools.IOSequences):
                    lines.extend(self.iosequence(seq))
            if isinstance(subseqs, sequencetools.InputSequences):
                lines.extend(self.load_data(subseqs))
            if isinstance(subseqs, sequencetools.IOSequences):
                lines.extend(self.open_files(subseqs))
                lines.extend(self.close_files(subseqs))
                if not isinstance(subseqs, sequencetools.InputSequence):
                    lines.extend(self.save_data(subseqs))
            if isinstance(subseqs, sequencetools.LinkSequences):
                lines.extend(self.set_pointer(subseqs))
        return lines