code: 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):
      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
      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))
          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 line[1] == '*':
                      continue
                  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
      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]