Run MSMS on this structure and store the residue depths/ca depths in the corresponding ChainProp SeqRecords
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))
735422
def _validate_residue_types(self):
        ''''''

        for chain_id, sequence_map in self.rosetta_to_atom_sequence_maps.iteritems():
            rosetta_sequence = self.rosetta_sequences[chain_id]
            atom_sequence = self.atom_sequences[chain_id]
            for rosetta_id, atom_id, _ in sequence_map:
                assert(rosetta_sequence[rosetta_id].ResidueAA == atom_sequence[atom_id].ResidueAA)

        for chain_id, sequence_map in self.atom_to_seqres_sequence_maps.iteritems():
            atom_sequence = self.atom_sequences[chain_id]
            seqres_sequence = self.seqres_sequences[chain_id]
            for atom_id, seqres_id, _ in sorted(sequence_map):
                assert(atom_sequence[atom_id].ResidueAA == seqres_sequence[seqres_id].ResidueAA)

        for chain_id, sequence_map in self.seqres_to_uniparc_sequence_maps.iteritems():
            if self.pdb_chain_to_uniparc_chain_mapping.get(chain_id):
                seqres_sequence = self.seqres_sequences[chain_id]
                uniparc_sequence = self.uniparc_sequences[self.pdb_chain_to_uniparc_chain_mapping[chain_id]]
                for seqres_id, uniparc_id_resid_pair, substitution_match in sequence_map:
                    uniparc_id = uniparc_id_resid_pair[1]
                    # Some of the matches may not be identical but all the '*' Clustal Omega matches should be identical
                    if substitution_match and substitution_match.clustal == 1:
                        assert(seqres_sequence[seqres_id].ResidueAA == uniparc_sequence[uniparc_id].ResidueAA)
61602
def sequence_dir(self):
        """"""
        if self.root_dir:
            return op.join(self.protein_dir, 'sequences')
        else:
            log.debug('Root directory not set')
            return None
1144467
def generate_top_psms(psms, protcol):
    """"""
    top_ms1_psms = {}
    for psm in psms:
        protacc = psm[protcol]
        precursor_amount = psm[mzidtsvdata.HEADER_PRECURSOR_QUANT]
        if ';' in protacc or precursor_amount == 'NA':
            continue
        precursor_amount = float(precursor_amount)
        psm_seq = psm[mzidtsvdata.HEADER_PEPTIDE]
        try:
            peptide_area = top_ms1_psms[protacc][psm_seq]
        except KeyError:
            try:
                top_ms1_psms[protacc][psm_seq] = precursor_amount
            except KeyError:
                top_ms1_psms[protacc] = {psm_seq: precursor_amount}
        else:
            if precursor_amount > peptide_area:
                top_ms1_psms[protacc][psm_seq] = precursor_amount
    return top_ms1_psms
567898
def _assemble_autophosphorylation(self, stmt):
        """"""
        sub_agent = deepcopy(stmt.enz)
        mc = stmt._get_mod_condition()
        sub_agent.mods.append(mc)
        # FIXME Ignore any bound conditions on the substrate!!!
        # This is because if they are included, a complex node will be returned,
        # which (at least currently) won't incorporate any protein
        # modifications.
        sub_agent.bound_conditions = []
        # FIXME
        self._add_nodes_edges(stmt.enz, sub_agent, pc.DIRECTLY_INCREASES,
                              stmt.evidence)
580998
def bam_to_fastq_pair(in_file, target_region, pair):
    
    space, start, end = target_region
    bam_file = pysam.Samfile(in_file, "rb")
    for read in bam_file:
        if (not read.is_unmapped and not read.mate_is_unmapped
                and bam_file.getrname(read.tid) == space
                and bam_file.getrname(read.mrnm) == space
                and read.pos >= start and read.pos <= end
                and read.mpos >= start and read.mpos <= end
                and not read.is_secondary
                and read.is_paired and getattr(read, "is_read%s" % pair)):
            seq = Seq.Seq(read.seq)
            qual = list(read.qual)
            if read.is_reverse:
                seq = seq.reverse_complement()
                qual.reverse()
            yield read.qname, str(seq), "".join(qual)
707707
def map_to_precursors_on_fly(seqs, names, loci, args):
    """"""
    precursor = precursor_sequence(loci, args.ref).upper()
    dat = dict()
    for s, n in itertools.izip(seqs, names):
        res = pyMatch.Match(precursor, str(s), 1, 3)
        if res > -1:
            dat[n] = [res, res + len(s)]
    logger.debug("mapped in %s: %s out of %s" % (loci, len(dat), len(seqs)))
    return dat
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
767466
def create_pymol_selection_from_PDB_residue_ids(residue_list):
    ''''''
    residues_by_chain = {}
    for residue_id in residue_list:
        chain_id = residue_id[0]
        pruned_residue_id = residue_id[1:]
        residues_by_chain[chain_id] = residues_by_chain.get(chain_id, [])
        residues_by_chain[chain_id].append(pruned_residue_id)

    str = []
    for chain_id, residue_list in sorted(residues_by_chain.iteritems()):
        str.append('(chain %s and resi %s)' % (chain_id, '+'.join(map(string.strip, sorted(residue_list)))))
    return ' or '.join(str)
250309
def _assemble_transphosphorylation(self, stmt):
        """"""
        # Check our assumptions about the bound condition of the enzyme
        assert len(stmt.enz.bound_conditions) == 1
        assert stmt.enz.bound_conditions[0].is_bound
        # Create a modified protein node for the bound target
        sub_agent = deepcopy(stmt.enz.bound_conditions[0].agent)
        sub_agent.mods.append(stmt._get_mod_condition())
        self._add_nodes_edges(stmt.enz, sub_agent, pc.DIRECTLY_INCREASES,
                              stmt.evidence)
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
200390
def from_sequence(chain, list_of_residues, sequence_type = None):
        ''''''
        s = Sequence(sequence_type)
        count = 1
        for ResidueAA in list_of_residues:
            s.add(Residue(chain, count, ResidueAA, sequence_type))
            count += 1
        return s
712932
def create_isobaric_quant_lookup(fns, tablefn_map, featmap, pqdb,
                                 featcolnr, allquantcols, psmcolmap):
    """"""
    pqdb.store_quant_channels(
            map_psmnrcol_to_quantcol(allquantcols, psmcolmap, tablefn_map),
            psmcolmap)
    quantmap = pqdb.get_quantchannel_map()
    to_store = []
    for fn, header, pquant in tsvreader.generate_tsv_pep_protein_quants(fns):
        pqdata = get_isob_quant_data(pquant, header, tablefn_map[fn],
                                     featmap, featcolnr, quantmap)
        to_store.extend(pqdata)
        if len(to_store) > 10000:
            pqdb.store_isobaric_quants(to_store, psmcolmap)
            to_store = []
    pqdb.store_isobaric_quants(to_store, psmcolmap)
39877
def sequence_prep(self):
        
        # Create a sorted list of all the FASTA files in the sequence path
        strains = sorted(glob(os.path.join(self.fastapath, '*.fa*'.format(self.fastapath))))
        for sample in strains:
            # Create the object
            metadata = MetadataObject()
            # Set the sample name to be the file name of the sequence by removing the path and file extension
            sample_name = os.path.splitext(os.path.basename(sample))[0]
            if sample_name in self.strainset:
                # Extract the OLNID from the dictionary using the SEQID
                samplename = self.straindict[sample_name]
                # samplename = sample_name
                # Set and create the output directory
                outputdir = os.path.join(self.path, samplename)
                make_path(outputdir)
                # Set the name of the JSON file
                json_metadata = os.path.join(outputdir, '{name}.json'.format(name=samplename))
                if not os.path.isfile(json_metadata):
                    # Create the name and output directory attributes
                    metadata.name = samplename
                    metadata.seqid = sample_name
                    metadata.outputdir = outputdir
                    metadata.jsonfile = json_metadata
                    # Set the name of the FASTA file to use in the analyses
                    metadata.bestassemblyfile = os.path.join(metadata.outputdir,
                                                             '{name}.fasta'.format(name=metadata.name))
                    # Symlink the original file to the output directory
                    relative_symlink(sample, outputdir, '{sn}.fasta'.format(sn=metadata.name))
                    # Associate the corresponding FASTQ files with the assembly
                    metadata.fastqfiles = sorted(glob(os.path.join(self.fastqpath,
                                                                   '{name}*.gz'.format(name=metadata.name))))
                    metadata.forward_fastq, metadata.reverse_fastq = metadata.fastqfiles
                    # Write the object to file
                    self.write_json(metadata)
                else:
                    metadata = self.read_json(json_metadata)
                # Add the metadata object to the list of objects
                self.metadata.append(metadata)
450820
def get_combo_subseq_within_2_5D(protein, props, within, filter_resnums=None):
        """"""
        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': []}

        all_resnums = []
        for prop in props:
            tmp_results = protein.get_generic_subseq_within_2_5D(prop_name=prop, within=within,
                                                         filter_resnums=filter_resnums)
            all_resnums.extend(tmp_results['subseq_resnums'])
        final_resnums = list(set(all_resnums))
        sub_id = '{}-{}->{}_within_{}_{}_extracted'.format(protein.representative_structure.id,
                                                           protein.representative_chain,
                                                           protein.representative_sequence.id,
                                                           within, props)
        new_sp = protein.representative_sequence.get_subsequence(resnums=final_resnums, new_id=sub_id,
                                                                 copy_letter_annotations=False)
        return {'subseq_len': len(final_resnums), 'subseq': new_sp, 'subseq_resnums': final_resnums}
523769
def create_merge_psm_map(peptides, ns):
    """"""
    psmmap = {}
    for peptide in peptides:
        seq = reader.get_peptide_seq(peptide, ns)
        psm_ids = reader.get_psm_ids_from_peptide(peptide, ns)
        for psm_id in psm_ids:
            try:
                psmmap[seq][psm_id.text] = 1
            except KeyError:
                psmmap[seq] = {psm_id.text: 2}
    for seq, psm_id_dict in psmmap.items():
        psmmap[seq] = [x for x in psm_id_dict]
    return psmmap
491025
def generate_master_proteins(psms, protcol):
    """"""
    master_proteins = {}
    if not protcol:
        protcol = mzidtsvdata.HEADER_MASTER_PROT
    for psm in psms:
        protacc = psm[protcol]
        if ';' in protacc:
            continue
        master_proteins[protacc] = 1
    if 'NA' in master_proteins:
        master_proteins.pop('NA')
    if '' in master_proteins:
        master_proteins.pop('')
    for protacc in master_proteins:
        yield {prottabledata.HEADER_PROTEIN: protacc}