assigns secondary structure information based on a simple and very fast algorithm published by Zhang and Skolnick in their TM-align paper. Reference:  TM-align: a protein structure alignment algorithm based on the Tm-sutils (2005) NAR, 33(7) 2302-2309
967779
function assignHelixSheet(structure) {
  console.time('mol.assignHelixSheet');
  var chains = structure.chains();
  for (var ci = 0; ci < chains.length; ++ci) {
    var chain = chains[ci];
    chain.eachBackboneTrace(traceAssignHelixSheet);
  }
  console.timeEnd('mol.assignHelixSheet');
}
202958
def align_to_sort_bam(fastq1, fastq2, aligner, data):
    
    names = data["rgnames"]
    align_dir_parts = [data["dirs"]["work"], "align", names["sample"]]
    if data.get("disambiguate"):
        align_dir_parts.append(data["disambiguate"]["genome_build"])
    aligner_index = _get_aligner_index(aligner, data)
    align_dir = utils.safe_makedir(os.path.join(*align_dir_parts))
    ref_file = tz.get_in(("reference", "fasta", "base"), data)
    if fastq1.endswith(".bam"):
        data = _align_from_bam(fastq1, aligner, aligner_index, ref_file,
                               names, align_dir, data)
    else:
        data = _align_from_fastq(fastq1, fastq2, aligner, aligner_index, ref_file,
                                 names, align_dir, data)
    if data["work_bam"] and utils.file_exists(data["work_bam"]):
        if data.get("align_split") and dd.get_mark_duplicates(data):
            # If merging later with with bamsormadup need query sorted inputs
            # but CWL requires a bai file. Create a fake one to make it happy.
            bam.fake_index(data["work_bam"], data)
        else:
            bam.index(data["work_bam"], data["config"])
        for extra in ["-sr", "-disc"]:
            extra_bam = utils.append_stem(data['work_bam'], extra)
            if utils.file_exists(extra_bam):
                bam.index(extra_bam, data["config"])
    return data
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.')
512801
def concat_align(fastas):
    
    # read in sequences
    fa2len = {}
    seqs = {}
    IDs = []
    for fasta in fastas:
        seqs[fasta] = {}
        for seq in parse_fasta(fasta):
            ID = seq[0].split('>')[1].split()[0]
            IDs.append(ID)
            seqs[fasta][ID] = seq[1]
        fa2len[fasta] = len(seq[1])
    # concat sequences
    IDs = set(IDs)
    concat = {}
    for fasta in fastas:
        for ID in IDs:
            if ID not in concat:
                concat[ID] = []
            if ID not in seqs[fasta]:
                concat[ID].append('-'*fa2len[fasta])
            else:
                concat[ID].append(seqs[fasta][ID])
    return concat
1149178
def build_bam_tags():
    ''''''
    #pylint: disable=unused-argument
    def _combine_filters(fam, paired_align, align):
        filters = [x.filter_value for x in [fam, align] if x and x.filter_value]
        if filters:
            return ";".join(filters).replace('; ', ';')
        return None

    boolean_tag_value = {True:1}
    tags = [
        BamTag("X0", "Z",
               ("filter (why the alignment was excluded)"),
               _combine_filters),
        BamTag("X1", "Z",
               ("leftmost~rightmost matched pair positions"),
               lambda fam, pair, align: pair.positions('{left}~{right}')),
        BamTag("X2", "Z",
               ("L~R CIGARs"),
               lambda fam, pair, align: pair.cigars('{left}~{right}')),
        BamTag("X3", "i",
               "unique identifier for this alignment family",
               lambda fam, pair, align: fam.umi_sequence),
        BamTag("X4", "Z",
               ("L~R UMT barcodes for this alignment family; because "
                "of fuzzy matching the family UMT may be distinct "
                "from the UMT of the original alignment"),
               lambda fam, pair, align: fam.umt('{left}~{right}')),
        BamTag("X5", "i",
               "family size (number of align pairs in this family)",
               lambda fam, pair, align: fam.included_pair_count),
        BamTag("X6", "i",
               ("presence of this tag signals that this alignment "
                "would be the template for the consensus alignment"),
               lambda fam, pair, align: boolean_tag_value.get(fam.is_consensus_template(align), None))]
    return tags
176821
protected void align() {
		if (!isReady()) {
			return;
		}

		long timeStart = System.nanoTime();

		int[] dim = getScoreMatrixDimensions();
		if (storingScoreMatrix) {
			scores = new int[dim[0]][dim[1]][dim[2]];
		} else {
			scores = new int[dim[0]][][];
			scores[0] = new int[dim[1]][dim[2]];
			scores[1] = new int[dim[1]][dim[2]];
		}
		boolean linear = (gapPenalty.getType() == GapPenalty.Type.LINEAR);
		Last[][][] traceback = new Last[dim[0]][][];
		List<Step> sx = new ArrayList<Step>(), sy = new ArrayList<Step>();

		if (!local) {
			xyMax = new int[] { dim[0] - 1, dim[1] - 1 };
			xyStart = new int[] { 0, 0 };
			score = 0;
			List<Subproblem> problems = Subproblem.getSubproblems(anchors, xyMax[0], xyMax[1]);
			assert problems.size() == anchors.size() + 1;
			for (int i = 0; i < problems.size(); i++) {
				Subproblem subproblem = problems.get(i);
				for (int x = subproblem.getQueryStartIndex(); x <= subproblem.getQueryEndIndex(); x++) {

					traceback[x] =

					  linear ?
						setScoreVector(x, subproblem, gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x, subproblem), storingScoreMatrix, scores) :

						setScoreVector(x, subproblem, gapPenalty.getOpenPenalty(), gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x, subproblem), storingScoreMatrix, scores);
				}
			}
			setSteps(traceback, scores, sx, sy);
			score = Integer.MIN_VALUE;
			int[] finalScore = scores[xyMax[0]][xyMax[1]];
			for (int z = 0; z < finalScore.length; z++) {
				score = Math.max(score, finalScore[z]);
			}
		} else {
			for (int x = 0; x < dim[0]; x++) {

				traceback[x] =

				  linear ?
					setScoreVector(x, gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x), storingScoreMatrix, scores, xyMax, score) :

					setScoreVector(x, gapPenalty.getOpenPenalty(), gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x), storingScoreMatrix, scores, xyMax, score);

				if (xyMax[0] == x) {
					score = scores[x][xyMax[1]][0];
				}
			}
			xyStart = local ? setSteps(traceback, xyMax, sx, sy) : setSteps(traceback, scores, sx, sy);
		}

		setProfile(sx, sy);

		if (!storingScoreMatrix) {
			scores = null;
		}

		time = System.nanoTime() - timeStart;
	}
1348302
public void identify(final List<Chain> chains,
			final Set<ProteinModification> potentialModifications) {

		if (chains==null) {
			throw new IllegalArgumentException("Null structure.");
		}

		if (potentialModifications==null) {
			throw new IllegalArgumentException("Null potentialModifications.");
		}


		reset();

		if (potentialModifications.isEmpty()) {
			return;
		}


		residues = new ArrayList<Group>();
		List<Group> ligands = new ArrayList<Group>();
		Map<Component, Set<Group>> mapCompGroups = new HashMap<Component, Set<Group>>();

		for (Chain chain : chains) {

			List<Group> ress = StructureUtil.getAminoAcids(chain);

			//List<Group> ligs = chain.getAtomLigands();
			List<Group> ligs = StructureTools.filterLigands(chain.getAtomGroups());
			residues.addAll(ress);
			residues.removeAll(ligs);
			ligands.addAll(ligs);
			addModificationGroups(potentialModifications, ress, ligs, mapCompGroups);
		}

		if (residues.isEmpty()) {
			String pdbId = "?";
			if ( chains.size() > 0) {
				Structure struc = chains.get(0).getStructure();
				if ( struc != null)
					pdbId = struc.getPDBCode();
			}
			logger.warn("No amino acids found for {}. Either you did not parse the PDB file with alignSEQRES records, or this record does not contain any amino acids.", pdbId);
		}
		List<ModifiedCompound> modComps = new ArrayList<ModifiedCompound>();

		for (ProteinModification mod : potentialModifications) {
			ModificationCondition condition = mod.getCondition();
			List<Component> components = condition.getComponents();
			if (!mapCompGroups.keySet().containsAll(components)) {
				// not all components exist for this mod.
				continue;
			}

			int sizeComps = components.size();
			if (sizeComps==1) {

				processCrosslink1(mapCompGroups, modComps, mod, components);

			} else {

				processMultiCrosslink(mapCompGroups, modComps, mod, condition);
			}
		}

		if (recordAdditionalAttachments) {
			// identify additional groups that are not directly attached to amino acids.
			for (ModifiedCompound mc : modComps) {
				identifyAdditionalAttachments(mc, ligands, chains);
			}
		}

		mergeModComps(modComps);

		identifiedModifiedCompounds.addAll(modComps);


		// record unidentifiable linkage
		if (recordUnidentifiableModifiedCompounds) {
			recordUnidentifiableAtomLinkages(modComps, ligands);
			recordUnidentifiableModifiedResidues(modComps);
		}
	}
1071843
def subalignment(alnfle, subtype, alntype="fasta"):
    
    aln = AlignIO.read(alnfle, alntype)
    alnlen = aln.get_alignment_length()
    nseq = len(aln)
    subaln = None
    subalnfile = alnfle.rsplit(".", 1)[0] + "_{0}.{1}".format(subtype, alntype)

    if subtype == "synonymous":
        for j in range( 0, alnlen, 3 ):
            aa = None
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in CODON_TRANSLATION:
                    break
                if aa and CODON_TRANSLATION[codon] != aa:
                    break
                else:
                    aa = CODON_TRANSLATION[codon]
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subtype == "fourfold":
        for j in range( 0, alnlen, 3 ):
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in FOURFOLD:
                    break
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subaln:
        AlignIO.write(subaln, subalnfile, alntype)
        return subalnfile
    else:
        print("No sites {0} selected.".format(subtype), file=sys.stderr)
        return None
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
652901
def store_psm_protein_relations(fn, header, pgdb, proteins):
    
    # TODO do we need an OrderedDict or is regular dict enough?
    # Sorting for psm_id useful?
    allpsms = OrderedDict()
    last_id, psmids_to_store = None, set()
    store_soon = False
    for psm in tsvreader.generate_tsv_psms(fn, header):
        psm_id, prots = tsvreader.get_pepproteins(psm)
        prots = [x for x in prots if x in proteins]
        try:
            # In case the PSMs are presented unrolled
            allpsms[psm_id].extend(prots)
        except KeyError:
            allpsms[psm_id] = prots
        if len(psmids_to_store) % DB_STORE_CHUNK == 0:
            store_soon = True
        if store_soon and last_id != psm_id:
            pgdb.store_peptides_proteins(allpsms, psmids_to_store)
            store_soon = False
            psmids_to_store = set()
        psmids_to_store.add(psm_id)
        last_id = psm_id
    if len(psmids_to_store) > 0:
        pgdb.store_peptides_proteins(allpsms, psmids_to_store)
    pgdb.index_protein_peptides()
    return allpsms
1089341
def align(fastq_file, pair_file, index_dir, names, align_dir, data):
    
    umi_ext = "-cumi" if "umi_bam" in data else ""
    out_file = os.path.join(align_dir, "{0}-sort{1}.bam".format(dd.get_sample_name(data), umi_ext))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    rg_info = novoalign.get_rg_info(names)
    preset = "sr"

    pair_file = pair_file if pair_file else ""
    if data.get("align_split"):
        final_file = out_file
        out_file, data = alignprep.setup_combine(final_file, data)
        fastq_file, pair_file = alignprep.split_namedpipe_cls(fastq_file, pair_file, data)
    else:
        final_file = None

    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
            index_file = None
            # Skip trying to use indices now as they provide only slight speed-ups
            # and give inconsitent outputs in BAM headers
            # If a single index present, index_dir points to that
            # if index_dir and os.path.isfile(index_dir):
            #     index_dir = os.path.dirname(index_dir)
            #     index_file = os.path.join(index_dir, "%s-%s.mmi" % (dd.get_genome_build(data), preset))
            if not index_file or not os.path.exists(index_file):
                index_file = dd.get_ref_file(data)
            cmd = ("minimap2 -a -x {preset} -R '{rg_info}' -t {num_cores} {index_file} "
                   "{fastq_file} {pair_file} | ")
            do.run(cmd.format(**locals()) + tobam_cl, "minimap2 alignment: %s" % dd.get_sample_name(data))
    data["work_bam"] = out_file
    return data
1160113
def dM(self, t, param, Mt, tips=None, gaps=None):
        """"""
        assert isinstance(t, float) and t > 0, "Invalid t: {0}".format(t)
        assert (param == 't') or (param in self.freeparams), (
                "Invalid param: {0}".format(param))

        if Mt is None:
            Mt = self.M(t, tips=tips, gaps=gaps)

        if (param == 'mu') or (param == 't'):
            if param == 'mu':
                alpha = t
            else:
                alpha = self.mu
            if tips is None:
                dM_param = broadcastMatrixMultiply(self.Prxy, Mt, alpha=alpha)
            else:
                dM_param = broadcastMatrixVectorMultiply(self.Prxy, Mt, alpha=alpha)
                if gaps is not None:
                    dM_param[gaps] = scipy.zeros(N_CODON, dtype='float')
            return dM_param

        paramval = getattr(self, param)
        if isinstance(paramval, float):
            paramisvec = False
        else:
            assert isinstance(paramval, numpy.ndarray) and paramval.ndim == 1
            paramisvec = True
            paramlength = paramval.shape[0]

        if ('expD', t) not in self._cached:
            self._cached[('expD', t)] = scipy.exp(self.D * self.mu * t)
        expD = self._cached[('expD', t)]

        if ('V', t) not in self._cached:
            if 'Dxx_Dyy' not in self._cached:
                Dyy = scipy.tile(self.D, (1, N_CODON)).reshape(
                        self.nsites, N_CODON, N_CODON)
                Dxx = scipy.array([Dyy[r].transpose() for r in
                        range(self.nsites)])
                self._cached['Dxx_Dyy'] = Dxx - Dyy
            Dxx_Dyy = self._cached['Dxx_Dyy']
            if 'Dxx_Dyy_lt_ALMOST_ZERO' not in self._cached:
                self._cached['Dxx_Dyy_lt_ALMOST_ZERO'] = scipy.fabs(
                        Dxx_Dyy) < ALMOST_ZERO
            Dxx_Dyy_lt_ALMOST_ZERO = self._cached['Dxx_Dyy_lt_ALMOST_ZERO']
            with scipy.errstate(divide='raise', under='ignore',
                    over='raise', invalid='ignore'):
                expDyy = scipy.tile(expD,(1, N_CODON)).reshape(
                        self.nsites, N_CODON, N_CODON)
                expDxx = scipy.array([expDyy[r].transpose() for r in
                        range(self.nsites)])
                V = (expDxx - expDyy) / Dxx_Dyy
            with scipy.errstate(under='ignore'): # OK if some values 0
                scipy.copyto(V, self.mu * t * expDxx, where=
                        Dxx_Dyy_lt_ALMOST_ZERO)
            self._cached[('V', t)] = V
        V = self._cached[('V', t)]

        with scipy.errstate(under='ignore'): # don't worry if some values 0
            if tips is None:
                if not paramisvec:
                    dM_param = broadcastMatrixMultiply(self.A,
                            broadcastMatrixMultiply(self.B[param]
                            * V, self.Ainv))
                else:
                    dM_param = scipy.ndarray((paramlength, self.nsites,
                            N_CODON, N_CODON), dtype='float')
                    for j in range(paramlength):
                        dM_param[j] = broadcastMatrixMultiply(self.A,
                                broadcastMatrixMultiply(self.B[param][j]
                                * V, self.Ainv))
            else:
                if not paramisvec:
                    dM_param = broadcastMatrixVectorMultiply(self.A,
                            broadcastGetCols(broadcastMatrixMultiply(
                            self.B[param] * V, self.Ainv), tips))
                else:
                    dM_param = scipy.ndarray((paramlength, self.nsites,
                            N_CODON), dtype='float')
                    for j in range(paramlength):
                        dM_param[j] = broadcastMatrixVectorMultiply(self.A,
                            broadcastGetCols(broadcastMatrixMultiply(
                            self.B[param][j] * V, self.Ainv), tips))
                if gaps is not None:
                    if not paramisvec:
                        dM_param[gaps] = scipy.zeros(N_CODON, dtype='float')
                    else:
                        dM_param[:, gaps] = scipy.zeros(N_CODON, dtype='float')
        return dM_param
1373223
private void identifyAdditionalAttachments(ModifiedCompound mc,
			List<Group> ligands, List<Chain> chains) {
		if (ligands.isEmpty()) {
			return;
		}

		// TODO: should the additional groups only be allowed to the identified
		// ligands or both amino acids and ligands? Currently only on ligands
		// ligands to amino acid bonds for same modification of unknown category
		// will be combined in mergeModComps()
		// TODO: how about chain-chain links?
		List<Group> identifiedGroups = new ArrayList<Group>();
		for (StructureGroup num : mc.getGroups(false)) {
			Group group;
			try {
				//String numIns = "" + num.getResidueNumber();
				//if (num.getInsCode() != null) {
				//	numIns += num.getInsCode();
				//}
				ResidueNumber resNum = new ResidueNumber();
				resNum.setChainName(num.getChainId());
				resNum.setSeqNum(num.getResidueNumber());
				resNum.setInsCode(num.getInsCode());
				//group = chain.getGroupByPDB(numIns);

				group = getGroup(num,chains);
				//group = mapChainIdChain.get(num.getChainId()).getGroupByPDB(resNum);
			} catch (StructureException e) {
				logger.error("Exception: ", e);
				// should not happen
				continue;
			}
			identifiedGroups.add(group);
		}

		int start = 0;

		int n = identifiedGroups.size();
		while (n > start) {
			for (Group group1 : ligands) {
				for (int i=start; i<n; i++) {
					Group group2 = identifiedGroups.get(i);
					if (!identifiedGroups.contains(group1)) {
						List<Atom[]> linkedAtoms = StructureUtil.findAtomLinkages(
								group1, group2, false, bondLengthTolerance);
						if (!linkedAtoms.isEmpty()) {
							for (Atom[] atoms : linkedAtoms) {
								mc.addAtomLinkage(StructureUtil.getStructureAtomLinkage(atoms[0],
										false, atoms[1], false));
							}
							identifiedGroups.add(group1);
							break;
						}
					}
				}
			}

			start = n;
			n = identifiedGroups.size();
		}
	}
12299
def map_alignment(self, prot_seq, nucl_seq):
        
        if prot_seq.id != nucl_seq.id:
            logging.warning(
                'ID mismatch: %s != %s. Are the sequences in the same order?',
                prot_seq.id, nucl_seq.id)

        # Ungap nucleotides
        codons = batch(str(nucl_seq.seq.ungap('-')), 3)
        codons = [''.join(i) for i in codons]
        codon_iter = iter(codons)

        ungapped_prot = str(prot_seq.seq).replace('-', '')

        if len(ungapped_prot) != len(codons):
            table = self.translation_table.forward_table
            prot_str = ' '.join(' ' + p + ' ' for p in ungapped_prot)
            codon_str = ' '.join(codons)
            trans_str = ' '.join(' ' + table.get(codon, 'X') + ' '
                                 for codon in codons)
            raise ValueError("""Length of codon sequence ({0}) does not match \
length of protein sequence ({1}) for {2}
Protein:       {3}
Codons:        {4}
Trans. Codons: {5}""".format(len(codons), len(ungapped_prot), nucl_seq.id, prot_str,
                            codon_str, trans_str))

        try:
            nucl_align = ['---' if p == '-' else next(codon_iter)
                          for p in str(prot_seq.seq)]
        except StopIteration:
            assert False  # Should be checked above

        result = SeqRecord(Seq(''.join(nucl_align)), id=nucl_seq.id,
                           description=nucl_seq.description)

        # Validate
        self._validate_translation(prot_seq.seq.upper(), result.seq.upper())

        return result
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
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)
462912
def find_local_alignment(hit, prediction, hsp)
      # indexing in blast starts from 1
      hit_local   = hit.raw_sequence[hsp.hit_from - 1..hsp.hit_to - 1]
      query_local = prediction.raw_sequence[hsp.match_query_from -
                                            1..hsp.match_query_to - 1]

      # in case of nucleotide prediction sequence translate into protein
      # use translate with reading frame 1 because
      # to/from coordinates of the hsp already correspond to the
      # reading frame in which the prediction was read to match this hsp
      if @type == :nucleotide
        s = Bio::Sequence::NA.new(query_local)
        query_local = s.translate
      end

      opt = ['--maxiterate', '1000', '--localpair', '--anysymbol', '--quiet',
             '--thread', @num_threads.to_s]
      mafft = Bio::MAFFT.new('mafft', opt)

      # local alignment for hit and query
      seqs = [hit_local, query_local]
      report = mafft.query_align(seqs)
      report.alignment.map(&:to_s)
    rescue StandardError
      raise NoMafftInstallationError
    end