docs: def remove_artifacts_seqs(seqs_fp,
                            ref_fp,
                            working_dir,
                            ref_db_fp,
                            negate=False,
                            threads=1,
                            verbose=False,
                            sim_thresh=None,
                            coverage_thresh=None):
      """Remove artifacts from FASTA file using SortMeRNA.
  
      Parameters
      ----------
      seqs_fp: string
          file path to FASTA input sequence file
      ref_fp: tuple
          file path(s) to FASTA database file
      working_dir: string
          working directory path
      ref_db_fp: tuple
          file path(s) to indexed FASTA database
      negate: boolean, optional
          if True, discard all input sequences aligning
          to reference database
      threads: integer, optional
          number of threads to use for SortMeRNA
      verbose: boolean, optional
          If true, output SortMeRNA errors
      sim_thresh: float, optional
          The minimal similarity threshold (between 0 and 1)
          for keeping the sequence
          if None, the default values used are 0.65 for negate=False,
          0.95 for negate=True
      coverage_thresh: float, optional
          The minimal coverage threshold (between 0 and 1)
          for alignments for keeping the sequence
          if None, the default values used are 0.5 for negate=False,
          0.95 for negate=True
  
      Returns
      -------
      output_fp : str
          Name of the artifact removed fasta file
      okseqs : int
          The number of sequences left after artifact removal
      tmp_files : list of str
          Names of the tmp files created
      """
      logger = logging.getLogger(__name__)
      logger.info('remove_artifacts_seqs file %s' % seqs_fp)
  
      if stat(seqs_fp).st_size == 0:
          logger.warn('file %s has size 0, continuing' % seqs_fp)
          return None, 0, []
  
      if coverage_thresh is None:
          if negate:
              coverage_thresh = 0.95 * 100
          else:
              coverage_thresh = 0.5 * 100
  
      if sim_thresh is None:
          if negate:
              sim_thresh = 0.95 * 100
          else:
              sim_thresh = 0.65 * 100
  
      # the minimal average bitscore per nucleotide
      bitscore_thresh = 0.65
  
      output_fp = join(working_dir,
                       "%s.no_artifacts" % basename(seqs_fp))
      blast_output = join(working_dir,
                          '%s.sortmerna' % basename(seqs_fp))
      aligned_seq_ids = set()
      for i, db in enumerate(ref_fp):
          logger.debug('running on ref_fp %s working dir %s refdb_fp %s seqs %s'
                       % (db, working_dir, ref_db_fp[i], seqs_fp))
          # run SortMeRNA
          # we use -e 100 to remove E-value based filtering by sortmerna
          # since we use bitscore/identity/coverage filtering instead
          params = ['sortmerna', '--reads', seqs_fp, '--ref', '%s,%s' %
                    (db, ref_db_fp[i]),
                    '--aligned', blast_output, '--blast', '3', '--best', '1',
                    '--print_all_reads', '-v', '-e', '100']
  
          sout, serr, res = _system_call(params)
          if not res == 0:
              logger.error('sortmerna error on file %s' % seqs_fp)
              logger.error('stdout : %s' % sout)
              logger.error('stderr : %s' % serr)
              return output_fp, 0, []
  
          blast_output_filename = '%s.blast' % blast_output
          with open(blast_output_filename, 'r') as bfl:
              for line in bfl:
                  line = line.strip().split('\t')
                  # if * means no match
                  if line[1] == '*':
                      continue
                  # check if % identity[2] and coverage[13] are large enough
                  if (float(line[2]) >= sim_thresh) and \
                     (float(line[13]) >= coverage_thresh) and \
                     (float(line[11]) >= bitscore_thresh * len(line[0])):
                      aligned_seq_ids.add(line[0])
  
      if negate:
          def op(x): return x not in aligned_seq_ids
      else:
          def op(x): return x in aligned_seq_ids
  
      # if negate = False, only output sequences
      # matching to at least one of the databases
      totalseqs = 0
      okseqs = 0
      badseqs = 0
      with open(output_fp, 'w') as out_f:
          for label, seq in sequence_generator(seqs_fp):
              totalseqs += 1
              label = label.split()[0]
              if op(label):
                  out_f.write(">%s\n%s\n" % (label, seq))
                  okseqs += 1
              else:
                  badseqs += 1
      logger.info('total sequences %d, passing sequences %d, '
                  'failing sequences %d' % (totalseqs, okseqs, badseqs))
      return output_fp, okseqs, [blast_output_filename]