reads in an amino acid substitution matrix, if necessary
185850
private static SubstitutionMatrix<AminoAcidCompound> getAminoAcidMatrix(String file) {
		if (!aminoAcidMatrices.containsKey(file)) {
			aminoAcidMatrices.put(file, new SimpleSubstitutionMatrix<AminoAcidCompound>(
					AminoAcidCompoundSet.getAminoAcidCompoundSet(), getReader(file), file));
		}
		return aminoAcidMatrices.get(file);
	}
450589
def read_from_file(cls, filename):
        """"""
        with open(filename) as f:
            lines = list(line for line in f if line[0] != '#')
        r = []
        t = []
        for line in lines[:3]:
            values = list(float(word) for word in line.split())
            r.append(values[:3])
            t.append(values[3])
        transformation = Complete(r, t)
        affected_atoms = set(int(word) for word in lines[3].split())
        return cls(affected_atoms, transformation)
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]
1207364
def _read_gtf(gtf):
    
    if not gtf:
        return gtf
    db = defaultdict(list)
    with open(gtf) as in_handle:
        for line in in_handle:
            if line.startswith("#"):
                continue
            cols = line.strip().split("\t")
            name = [n.split("=")[1] for n in cols[-1].split(";") if n.startswith("Name")]
            chrom, start, end, strand = cols[0], cols[3], cols[4], cols[6]
            if cols[2] == "miRNA_primary_transcript":
                db[name[0]].append([chrom, int(start), int(end), strand])
    return db
62596
def alignprotein(self, sample, target):
        
        # Initialise dictionaries
        sample[self.analysistype].dnaseq = dict()
        sample[self.analysistype].protseq = dict()
        sample[self.analysistype].ntindex = dict()
        sample[self.analysistype].aaindex = dict()
        sample[self.analysistype].ntalign = dict()
        sample[self.analysistype].aaalign = dict()
        sample[self.analysistype].aaidentity = dict()
        # Remove any gaps incorporated into the sequence
        sample[self.analysistype].targetsequence[target] = \
            sample[self.analysistype].targetsequence[target].replace('-', '')
        # In order to properly translate the nucleotide sequence, BioPython requests that the sequence is a multiple of
        # three - not partial codons. Trim the sequence accordingly
        remainder = 0 - len(sample[self.analysistype].targetsequence[target]) % 3
        seq = sample[self.analysistype].targetsequence[target] if remainder == 0 \
            else sample[self.analysistype].targetsequence[target][:remainder]
        # Set the DNA and protein sequences of the target in the sample
        sample[self.analysistype].dnaseq[target] = Seq(seq, IUPAC.unambiguous_dna)
        # Translate the nucleotide sequence
        sample[self.analysistype].protseq[target] = str(sample[self.analysistype].dnaseq[target].translate())
        for targetfile in self.targetfiles:
            # Trim the reference sequence to multiples of three
            refremainder = 0 - len(self.records[targetfile][target].seq) % 3
            refseq = str(self.records[targetfile][target].seq) if refremainder % 3 == 0 \
                else str(self.records[targetfile][target].seq)[:refremainder]
            # Translate the nucleotide sequence of the reference sequence
            refdna = Seq(refseq, IUPAC.unambiguous_dna)
            refprot = str(refdna.translate())
            # Use pairwise2 to perform a local alignment with the following parameters:
            # x     No match parameters. Identical characters have score of 1, otherwise 0.
            # s     Same open (-1)  and extend (-.1) gap penalties for both sequences
            ntalignments = pairwise2.align.localxs(seq, refseq, -1, -.1)
            # Use format_alignment to create a formatted alignment that is subsequently split on newlines e.g.
            '''
            ACCGT
            | ||
            A-CG-
            Score=3
            '''
            ntformat = (str(format_alignment(*ntalignments[0])).split('\n'))
            # Align the nucleotide sequence of the reference (ntalignments[2]) to the sample (ntalignments[0]).
            # If the corresponding bases match, add a |, otherwise a space
            ntalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(ntformat[0], ntformat[2])))
            # Create the nucleotide alignment: the sample sequence, the (mis)matches, and the reference sequence
            sample[self.analysistype].ntalign[target] = self.interleaveblastresults(ntformat[0], ntformat[2])
            # Regex to determine location of mismatches in the sequences
            count = 0
            sample[self.analysistype].ntindex[target] = str()
            for snp in re.finditer(' ', ntalignment):
                # If there are many SNPs, then insert line breaks for every 10 SNPs
                if count <= 10:
                    sample[self.analysistype].ntindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].ntindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1
            # Perform the same steps, except for the amino acid sequence
            aaalignments = pairwise2.align.localxs(sample[self.analysistype].protseq[target], refprot, -1, -.1)
            aaformat = (str(format_alignment(*aaalignments[0])).split('\n'))
            aaalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(aaformat[0], aaformat[2])))
            sample[self.analysistype].aaidentity[target] = '{:.2f}'\
                .format(float(aaalignment.count('|')) / float(len(aaalignment)) * 100)
            sample[self.analysistype].aaalign[target] = self.interleaveblastresults(aaformat[0], aaformat[2])
            count = 0
            sample[self.analysistype].aaindex[target] = str()
            for snp in re.finditer(' ', aaalignment):
                if count <= 10:
                    sample[self.analysistype].aaindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].aaindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1
1078224
def _load_data_index(self, chain1, chain2):
        
        parser = PDB.PDBParser(QUIET=True)

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

        residues1 = list(structure1.get_residues())
        residues2 = list(structure2.get_residues())

        indexes1 = set(r.id[1] for r in residues1)
        indexes2 = set(r.id[1] for r in residues2)

        indexes = indexes1.intersection(indexes2)
        self.indexes = indexes.copy()
        self.N = len(indexes)

        coord1 = []
        indexes1 = indexes.copy()
        for r in residues1:
            if r.id[1] in indexes1 and 'CA' in r:
                coord1.append(np.concatenate((r['CA'].get_coord(), (1,)))[:, None])
                # Remove from index to avoid repeated residues
                indexes1.remove(r.id[1])
        coord1 = np.hstack(coord1).astype(DTYPE, copy=False)

        coord2 = []
        for r in residues2:
            if r.id[1] in indexes and 'CA' in r:
                coord2.append(np.concatenate((r['CA'].get_coord(), (1,)))[:, None])
                indexes.remove(r.id[1])
        coord2 = np.hstack(coord2).astype(DTYPE, copy=False)

        self.coord1 = coord1
        self.coord2 = coord2
429033
def read_RNAplfold(tmpdir, maxlen=None, seq_align="start", pad_with="E"):
    
    assert pad_with in {"P", "H", "I", "M", "E"}

    def read_profile(tmpdir, P):
        return [values.strip().split("\t")
                for seq_name, values in iter_fasta("{tmp}/{P}_profile.fa".format(tmp=tmpdir, P=P))]

    def nelem(P, pad_width):
        """get the right neutral element
        """
        return 1 if P is pad_with else 0

    arr_hime = np.array([pad_sequences(read_profile(tmpdir, P),
                                       value=[nelem(P, pad_with)],
                                       align=seq_align,
                                       maxlen=maxlen)
                         for P in RNAplfold_PROFILES_EXECUTE], dtype="float32")

    # add the pairness column
    arr_p = 1 - arr_hime.sum(axis=0)[np.newaxis]
    arr = np.concatenate((arr_p, arr_hime))

    # reshape to: seq, seq_length, num_channels
    arr = np.moveaxis(arr, 0, 2)
    return arr
1200284
def _read_amino_acids():
    """"""
    this_dir = os.path.dirname(os.path.abspath(__file__))
    aa_file = os.path.join(this_dir, os.pardir, 'resources', 'amino_acids.tsv')
    amino_acids = {}
    amino_acids_reverse = {}
    with open(aa_file, 'rt') as fh:
        lines = fh.readlines()
    for lin in lines[1:]:
        terms = lin.strip().split('\t')
        key = terms[2]
        val = {'full_name': terms[0],
               'short_name': terms[1],
               'indra_name': terms[3]}
        amino_acids[key] = val
        for v in val.values():
            amino_acids_reverse[v] = key
    return amino_acids, amino_acids_reverse
245706
def parseFullScan(self, i, modifications=False):
        
        scanObj = PeptideObject()
        peptide = str(i[1])
        pid=i[2]
        scanObj.acc = self.protein_map.get(i[4], i[4])
        if pid is None:
            return None
        if modifications:
            sql = 'select aam.ModificationName,pam.Position,aam.DeltaMass from peptidesaminoacidmodifications pam left join aminoacidmodifications aam on (aam.AminoAcidModificationID=pam.AminoAcidModificationID) where pam.PeptideID=%s'%pid
            for row in self.conn.execute(sql):
                scanObj.addModification(peptide[row[1]], str(row[1]), str(row[2]), row[0])
        else:
            mods = self.mods.get(int(pid))
            if mods is not None:
                for modId, modPosition in zip(mods[0].split(','),mods[1].split(',')):
                    modEntry = self.modTable[str(modId)]
                    scanObj.addModification(peptide[int(modPosition)], modPosition, modEntry[1], modEntry[0])
            tmods = self.tmods.get(int(pid))
            if tmods is not None:
                for modIds in tmods:
                    for modId in modIds.split(','):
                        modEntry = self.modTable[str(modId)]
                        scanObj.addModification('[', 0, modEntry[1], modEntry[0])
        scanObj.peptide = peptide
        if self.decompressScanInfo(scanObj, i[0]):
            return scanObj
        return None
68613
private void initMatrix(String line) {
		String[] spl = line.split(" ");

		// trim off the final , character
		currentRows = spl[3].substring(0, spl[3].length()-1);
		currentCols = spl[6];
		currentRowPos = -1;

		int nrRows = currentRows.length();
		int nrCols = currentCols.length();

		matrix = new short[nrRows][nrCols];

		rows = new ArrayList<AminoAcidCompound>();
		cols = new ArrayList<AminoAcidCompound>();


		//System.out.println(">" + currentRows+"<");
		AminoAcidCompoundSet compoundSet = AminoAcidCompoundSet.getAminoAcidCompoundSet();
		for ( int i = 0 ; i < currentRows.length() ; i ++){
			char c = currentRows.charAt(i);
			AminoAcidCompound aa = compoundSet.getCompoundForString(String.valueOf(c));

			rows.add(aa);
		}

		for ( int i = 0 ; i < currentCols.length() ; i ++){
			char c = currentRows.charAt(i);
			AminoAcidCompound aa = compoundSet.getCompoundForString(String.valueOf(c));

			cols.add(aa);
		}





		currentMatrix.setScale(scale);
	}
1064830
def _fileToMatrix(file_name):
    """"""
    # TODO: np.loadtxt() might be an alternative
    #     try:
    if 1 < 3:
        lres = []
        for line in open(file_name, 'r').readlines():
            if len(line) > 0 and line[0] not in ('%', '#'):
                lres.append(list(map(float, line.split())))
        res = lres
    while res != [] and res[0] == []:  # remove further leading empty lines
        del res[0]
    return res
    #     except:
    print('could not read file ' + file_name)
1064025
private void assembleLinkages(List<List<Atom[]>> matchedAtomsOfLinkages,
			ProteinModification mod, List<ModifiedCompound> ret) {
		ModificationCondition condition = mod.getCondition();
		List<ModificationLinkage> modLinks = condition.getLinkages();

		int nLink = matchedAtomsOfLinkages.size();
		int[] indices = new int[nLink];
		Set<ModifiedCompound> identifiedCompounds = new HashSet<ModifiedCompound>();
		while (indices[0]<matchedAtomsOfLinkages.get(0).size()) {
			List<Atom[]> atomLinkages = new ArrayList<Atom[]>(nLink);
			for (int iLink=0; iLink<nLink; iLink++) {
				Atom[] atoms = matchedAtomsOfLinkages.get(iLink).get(indices[iLink]);
				atomLinkages.add(atoms);
			}
			if (matchLinkages(modLinks, atomLinkages)) {
				// matched

				int n = atomLinkages.size();
				List<StructureAtomLinkage> linkages = new ArrayList<StructureAtomLinkage>(n);
				for (int i=0; i<n; i++) {
					Atom[] linkage = atomLinkages.get(i);
					StructureAtomLinkage link = StructureUtil.getStructureAtomLinkage(
							linkage[0], residues.contains(linkage[0].getGroup()),
							linkage[1], residues.contains(linkage[1].getGroup()));
					linkages.add(link);
				}

				ModifiedCompound mc = new ModifiedCompoundImpl(mod, linkages);
				if (!identifiedCompounds.contains(mc)) {
					ret.add(mc);
					identifiedCompounds.add(mc);
				}
			}

			// indices++ (e.g. [0,0,1]=>[0,0,2]=>[1,2,0])
			int i = nLink-1;
			while (i>=0) {
				if (i==0 || indices[i]<matchedAtomsOfLinkages.get(i).size()-1) {
					indices[i]++;
					break;
				} else {
					indices[i] = 0;
					i--;
				}
			}
		}
	}
1046092
def _incorporate_delins(self):
        """"""
        seq, cds_start, cds_stop, start, end = self._setup_incorporate()

        ref = self._var_c.posedit.edit.ref
        alt = self._var_c.posedit.edit.alt
        ref_length = end - start if ref is not None else 0    # can't just get from ref since ref isn't always known
        alt_length = len(
            self._var_c.posedit.edit.alt) if self._var_c.posedit.edit.alt is not None else 0
        net_base_change = alt_length - ref_length
        cds_stop += net_base_change

        # incorporate the variant into the sequence (depending on the type)
        is_substitution = False
        if ref is not None and alt is not None:    # delins or SNP
            seq[start:end] = list(alt)
            if len(ref) == 1 and len(alt) == 1:
                is_substitution = True
        elif ref is not None:    # deletion
            del seq[start:end]
        else:    # insertion
            seq[start + 1:start + 1] = list(alt)    # insertion in list before python list index

        if DBG:
            print("net base change: {}".format(net_base_change))
        is_frameshift = net_base_change % 3 != 0
        # use max of mod 3 value and 1 (in event that indel starts in the 5'utr range)
        variant_start_aa = max(int(math.ceil((self._var_c.posedit.pos.start.base) / 3.0)), 1)

        alt_data = AltTranscriptData(
            seq,
            cds_start,
            cds_stop,
            is_frameshift,
            variant_start_aa,
            self._transcript_data.protein_accession,
            is_substitution=is_substitution,
            is_ambiguous=self._ref_has_multiple_stops)
        return alt_data
1041065
def read_alignments(self, input_dir, file_format, header_grep=None):
        

        compression = ['', 'gz', 'bz2']

        if file_format == 'fasta':
            extensions = ['fa', 'fas', 'fasta']

        elif file_format == 'phylip':
            extensions = ['phy']

        else:
            extensions = []

        extensions = list('.'.join([x,y]) if y else x for x,y in itertools.product(extensions, compression))

        files = fileIO.glob_by_extensions(input_dir, extensions)
        files.sort(key=SORT_KEY)
        self._input_files = files
        records = []

        if self.show_progress:
            pbar = setup_progressbar("Loading files", len(files), simple_progress=True)
            pbar.start()

        for i, f in enumerate(files):
            if f.endswith('.gz') or f.endswith('.bz2'):
                fd, tmpfile = tempfile.mkstemp()
                with fileIO.freader(f, f.endswith('.gz'), f.endswith('.bz2')) as reader,\
                     fileIO.fwriter(tmpfile) as writer:
                    for line in reader:
                        if ISPY3:
                            line = line.decode()
                        writer.write(line)
                try:
                    record = Alignment(tmpfile, file_format, True)
                except ValueError:
                    record = Alignment(tmpfile, file_format, False)
                finally:
                    os.close(fd)
                    os.unlink(tmpfile)

            else:
                try:
                    record = Alignment(f, file_format, True)
                except RuntimeError:
                    record = Alignment(f, file_format, False)

            if header_grep:
                try:
                    datatype = 'dna' if record.is_dna() else 'protein'

                    record = Alignment([(header_grep(x), y) for (x, y) in record.get_sequences()], datatype)

                except TypeError:
                    raise TypeError("Couldn't apply header_grep to header\n"
                                    "alignment number={}, name={}\n"
                                    "header_grep={}".format(i, fileIO.strip_extensions(f), header_grep))
                except RuntimeError:
                    print('RuntimeError occurred processing alignment number={}, name={}'
                          .format(i, fileIO.strip_extensions(f)))
                    raise

            record.name = (fileIO.strip_extensions(f))
            records.append(record)
            if self.show_progress:
                pbar.update(i)
        if self.show_progress:
            pbar.finish()
        return records
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file
573221
public SubstitutionMatrix<AminoAcidCompound> getSubstitutionMatrix() {
		if ( substitutionMatrix == null){
			substitutionMatrix = SubstitutionMatrixHelper.getMatrixFromAAINDEX(DEFAULT_SUBSTITUTION_MATRIX);

		}
		return substitutionMatrix;
	}
1323543
def run_all2(protgroup, memornot, subsequences, base_outdir,
            protgroup_dict, protein_feathers_dir, date, errfile, impute_counts=True,
            cutoff_num_proteins=0, core_only_genes=None,
            length_filter_pid=.8, remove_correlated_feats=True,
            force_rerun_counts=False, force_rerun_percentages=False, force_rerun_pca=False):
    """"""
    import ssbio.utils

    # Need to set multiprocessing limit for scipy/numpy stuff if parallelizing anything
    import os
    os.environ['OMP_NUM_THREADS'] = '1'

    # First, filter down the protein group to the membrane/nonmembrane definition
    prots_filtered_feathers = get_protein_feather_paths(protgroup=protgroup, memornot=memornot,
                                                        protgroup_dict=protgroup_dict,
                                                        protein_feathers_dir=protein_feathers_dir,
                                                        core_only_genes=core_only_genes)
    num_proteins = len(prots_filtered_feathers)
    if num_proteins <= cutoff_num_proteins:
        return

    # Make output directories
    protscale = 'proteome_unscaled'
    outdir_d0 = ssbio.utils.make_dir(op.join(base_outdir, protscale))
    outdir_d1 = ssbio.utils.make_dir(op.join(outdir_d0, '-'.join(memornot)))
    outdir_final = ssbio.utils.make_dir(op.join(outdir_d1, '-'.join(protgroup)))

    if impute_counts:
        big_strain_counts_df = get_proteome_counts_impute_missing(prots_filtered_feathers=prots_filtered_feathers,
                                                                  outpath=op.join(outdir_final,
                                                                                  '{}-subsequence_proteome_IMP.fthr'.format(
                                                                                      date)),
                                                                  length_filter_pid=length_filter_pid,
                                                                  force_rerun=force_rerun_counts)

        big_strain_percents_df = get_proteome_percentages(counts_df=big_strain_counts_df,
                                                          outpath=op.join(outdir_final,
                                                                          '{}-subsequence_proteome_perc_IMP.fthr'.format(
                                                                                  date)),
                                                          force_rerun=force_rerun_percentages)
        pca_pickle = op.join(outdir_final, '{}-subsequence_pca.pckl'.format(date))
    # Divide by totals to get percentages in a new dataframe
    else:
        try:
            big_strain_percents_df = get_proteome_correct_percentages(prots_filtered_feathers=prots_filtered_feathers,
                                                                      outpath=op.join(outdir_final,
                                                                                      '{}-subsequence_proteome_perc_AVG.fthr'.format(
                                                                                              date)),
                                                                      length_filter_pid=length_filter_pid,
                                                                      force_rerun=force_rerun_percentages)
            pca_pickle = op.join(outdir_final, '{}-subsequence_pca_AVG.pckl'.format(date))
        except:
            with open(errfile, "a") as myfile:
                myfile.write('PERCENTAGES ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return



    if ssbio.utils.force_rerun(flag=force_rerun_pca, outfile=pca_pickle):

        # Then, get filters for rows of the loaded feathers for interested subsequences
        keep_subsequences = get_interested_subsequences(subsequences=subsequences)

        # Some numbers: number of features
        num_feats = len(big_strain_percents_df)

        # Make an unwieldy title
        big_title = 'LOC={0}; PROTGROUP={1};\n' \
                    'NUMPROTS={2}; NUMFEATS={3}'.format('-'.join(memornot),
                                                                        '-'.join(protgroup),
                                                                        num_proteins,
                                                                        num_feats)

        # Run PCA and make plots
        runner = PCAMultiROS(features_df=big_strain_percents_df, observations_df=pd.DataFrame(), plot_title=big_title)
        try:
            runner.clean_data(keep_features=keep_subsequences, remove_correlated_feats=remove_correlated_feats)
        except:
            with open(errfile, "a") as myfile:
                myfile.write(
                    'CLEAN ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return
        # try:
        runner.run_pca()
        # except:
        #     with open(errfile, "a") as myfile:
        #         myfile.write(
        #             'PCA ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
        #     return
        with open(pca_pickle, 'wb') as f:
            pickle.dump(runner, f)
    else:
        with open(pca_pickle, 'rb') as f:
            runner = pickle.load(f)