Get the sequence alignment information for a sequence to a structure's chain.
961759
def get_seqprop_to_seqprop_alignment_stats(self, seqprop1, seqprop2):
        """"""
        alignment = self._get_seqprop_to_seqprop_alignment(seqprop1=seqprop1, seqprop2=seqprop2)
        return ssbio.protein.sequence.utils.alignment.pairwise_alignment_stats(reference_seq_aln=alignment[0],
                                                                               other_seq_aln=alignment[1])
1110501
def align_two_alignments(aln1, aln2, params=None):
    
    if not params:
        params = {}

    #create SequenceCollection object from aln1
    aln1_collection = SequenceCollection(aln1)
    #Create mapping between abbreviated IDs and full IDs
    aln1_int_map, aln1_int_keys = aln1_collection.getIntMap(prefix='aln1_')
    #Create SequenceCollection from int_map.
    aln1_int_map = SequenceCollection(aln1_int_map)

    #create SequenceCollection object from aln2
    aln2_collection = SequenceCollection(aln2)
    #Create mapping between abbreviated IDs and full IDs
    aln2_int_map, aln2_int_keys = aln2_collection.getIntMap(prefix='aln2_')
    #Create SequenceCollection from int_map.
    aln2_int_map = SequenceCollection(aln2_int_map)

    #set output and profile options
    params.update({'-out':get_tmp_filename(), '-profile':True})

    #save aln1 to tmp file
    aln1_filename = get_tmp_filename()
    aln1_out = open(aln1_filename,'w')
    aln1_out.write(aln1_int_map.toFasta())
    aln1_out.close()

    #save aln2 to tmp file
    aln2_filename = get_tmp_filename()
    aln2_out = open(aln2_filename, 'w')
    aln2_out.write(aln2_int_map.toFasta())
    aln2_out.close()

    #Create Muscle app and get results
    app = Muscle(InputHandler='_input_as_multifile', params=params,
                 WorkingDir=tempfile.gettempdir())
    res = app((aln1_filename, aln2_filename))

    #Get alignment as dict out of results
    alignment = dict(parse_fasta(res['MuscleOut']))

    #Make new dict mapping original IDs
    new_alignment = {}
    for k,v in alignment.items():
        if k in aln1_int_keys:
            new_alignment[aln1_int_keys[k]] = v
        else:
            new_alignment[aln2_int_keys[k]] = v

    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment)

    #Clean up
    res.cleanUp()
    del(aln1_collection, aln1_int_map, aln1_int_keys)
    del(aln2_collection, aln2_int_map, aln2_int_keys)
    del(app, res, alignment, params)
    remove(aln1_filename)
    remove(aln2_filename)

    return new_alignment
1058663
def encode_data(self):
        """"""
        output_data = {}
        output_data["groupTypeList"] = encode_array(self.group_type_list, 4, 0)
        output_data["xCoordList"] = encode_array(self.x_coord_list, 10, 1000)
        output_data["yCoordList"] = encode_array(self.y_coord_list, 10, 1000)
        output_data["zCoordList"] = encode_array(self.z_coord_list, 10, 1000)
        output_data["bFactorList"] = encode_array(self.b_factor_list, 10, 100)
        output_data["occupancyList"] = encode_array(self.occupancy_list, 9, 100)
        output_data["atomIdList"] = encode_array(self.atom_id_list, 8, 0)
        output_data["altLocList"] = encode_array(self.alt_loc_list, 6, 0)
        output_data["insCodeList"] = encode_array(self.ins_code_list, 6, 0)
        output_data["groupIdList"] = encode_array(self.group_id_list, 8, 0)
        output_data["groupList"] = self.group_list
        output_data["sequenceIndexList"] = encode_array(self.sequence_index_list, 8, 0)
        output_data["chainNameList"] = encode_array(self.chain_name_list, 5, 4)
        output_data["chainIdList"] = encode_array(self.chain_id_list, 5, 4)
        output_data["bondAtomList"] = encode_array(self.bond_atom_list, 4, 0)
        output_data["bondOrderList"] = encode_array(self.bond_order_list, 2, 0)
        output_data["secStructList"] = encode_array(self.sec_struct_list, 2, 0)
        output_data["chainsPerModel"] = self.chains_per_model
        output_data["groupsPerChain"] = self.groups_per_chain
        output_data["spaceGroup"] = self.space_group
        output_data["mmtfVersion"] = self.mmtf_version
        output_data["mmtfProducer"] = self.mmtf_producer
        output_data["structureId"] = self.structure_id
        output_data["entityList"] = self.entity_list
        output_data["bioAssemblyList"] = self.bio_assembly
        output_data["rFree"] = self.r_free
        output_data["rWork"] = self.r_work
        output_data["resolution"] = self.resolution
        output_data["title"] = self.title
        output_data["experimentalMethods"] = self.experimental_methods
        output_data["depositionDate"] = self.deposition_date
        output_data["releaseDate"] = self.release_date
        output_data["unitCell"] = self.unit_cell
        output_data["numBonds"] = self.num_bonds
        output_data["numChains"] = self.num_chains
        output_data["numModels"] = self.num_models
        output_data["numAtoms"] = self.num_atoms
        output_data["numGroups"] = self.num_groups
        return output_data
176247
def load_sequences_to_reference(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _load_sequences_to_reference_gene_sc(g_id, outdir=self.sequences_by_gene_dir,
                                                 g_to_pickle=self.gene_protein_pickles,
                                                 strain_infodict=self.strain_infodict,
                                                 orth_matrix=self.df_orthology_matrix, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            from Bio import SeqIO
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
                protein_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_pickle_path)

                for strain, info in strain_infodict.items():
                    strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                    strain_gene_functional = info['functional_genes'][g_id]
                    if strain_gene_functional:
                        # Pull the gene ID of the strain from the orthology matrix
                        strain_gene_key = orth_matrix.at[g_id, strain]
                        new_id = '{}_{}'.format(g_id, strain)
                        if protein_pickle.sequences.has_id(new_id):
                            continue
                        protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                            ident=new_id,
                                                            set_as_representative=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_load_sequences_to_reference_gene_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        updated = []
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
            updated.append(g_id)
        not_updated = set(list(self.gene_protein_pickles.keys())).difference(updated)
        log.info('No change to {} genes, removing from gene_protein_pickles'.format(len(not_updated)))
        log.debug(not_updated)
        for rem in not_updated:
            del self.gene_protein_pickles[rem]
8632
def _load_sequences_to_strain(self, strain_id, force_rerun=False):
        """"""
        gp_seqs_path = op.join(self.model_dir, '{}_gp_withseqs.pckl'.format(strain_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=gp_seqs_path):
            gp_noseqs = ssbio.io.load_pickle(self.strain_infodict[strain_id]['gp_noseqs_path'])
            strain_sequences = SeqIO.index(self.strain_infodict[strain_id]['genome_path'], 'fasta')
            for strain_gene in gp_noseqs.functional_genes:
                # Pull the gene ID of the strain from the orthology matrix
                strain_gene_key = self.df_orthology_matrix.at[strain_gene.id, strain_id]
                # Load into the strain GEM-PRO
                new_id = '{}_{}'.format(strain_gene.id, strain_id)
                if strain_gene.protein.sequences.has_id(new_id):
                    continue
                strain_gene.protein.load_manual_sequence(seq=strain_sequences[strain_gene_key], ident=new_id,
                                                         set_as_representative=True)
            gp_noseqs.save_pickle(outfile=gp_seqs_path)

        return strain_id, gp_seqs_path
187839
private static ProteinSequence getProteinSequence(String str) {
		try {
			ProteinSequence s = new ProteinSequence(str);
			return s;
		} catch (CompoundNotFoundException e) {

			logger.error("Unexpected error when creating ProteinSequence",e);
		}
		return null;
	}
819060
def get_chain_sequence_string(self, chain_id, use_seqres_sequences_if_possible, raise_Exception_if_not_found = True):
        ''''''
        chain_pair = self.get_annotated_chain_sequence_string(chain_id, use_seqres_sequences_if_possible, raise_Exception_if_not_found = raise_Exception_if_not_found)
        if chain_pair:
            return chain_pair[1]
        return None
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)
1064170
def load_bioassembly_info_from_file(self, biomol_num):
        

        # current functionality is to take in a pre-assembled bioassembly file, parse it's MODELs and get info from that
        pass
        from Bio.PDB import PDBParser
        p = PDBParser(PERMISSIVE=True, QUIET=True)

        structure = StructureIO(self.structure_path, self.file_type)

        structure = p.get_structure('tmp', infile)
        bioass_to_chain_stoich = defaultdict(int)
        for model in structure:
            for chain in model:
                bioass_to_chain_stoich[chain.id] += 1
        return dict(bioass_to_chain_stoich)
1179580
def protein_statistics(self):
        """"""

        # TODO: can i use get_dict here instead
        d = {}

        d['id'] = self.id
        d['sequences'] = [x.id for x in self.sequences]
        d['num_sequences'] = self.num_sequences
        if self.representative_sequence:
            d['representative_sequence'] = self.representative_sequence.id
            d['repseq_gene_name'] = self.representative_sequence.gene_name
            d['repseq_uniprot'] = self.representative_sequence.uniprot
            d['repseq_description'] = self.representative_sequence.description
        d['num_structures'] = self.num_structures
        d['experimental_structures'] = [x.id for x in self.get_experimental_structures()]
        d['num_experimental_structures'] = self.num_structures_experimental
        d['homology_models'] = [x.id for x in self.get_homology_models()]
        d['num_homology_models'] = self.num_structures_homology
        if self.representative_structure:
            d['representative_structure'] = self.representative_structure.id
            d['representative_chain'] = self.representative_chain
            d['representative_chain_seq_coverage'] = self.representative_chain_seq_coverage
            d['repstruct_description'] = self.description
            if self.representative_structure.is_experimental:
                d['repstruct_resolution'] = self.representative_structure.resolution
        d['num_sequence_alignments'] = len(self.sequence_alignments)
        d['num_structure_alignments'] = len(self.structure_alignments)

        return d
113346
def seqrecord2sequence(record, qlen, index):
    
    # aligned_len = sum(map(str.isupper, str(record.seq)))
    # assert qlen == aligned_len, (
    #         "Aligned sequence length %s, query %s\n%s"
    #         % (aligned_len, qlen, str(record)))

    description = (record.description.split(' ', 1)[1]
            if ' ' in record.description
            else ' ')
    return dict(index=index,
                id=record.id,
                description=description,
                dbxrefs={},
                phylum='',
                taxchar='',
                head_len=None,
                tail_len=None,
                head_seq='',
                tail_seq='',
                length=len(record.seq) - record.seq.count('-'),
                seq=str(record.seq),
                )
9881
def get_surface_subseq_3D(protein,
                              depth_prop='RES_DEPTH-msms', depth_cutoff=2.5, depth_condition='<',
                              acc_prop='RSA_ALL-freesasa_het', acc_cutoff=25, acc_condition='>'):
        """"""
        empty = {'surface_3D': {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []},
                    'notdeep_3D': {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []},
                    'acc_3D'    : {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []}}
        if not protein.representative_structure:
            log.error('{}: no representative structure, cannot search for subseq'.format(protein.id))
            return empty

        notdeep_subseq, notdeep_subseq_resnums = protein.get_seqprop_subsequence_from_structchain_property(
            property_key=depth_prop,
            property_value=depth_cutoff,
            condition=depth_condition,
            use_representatives=True,
            return_resnums=True) or (None, [])

        acc_subseq, acc_subseq_resnums = protein.get_seqprop_subsequence_from_structchain_property(
            property_key=acc_prop,
            property_value=acc_cutoff,
            condition=acc_condition,
            use_representatives=True,
            return_resnums=True) or (None, [])

        surface_subseq_resnums = list(set(notdeep_subseq_resnums).intersection(acc_subseq_resnums))
        surface_subseq = protein.representative_sequence.get_subsequence(surface_subseq_resnums)

        all_info = {'surface_3D': {'subseq_len'    : len(surface_subseq_resnums), 'subseq': surface_subseq,
                                   'subseq_resnums': surface_subseq_resnums},
                    'notdeep_3D': {'subseq_len'    : len(notdeep_subseq_resnums), 'subseq': notdeep_subseq,
                                   'subseq_resnums': notdeep_subseq_resnums},
                    'acc_3D'    : {'subseq_len'    : len(acc_subseq_resnums), 'subseq': acc_subseq,
                                   'subseq_resnums': acc_subseq_resnums}}
        return all_info
1011226
def _get_ATOM_sequences_2(self):
        ''''''

        # Get a list of all residues with ATOM or HETATM records
        atom_sequences = {}
        structural_residue_IDs_set = set() # use a set for a quicker lookup
        ignore_HETATMs = True # todo: fix this if we need to deal with HETATMs
        for l in self.structure_lines:
            residue_type = l[17:20].strip()
            if l.startswith("HETATM"):
                if self.modified_residue_mapping_3.get(residue_type):
                    residue_type = self.modified_residue_mapping_3[residue_type]
                elif ignore_HETATMs:
                    continue

            residue_id = l[21:27]
            if residue_id not in structural_residue_IDs_set:
                chain_id = l[21]
                chain_type = self.chain_types[chain_id]
                atom_sequences[chain_id] = atom_sequences.get(chain_id, Sequence(chain_type))
                residue_type = l[17:20].strip()

                residue_type = self.modified_residue_mapping_3.get(residue_type, residue_type)
                short_residue_type = None
                if residue_type == 'UNK':
                    short_residue_type = 'X'
                elif chain_type == 'Unknown':
                    assert(False) # we should not reach here - Unknown chains should only contain UNK records
                elif chain_type == 'Protein' or chain_type == 'Protein skeleton':
                    short_residue_type = residue_type_3to1_map.get(residue_type) or protonated_residue_type_3to1_map.get(residue_type)
                elif chain_type == 'DNA':
                    short_residue_type = dna_nucleotides_2to1_map.get(residue_type) or non_canonical_dna.get(residue_type)
                elif chain_type == 'RNA':
                    short_residue_type = non_canonical_rna.get(residue_type) or residue_type
                elif not self.strict:
                    short_residue_type = 'X'
                else:
                    raise NonCanonicalResidueException("Unrecognized residue type %s in PDB file '%s', residue ID '%s'." % (residue_type, str(self.pdb_id), str(residue_id)))

                #structural_residue_IDs.append((residue_id, short_residue_type))
                atom_sequences[chain_id].add(PDBResidue(residue_id[0], residue_id[1:], short_residue_type, chain_type))
                structural_residue_IDs_set.add(residue_id)

        self.atom_sequences = atom_sequences
498594
def _load_sequences_to_reference_gene(self, g_id, force_rerun=False):
        """"""
        protein_seqs_pickle_path = op.join(self.sequences_by_gene_dir, '{}_protein_withseqs.pckl'.format(g_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
            protein_pickle_path = self.gene_protein_pickles[g_id]
            protein_pickle = ssbio.io.load_pickle(protein_pickle_path)

            for strain, info in self.strain_infodict.items():
                strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                strain_gene_functional = info['functional_genes'][g_id]
                if strain_gene_functional:
                    # Pull the gene ID of the strain from the orthology matrix
                    strain_gene_key = self.df_orthology_matrix.at[g_id, strain]
                    new_id = '{}_{}'.format(g_id, strain)
                    if protein_pickle.sequences.has_id(new_id):
                        continue
                    protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                        ident=new_id,
                                                        set_as_representative=False)
            protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

        return g_id, protein_seqs_pickle_path
370664
def convert_block_dicts_to_string(self, block_1st2nd, block_1st, block_2nd, block_3rd):
        """"""
        out = ""
        # We need 1st and 2nd positions
        if self.codon_positions in ['ALL', '1st-2nd']:
            for gene_code, seqs in block_1st2nd.items():
                out += '>{0}_1st-2nd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        elif self.codon_positions == '1st':
            for gene_code, seqs in block_1st.items():
                out += '>{0}_1st\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        elif self.codon_positions == '2nd':
            for gene_code, seqs in block_2nd.items():
                out += '>{0}_2nd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq

        # We also need 3rd positions
        if self.codon_positions in ['ALL', '3rd']:
            for gene_code, seqs in block_3rd.items():
                out += '\n>{0}_3rd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        return out
395665
def align(self, alignment_tool = 'clustalw', gap_opening_penalty = 0.2, ignore_bad_chains = False):
        ''''''
        if len(self.pdb_chains) > 1:
            sa = SequenceAligner(alignment_tool = alignment_tool, gap_opening_penalty = gap_opening_penalty)
            for pdb_chain in self.pdb_chains:
                sa.add_sequence('%s_%s' % (pdb_chain['pdb_id'], pdb_chain['chain_id']), pdb_chain['sequence'], ignore_bad_chains = ignore_bad_chains)
            best_matches = sa.align()
            return sa.alignment_output, best_matches
        else:
            raise Exception('Cannot align sequences - less than two chains were specified.')
455172
def from_block(cls, block):
        """"""
        rseqs = cma.realign_seqs(block)
        records = (SeqRecord(Seq(rseq, extended_protein),
                    id=bseq['id'],
                    description=bseq['description'],
                    dbxrefs=bseq['dbxrefs'].values(),   # list of strings
                    annotations=dict(
                        index=bseq['index'],
                        length=bseq['length'],
                        dbxrefs=bseq['dbxrefs'],
                        phylum=bseq['phylum'],
                        taxchar=bseq['taxchar'],
                        head_seq=bseq['head_seq'],
                        tail_seq=bseq['tail_seq'],
                        head_len=bseq['head_len'],
                        tail_len=bseq['tail_len'],      # dict
                        ),
                    # ENH: annotate with conservation levels
                    # letter_annotations=bseq['x'],
                    )
                for bseq, rseq in zip(block['sequences'], rseqs))
        return cls(records,
                # CMA attributes
                # block['one'],
                block['level'],
                block['name'],
                block['params'],
                block['query_length'],
                block['query_chars'],
                )