taxid=559292
    del_exp_syss=["Co-purification", "Co-fractionation", "Proximity Label-MS", 
                                  "Affinity Capture-RNA", "Protein-peptide"]
339448
def get_dbiogrid_intmap(taxid,dbiogridp,dbiogrid_intmap_symmp,dbiogrid_intlinp,dbiogrid_intmapp,
                        logf=None,
                        experimental_system_type='physical',
                        genefmt='name',
                        force=False,test=False,
                        keep_exp_syss=[],del_exp_syss=[],
                        filldiagonal_withna=False):
        
    gene_fmt2colns={'biogrid':{'id':'Systematic Name',
                 'name':'Official Symbol'},}
    if (not exists(dbiogrid_intmap_symmp)) or force:
        if not exists(dbiogrid_intmapp) or force:
            if not exists(dbiogrid_intlinp) or force:
                if dbiogridp.endswith('tab2.txt'):
                    print(f"converting to parquet")
                    dbiogrid=pd.read_csv(dbiogridp,
                       sep='\t',low_memory=False)
                    dbiogridp=dbiogridp+".pqt"
                    to_table_pqt(dbiogrid,dbiogridp)
                    print(f"use {dbiogridp}")
                elif dbiogridp.endswith('tab2.txt.pqt'):
                    dbiogrid=read_table_pqt(dbiogridp)
                # filter biogrid
                # taxonomic id of Scer
                dbiogrid=dbiogrid.loc[((dbiogrid['Experimental System Type']==experimental_system_type) \
                           & (dbiogrid['Organism Interactor A']==taxid) \
                           & (dbiogrid['Organism Interactor B']==taxid)),:]

                print('All physical Experimental Systems.')
                if test:
                    dlog=pd.DataFrame(dbiogrid['Experimental System'].value_counts())
                    print(dlog)
                    if not logf is None:
                        to_table(dlog,f'{logf}.all_int.tsv')                    
                if len(keep_exp_syss)!=0:
                    dbiogrid=dbiogrid.loc[((dbiogrid['Experimental System'].isin(keep_exp_syss))),:]
                elif len(del_exp_syss)!=0:
                    dbiogrid=dbiogrid.loc[((~dbiogrid['Experimental System'].isin(del_exp_syss))),:]
                elif len(del_exp_syss)!=0 and len(keep_exp_syss)!=0:
                    print('Error: either specify keep_exp_syss or del_exp_syss')
                    return False
                if test:
                    print('Experimental Systems used.')
                    dlog=pd.DataFrame(dbiogrid['Experimental System'].value_counts())
                    if not logf is None:
                        to_table(dlog,f'{logf}.kept_int.tsv')
                    print(dlog)
                to_table_pqt(dbiogrid,dbiogrid_intlinp)
            else:
                dbiogrid=read_table_pqt(dbiogrid_intlinp)
            # this cell makes a symmetric interaction map
            dbiogrid['count']=1
            if test:
                print(dbiogrid['count'].sum())
            # mean of counts of intearations
            ## higher is more confident interaction captured in multiple assays
            dbiogrid_grid=dbiogrid.pivot_table(values='count',
                                             index=f"{gene_fmt2colns['biogrid'][genefmt]} Interactor A",
                                            columns=f"{gene_fmt2colns['biogrid'][genefmt]} Interactor B",
                                            aggfunc='sum',)

            # make it symmetric
            if test:
                print('shape of non-symm intmap: ',dbiogrid_grid.shape)
            to_table_pqt(dbiogrid_grid,dbiogrid_intmapp)
        else:         
#             dbiogrid_grid=pd.read_table(dbiogrid_intmapp)
            dbiogrid_grid=read_table_pqt(dbiogrid_intmapp+'.pqt')
        dbiogrid_grid=set_index(dbiogrid_grid,f"{gene_fmt2colns['biogrid'][genefmt]} Interactor A")
        geneids=set(dbiogrid_grid.index).union(set(dbiogrid_grid.columns))
        if test:
            print('total number of genes',len(geneids))
        # missing rows to nan
        dbiogrid_intmap_symm=pd.DataFrame(columns=geneids,index=geneids)

        dbiogrid_intmap_symm.loc[dbiogrid_grid.index,:]=dbiogrid_grid.loc[dbiogrid_grid.index,:]
        dbiogrid_intmap_symm.loc[:,dbiogrid_grid.columns]=dbiogrid_grid.loc[:,dbiogrid_grid.columns]
        if test:
            print(dbiogrid_intmap_symm.shape)
        dbiogrid_intmap_symm=dbiogrid_intmap_symm.fillna(0)
        dbiogrid_intmap_symm=(dbiogrid_intmap_symm+dbiogrid_intmap_symm.T)/2
        dbiogrid_intmap_symm.index.name='Interactor A'
        dbiogrid_intmap_symm.columns.name='Interactor B'
#         if test:
#             dbiogrid_intmap_symm=dbiogrid_intmap_symm.iloc[:5,:5]
        if filldiagonal_withna:
            dbiogrid_intmap_symm=filldiagonal(dbiogrid_intmap_symm)
        to_table_pqt(dbiogrid_intmap_symm,dbiogrid_intmap_symmp)
        print('file saved at: ',dbiogrid_intmap_symmp)
        dbiogrid_intmap_symm_lin=get_degrees(dbiogrid_intmap_symm)
        dbiogrid_intmap_symm_lin.index.name=f'gene {genefmt}'
        to_table(dbiogrid_intmap_symm_lin,f'{dbiogrid_intmap_symmp}.degrees.tsv')
    else:
        dbiogrid_intmap_symm=read_table_pqt(dbiogrid_intmap_symmp)
    return dbiogrid_intmap_symm
984494
def fasttagcount(sam, out, genemap, positional, minevidence, cb_histogram, 
                 cb_cutoff, subsample, parse_tags, gene_tags, umi_matrix):
    ''' 
    '''
    from pysam import AlignmentFile

    from io import StringIO
    import pandas as pd

    from utils import weigh_evidence

    if sam.endswith(".sam"):
        logger.error("To use the fasttagcount subcommand, the alignment file must be a "
                     "coordinate sorted, indexed BAM file.")
        sys.exit(1)

    logger.info('Reading optional files')

    gene_map = None
    if genemap:
        with open(genemap) as fh:
            try:
                gene_map = dict(p.strip().split() for p in fh)
            except ValueError:
                logger.error('Incorrectly formatted gene_map, need to be tsv.')
                sys.exit()

    if positional:
        tuple_template = '{0},{1},{2},{3}'
    else:
        tuple_template = '{0},{1},{3}'

    if not cb_cutoff:
        cb_cutoff = 0

    if cb_histogram and cb_cutoff == "auto":
        cb_cutoff = guess_depth_cutoff(cb_histogram)

    cb_cutoff = int(cb_cutoff)

    cb_hist = None
    filter_cb = False
    if cb_histogram:
        cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
        total_num_cbs = cb_hist.shape[0]
        cb_hist = cb_hist[cb_hist > cb_cutoff]
        logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
        filter_cb = True

    parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')

    if subsample:
        logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
        start_sampling  = time.time()

        reservoir = collections.defaultdict(list)
        cb_hist_sampled = 0 * cb_hist
        cb_obs = 0 * cb_hist

        track = stream_bamfile(sam)
        current_read = 'none_observed_yet'
        for i, aln in enumerate(track):
            if aln.qname == current_read:
                continue

            current_read = aln.qname

            if parse_tags:
                CB = aln.get_tag('CR')
            else:
                match = parser_re.match(aln.qname)
                CB = match.group('CB')

            if CB not in cb_hist.index:
                continue

            cb_obs[CB] += 1
            if len(reservoir[CB]) < subsample:
                reservoir[CB].append(i)
                cb_hist_sampled[CB] += 1
            else:
                s = pd.np.random.randint(0, cb_obs[CB])
                if s < subsample:
                    reservoir[CB][s] = i

        index_filter = set(itertools.chain.from_iterable(reservoir.values()))
        sam_file.close()
        sampling_time = time.time() - start_sampling
        logger.info('Sampling done - {:.3}s'.format(sampling_time))

    evidence = collections.defaultdict(lambda: collections.defaultdict(float))
    bare_evidence = collections.defaultdict(float)
    logger.info('Tallying evidence')
    start_tally = time.time()

    sam_mode = 'r' if sam.endswith(".sam") else 'rb'
    sam_file = AlignmentFile(sam, mode=sam_mode)
    transcript_map = collections.defaultdict(set)
    sam_transcripts = [x["SN"] for x in sam_file.header["SQ"]]
    if gene_map:
        for transcript, gene in gene_map.items():
            if transcript in sam_transcripts:
                transcript_map[gene].add(transcript)
    else:
        for transcript in sam_transcripts:
            transcript_map[transcript].add(transcript)
    missing_transcripts = set()
    alignments_processed = 0
    unmapped = 0
    kept = 0
    nomatchcb = 0
    current_read = 'none_observed_yet'
    current_transcript = None
    count_this_read = True
    transcripts_processed = 0
    genes_processed = 0
    cells = list(cb_hist.index)
    targets_seen = set()

    if umi_matrix:
        bare_evidence_handle = open(umi_matrix, "w")
        bare_evidence_handle.write(",".join(["gene"] + cells) + "\n")

    with open(out, "w") as out_handle:
        out_handle.write(",".join(["gene"] + cells) + "\n")
        for gene, transcripts in transcript_map.items():
            for transcript in transcripts:
                for aln in sam_file.fetch(transcript):
                    alignments_processed += 1

                    if aln.is_unmapped:
                        unmapped += 1
                        continue

                    if gene_tags and not aln.has_tag('GX'):
                        unmapped += 1
                        continue

                    if aln.qname != current_read:
                        current_read = aln.qname
                        if subsample and i not in index_filter:
                            count_this_read = False
                            continue
                        else:
                            count_this_read = True
                    else:
                        if not count_this_read:
                            continue

                    if parse_tags:
                        CB = aln.get_tag('CR')
                    else:
                        match = parser_re.match(aln.qname)
                        CB = match.group('CB')

                    if filter_cb:
                        if CB not in cb_hist.index:
                            nomatchcb += 1
                            continue

                    if parse_tags:
                        MB = aln.get_tag('UM')
                    else:
                        MB = match.group('MB')

                    if gene_tags:
                        target_name = aln.get_tag('GX').split(',')[0]
                    else:
                        txid = sam_file.getrname(aln.reference_id)
                        if gene_map:
                            if txid in gene_map:
                                target_name = gene_map[txid]
                            else:
                                missing_transcripts.add(txid)
                                continue
                        else:
                            target_name = txid
                    targets_seen.add(target_name)

                    # Scale evidence by number of hits
                    evidence[CB][MB] += weigh_evidence(aln.tags)
                    bare_evidence[CB] += weigh_evidence(aln.tags)
                    kept += 1
                transcripts_processed += 1
                if not transcripts_processed % 1000:
                    logger.info("%d genes processed." % genes_processed)
                    logger.info("%d transcripts processed." % transcripts_processed)
                    logger.info("%d alignments processed." % alignments_processed)

            earray = []
            for cell in cells:
                umis = [1 for _, v in evidence[cell].items() if v >= minevidence]
                earray.append(str(sum(umis)))
            out_handle.write(",".join([gene] + earray) + "\n")
            earray = []
            if umi_matrix:
                for cell in cells:
                    earray.append(str(int(bare_evidence[cell])))
                bare_evidence_handle.write(",".join([gene] + earray) + "\n")

            evidence = collections.defaultdict(lambda: collections.defaultdict(int))
            bare_evidence = collections.defaultdict(int)
            genes_processed += 1

    if umi_matrix:
        bare_evidence_handle.close()

    # fill dataframe with missing values, sort and output
    df = pd.read_csv(out, index_col=0, header=0)
    targets = pd.Series(index=set(transcript_map.keys()))
    targets = targets.sort_index()
    df = df.reindex(targets.index.values, fill_value=0)
    df = df.sort_index()
    df.to_csv(out)

    if umi_matrix:
        df = pd.read_csv(umi_matrix, index_col=0, header=0)
        df = df.reindex(targets.index.values, fill_value=0)
        df = df.sort_index()
        df.to_csv(umi_matrix)
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
496476
def convert(data, _first)
      data
        .reject { |o| o['refseq_protein_ids'].empty? }
        .map do |o|
          "#{o['peptide']}\tref|#{o['refseq_protein_ids']}|\t100\t10\t0\t0\t0\t10\t0\t10\t1e-100\t100\n"
        end
        .join
    end
433646
def kpp(self) -> str:
        
        tax_codes = [
            '7700', '7800', '5000', '0100',
            '0200', '0300', '0500', '0600',
            '0700', '0800', '0900', '1000',
            '1100', '1200', '1300', '1400',
            '1500', '1600', '1700', '1800',
            '1900', '2000', '2100', '2200',
            '2300', '2400', '2500', '2600',
            '2700', '2800', '2900', '3000',
            '3100', '3200', '3300', '3400',
            '3500', '3600', '3700', '3800',
            '3900', '4000', '4100', '4900',
            '5100', '5200', '5300', '5400',
            '5500', '5600', '5700', '5800',
            '5900', '6000', '6100', '6200',
            '6300', '6400', '6500', '6600',
            '6700', '6800', '6900', '7000',
            '7100', '7200', '7300', '7400',
            '7500', '7600', '7900', '8600',
            '8700', '8900', '9100', '9200',
            '9800', '9900', '9901', '9951',
            '9952', '9953', '9954', '9955',
            '9956', '9957', '9958', '9959',
            '9961', '9962', '9965', '9966',
            '9971', '9972', '9973', '9974',
            '9975', '9976', '9977', '9979',
            '9998',
        ]

        tax_code = tax_codes[self.random.randint(0, len(tax_codes) - 1)]
        reg_code = '{:02}'.format(self.random.randint(1, 99))
        reg_number = '{:03}'.format(self.random.randint(1, 999))
        kpp = tax_code + reg_code + reg_number
        return kpp
191635
def find_databases(databases):
    
    # 16 ribosomal proteins in their expected order
    proteins = ['L15', 'L18', 'L6', 'S8', 'L5', 'L24', 'L14',
            'S17', 'L16', 'S3', 'L22', 'S19', 'L2', 'L4', 'L3', 'S10']
    # curated databases
    protein_databases = {
                'L14': 'rpL14_JGI_MDM.filtered.faa',
                'L15': 'rpL15_JGI_MDM.filtered.faa',
                'L16': 'rpL16_JGI_MDM.filtered.faa',
                'L18': 'rpL18_JGI_MDM.filtered.faa',
                'L22': 'rpL22_JGI_MDM.filtered.faa',
                'L24': 'rpL24_JGI_MDM.filtered.faa',
                'L2': 'rpL2_JGI_MDM.filtered.faa',
                'L3': 'rpL3_JGI_MDM.filtered.faa',
                'L4': 'rpL4_JGI_MDM.filtered.faa',
                'L5': 'rpL5_JGI_MDM.filtered.faa',
                'L6': 'rpL6_JGI_MDM.filtered.faa',
                'S10': 'rpS10_JGI_MDM.filtered.faa',
                'S17': 'rpS17_JGI_MDM.filtered.faa',
                'S19': 'rpS19_JGI_MDM.filtered.faa',
                'S3': 'rpS3_JGI_MDM.filtered.faa',
                'S8': 'rpS8_JGI_MDM.filtered.faa'}
    protein_databases = {key: '%s/%s' % (databases, database) \
            for key, database in list(protein_databases.items())}
    return proteins, protein_databases
519499
def protein_ids(self, contig=None, strand=None):
        
        protein_ids = self._all_feature_values(
            column="protein_id",
            feature="CDS",
            contig=contig,
            strand=strand,
            distinct=True)
        # drop None values
        return [protein_id for protein_id in protein_ids if protein_id]
865140
def calcThermo(seq1, seq2, calc_type='ANY', mv_conc=50, dv_conc=0,
                 dntp_conc=0.8, dna_conc=50, temp_c=37, max_loop=30,
                 temp_only=False):
    
    args = [pjoin(PRIMER3_HOME, 'ntthal'),
            '-a',       str(calc_type),
            '-mv',      str(mv_conc),
            '-dv',      str(dv_conc),
            '-n',       str(dntp_conc),
            '-d',       str(dna_conc),
            '-t',       str(temp_c),
            '-maxloop', str(max_loop),
            '-path',    THERMO_PATH,
            '-s1',      seq1,
            '-s2',      seq2]
    if temp_only:
        args += ['-r']
    out = subprocess.check_output(args, stderr=DEV_NULL,
                                  env=os.environ)
    return _parse_ntthal(out)
580744
def parse_multi_object_delete_response(data):
    
    root = S3Element.fromstring('MultiObjectDeleteResult', data)
    return [
        MultiDeleteError(errtag.get_child_text('Key'),
                         errtag.get_child_text('Code'),
                         errtag.get_child_text('Message'))
        for errtag in root.findall('Error')
    ]
602806
def get_proteins_for_db(fastafn):
    
    objects = {}
    for record in parse_fasta(fastafn):
        objects[parse_protein_identifier(record)] = record
    return (((acc,) for acc in list(objects)),
            ((acc, str(record.seq)) for acc, record in objects.items()),
            ((acc, get_uniprot_evidence_level(record.description))
             for acc, record in objects.items()))
1129588
@Override
    public DescriptorValue calculate(IAtomContainer container) {
        if (taeParams == null) return getDummyDescriptorValue(new CDKException("TAE parameters were not initialized"));
        if (!(container instanceof IBioPolymer))
            return getDummyDescriptorValue(new CDKException("The molecule should be of type IBioPolymer"));

        IBioPolymer peptide = (IBioPolymer) container;

        // I assume that we get single letter names
        //Collection aas = peptide.getMonomerNames();

        double[] desc = new double[ndesc];
        for (int i = 0; i < ndesc; i++)
            desc[i] = 0.0;

        List<IMonomer> monomers = getMonomers(peptide);

        for (Iterator<IMonomer> iterator = monomers.iterator(); iterator.hasNext();) {
            IMonomer monomer = iterator.next();

            String o = monomer.getMonomerName();

            if (o.length() == 0) continue;

            String olc = String.valueOf(o.toLowerCase().charAt(0));
            String tlc = (String) nametrans.get(olc);

            logger.debug("Converted " + olc + " to " + tlc);

            // get the params for this AA
            Double[] params = (Double[]) taeParams.get(tlc);

            for (int i = 0; i < ndesc; i++)
                desc[i] += params[i];
        }

        DoubleArrayResult retval = new DoubleArrayResult(ndesc);
        for (int i = 0; i < ndesc; i++)
            retval.add(desc[i]);

        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                getDescriptorNames());
    }
509737
def blastp(seqs, blast_db="nr", e_value="1e-20", max_hits=200,
           working_dir=tempfile.gettempdir(), blast_mat_root=None,
           extra_params={}):
    

    # set up params to use with blastp
    params = {
        # matrix
        "-M":"BLOSUM62",

        # max procs
        "-a":"1",

        # expectation
        "-e":e_value,

        # max seqs to show
        "-b":max_hits,

        # max one line descriptions
        "-v":max_hits,

        # program
        "-p":"blastp"
    }
    params.update(extra_params)

    # blast
    blast_res =  blast_seqs(seqs,
        Blastall,
        blast_mat_root=blast_mat_root,
        blast_db=blast_db,
        params=params,
        add_seq_names=False,
        WorkingDir=working_dir
        )

    # get prot id map
    if blast_res['StdOut']:
        lines = [x for x in blast_res['StdOut']]
        return BlastResult(lines)

    return None
256456
def tag(ctx, corpus, output):
    """"""
    click.echo('chemdataextractor.chemdner.tag')
    for line in corpus:
        pmid, title, abstract = line.strip().split(u'\t')
        # print(pmid)
        counter = 1
        d = Document(Title(title), Paragraph(abstract))
        for t, section in [(d.elements[0], u'T'), (d.elements[1], u'A')]:
            for cem in t.cems:
                code = u'%s:%s:%s' % (section, cem.start, cem.end)
                output.write(u'\t'.join([pmid, code, six.text_type(counter), u'1']))
                output.write(u'\n')
                counter += 1
198873
def batch_taxonomy(list_of_taxids):
    
    for taxid in list_of_taxids:
        handle = Entrez.efetch(db='Taxonomy', id=taxid, retmode="xml")
        records = Entrez.read(handle)
        yield records[0]["ScientificName"]
1051162
def _extract_number_of_taxa(self):
        
        n_taxa = dict()
        for i in self.seq_records:
            if i.gene_code not in n_taxa:
                n_taxa[i.gene_code] = 0
            n_taxa[i.gene_code] += 1
        number_taxa = sorted([i for i in n_taxa.values()], reverse=True)[0]
        self.number_taxa = str(number_taxa)
749953
def parse_psimitab(content, fmt='tab27'):
    
    columns = [
        'Unique identifier for interactor A',
        'Unique identifier for interactor B',
        'Alternative identifier for interactor A',
        'Alternative identifier for interactor B',
        'Aliases for A',
        'Aliases for B',
        'Interaction detection methods',
        'First author',
        'Identifier of the publication',
        'NCBI Taxonomy identifier for interactor A',
        'NCBI Taxonomy identifier for interactor B',
        'Interaction types',
        'Source databases',
        'Interaction identifier(s)',
        'Confidence score']
    columns += [
        'Complex expansion',
        'Biological role A', 'Biological role B',
        'Experimental role A', 'Experimental role B',
        'Interactor type A', 'Interactor type B',
        'Xref for interactor A', 'Xref for interactor B',
        'Xref for the interaction',
        'Annotaions for interactor A', 'Annotations for interactor B',
        'Annotations for the interaction',
        'NCBI Taxonomy identifier for the host organism',
        'Prameters of the interaction',
        'Creaction date', 'Update date',
        'Checksum for the interactor A', 'Checksum for the interactor B',
        'Checksum for the interaction',
        'negative',
        'Feature(s) for interactor A', 'Feature(s) for interactor B',
        'Stoichiometry for interactor A', 'Stoichiometroy for interactor B',
        'Participant identification method for interactor A',
        'Participant identification method for interactor B'
        ]
    if fmt == 'tab25':
        columns = columns[: 15]

    rexp = re.compile(r"(?P<fields>((\"([^\"]|((?<=\\)\"))*\")|([^\t\"])|((?<=\\)\"))+)(\t|$)")

    retval = []
    for line in content.split('\n'):
        line = line.strip()
        if line == '' or line[0] == '#':
            continue

        start = 0
        tmp = []
        for mobj in rexp.finditer(line):
            if mobj.start() != start:
                print(repr(line))
            assert mobj.start() == start
            start = mobj.end()
            tmp.append(mobj.group('fields'))
        assert len(tmp) == len(columns)
        retval.append(dict(zip(columns, tmp)))
    return retval
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)