docs: def calculate_residue_counts_perstrain(protein_pickle_path, outdir, pdbflex_keys_file, wt_pid_cutoff=None, force_rerun=False):
      """Writes out a feather file for a PROTEIN counting amino acid occurences for ALL STRAINS along with SUBSEQUENCES"""
      from collections import defaultdict
      from ssbio.protein.sequence.seqprop import SeqProp
      from ssbio.protein.sequence.properties.residues import _aa_property_dict_one
      log = logging.getLogger(__name__)
  
      protein_id = op.splitext(op.basename(protein_pickle_path))[0].split('_')[0]
  
      protein_df_outfile = op.join(outdir, '{}_protein_strain_properties.fthr'.format(protein_id))
  
      if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_df_outfile):
  
          protein = ssbio.io.load_pickle(protein_pickle_path)
  
          # First calculate disorder cuz i forgot to
          protein.get_all_disorder_predictions(representative_only=True)
  
          # Then get all subsequences
          all_protein_subseqs = protein.get_all_subsequences(pdbflex_keys_file=pdbflex_keys_file)
          if not all_protein_subseqs:
              log.error('{}: cannot run subsequence calculator'.format(protein.id))
              return
  
          # Each strain gets a dictionary
          strain_to_infodict = defaultdict(dict)
  
          for seqprop_to_analyze in protein.sequences:
              if seqprop_to_analyze.id == protein.representative_sequence.id:
                  strain_id = 'K12'
              elif type(seqprop_to_analyze) == SeqProp and seqprop_to_analyze.id != protein.id:  # This is to filter out other KEGGProps or UniProtProps
                  strain_id = seqprop_to_analyze.id.split('_', 1)[1]  # This split should work for all strains
              else:
                  continue
  
              ## Additional filtering for genes marked as orthologous but actually have large deletions or something
              ## TODO: experiment with other cutoffs?
              if wt_pid_cutoff:
                  aln = protein.sequence_alignments.get_by_id('{0}_{0}_{1}'.format(protein.id, seqprop_to_analyze.id))
                  if aln.annotations['percent_identity'] < wt_pid_cutoff:
                      continue
  
              ###### Calculate "all" properties ######
              seqprop_to_analyze.get_biopython_pepstats()
  
              # [ALL] aa_count
              if 'amino_acids_percent-biop' not in seqprop_to_analyze.annotations:  # May not run if weird amino acids in the sequence
                  log.warning('Protein {}, sequence {}: skipping, unable to run Biopython ProteinAnalysis'.format(protein.id,
                                                                                                              seqprop_to_analyze.id))
                  continue
              strain_to_infodict[strain_id].update({'aa_count_{}'.format(k): v for k, v in seqprop_to_analyze.annotations['amino_acids_content-biop'].items()})
  
              # [ALL] aa_count_total
              strain_to_infodict[strain_id]['aa_count_total'] = seqprop_to_analyze.seq_len
  
              ###### Calculate subsequence properties ######
              for prop, propdict in all_protein_subseqs.items():
                  strain_to_infodict[strain_id].update(protein.get_subseq_props(property_dict=propdict, property_name=prop,
                                                                                seqprop=seqprop_to_analyze))
  
          protein_df = pd.DataFrame(strain_to_infodict)
          protein_df.reset_index().to_feather(protein_df_outfile)
  
      return protein_pickle_path, protein_df_outfile