nodes: 0, module; 1, function_definition; 2, function_name:remove_artifacts_seqs; 3, parameters; 4, block; 5, identifier:seqs_fp; 6, identifier:ref_fp; 7, identifier:working_dir; 8, identifier:ref_db_fp; 9, default_parameter; 10, default_parameter; 11, default_parameter; 12, default_parameter; 13, default_parameter; 14, expression_statement; 15, expression_statement; 16, expression_statement; 17, if_statement; 18, if_statement; 19, if_statement; 20, comment:# the minimal average bitscore per nucleotide; 21, expression_statement; 22, expression_statement; 23, expression_statement; 24, expression_statement; 25, for_statement; 26, if_statement; 27, comment:# if negate = False, only output sequences; 28, comment:# matching to at least one of the databases; 29, expression_statement; 30, expression_statement; 31, expression_statement; 32, with_statement; 33, expression_statement; 34, return_statement; 35, identifier:negate; 36, False; 37, identifier:threads; 38, integer:1; 39, identifier:verbose; 40, False; 41, identifier:sim_thresh; 42, None; 43, identifier:coverage_thresh; 44, None; 45, comment:"""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
      """; 46, assignment; 47, call; 48, comparison_operator:stat(seqs_fp).st_size == 0; 49, block; 50, comparison_operator:coverage_thresh is None; 51, block; 52, comparison_operator:sim_thresh is None; 53, block; 54, assignment; 55, assignment; 56, assignment; 57, assignment; 58, pattern_list; 59, call; 60, block; 61, identifier:negate; 62, block; 63, else_clause; 64, assignment; 65, assignment; 66, assignment; 67, with_clause; 68, block; 69, call; 70, expression_list; 71, identifier:logger; 72, call; 73, attribute; 74, argument_list; 75, attribute; 76, integer:0; 77, expression_statement; 78, return_statement; 79, identifier:coverage_thresh; 80, None; 81, if_statement; 82, identifier:sim_thresh; 83, None; 84, if_statement; 85, identifier:bitscore_thresh; 86, float:0.65; 87, identifier:output_fp; 88, call; 89, identifier:blast_output; 90, call; 91, identifier:aligned_seq_ids; 92, call; 93, identifier:i; 94, identifier:db; 95, identifier:enumerate; 96, argument_list; 97, expression_statement; 98, comment:# run SortMeRNA; 99, comment:# we use -e 100 to remove E-value based filtering by sortmerna; 100, comment:# since we use bitscore/identity/coverage filtering instead; 101, expression_statement; 102, expression_statement; 103, if_statement; 104, expression_statement; 105, with_statement; 106, function_definition; 107, block; 108, identifier:totalseqs; 109, integer:0; 110, identifier:okseqs; 111, integer:0; 112, identifier:badseqs; 113, integer:0; 114, with_item; 115, for_statement; 116, attribute; 117, argument_list; 118, identifier:output_fp; 119, identifier:okseqs; 120, list; 121, attribute; 122, argument_list; 123, identifier:logger; 124, identifier:info; 125, binary_operator:'remove_artifacts_seqs file %s' % seqs_fp; 126, call; 127, identifier:st_size; 128, call; 129, expression_list; 130, identifier:negate; 131, block; 132, else_clause; 133, identifier:negate; 134, block; 135, else_clause; 136, identifier:join; 137, argument_list; 138, identifier:join; 139, argument_list; 140, identifier:set; 141, argument_list; 142, identifier:ref_fp; 143, call; 144, assignment; 145, assignment; 146, not_operator; 147, block; 148, assignment; 149, with_clause; 150, block; 151, function_name:op; 152, parameters; 153, block; 154, function_definition; 155, as_pattern; 156, pattern_list; 157, call; 158, block; 159, identifier:logger; 160, identifier:info; 161, binary_operator:'total sequences %d, passing sequences %d, '
                  'failing sequences %d' % (totalseqs, okseqs, badseqs); 162, identifier:blast_output_filename; 163, identifier:logging; 164, identifier:getLogger; 165, identifier:__name__; 166, string; 167, identifier:seqs_fp; 168, identifier:stat; 169, argument_list; 170, attribute; 171, argument_list; 172, None; 173, integer:0; 174, list; 175, expression_statement; 176, block; 177, expression_statement; 178, block; 179, identifier:working_dir; 180, binary_operator:"%s.no_artifacts" % basename(seqs_fp); 181, identifier:working_dir; 182, binary_operator:'%s.sortmerna' % basename(seqs_fp); 183, attribute; 184, argument_list; 185, identifier:params; 186, list; 187, pattern_list; 188, call; 189, comparison_operator:res == 0; 190, expression_statement; 191, expression_statement; 192, expression_statement; 193, return_statement; 194, identifier:blast_output_filename; 195, binary_operator:'%s.blast' % blast_output; 196, with_item; 197, for_statement; 198, identifier:x; 199, return_statement; 200, function_name:op; 201, parameters; 202, block; 203, call; 204, as_pattern_target; 205, identifier:label; 206, identifier:seq; 207, identifier:sequence_generator; 208, argument_list; 209, expression_statement; 210, expression_statement; 211, if_statement; 212, concatenated_string; 213, tuple; 214, string_content:remove_artifacts_seqs file %s; 215, identifier:seqs_fp; 216, identifier:logger; 217, identifier:warn; 218, binary_operator:'file %s has size 0, continuing' % seqs_fp; 219, assignment; 220, expression_statement; 221, assignment; 222, expression_statement; 223, string:"%s.no_artifacts"; 224, call; 225, string; 226, call; 227, identifier:logger; 228, identifier:debug; 229, binary_operator:'running on ref_fp %s working dir %s refdb_fp %s seqs %s'
                       % (db, working_dir, ref_db_fp[i], seqs_fp); 230, string; 231, string; 232, identifier:seqs_fp; 233, string; 234, binary_operator:'%s,%s' %
                    (db, ref_db_fp[i]); 235, string; 236, identifier:blast_output; 237, string; 238, string; 239, string; 240, string; 241, string; 242, string; 243, string; 244, string; 245, identifier:sout; 246, identifier:serr; 247, identifier:res; 248, identifier:_system_call; 249, argument_list; 250, identifier:res; 251, integer:0; 252, call; 253, call; 254, call; 255, expression_list; 256, string; 257, identifier:blast_output; 258, as_pattern; 259, identifier:line; 260, identifier:bfl; 261, block; 262, comparison_operator:x not in aligned_seq_ids; 263, identifier:x; 264, return_statement; 265, identifier:open; 266, argument_list; 267, identifier:out_f; 268, identifier:seqs_fp; 269, augmented_assignment; 270, assignment; 271, call; 272, block; 273, else_clause; 274, string; 275, string; 276, identifier:totalseqs; 277, identifier:okseqs; 278, identifier:badseqs; 279, string; 280, identifier:seqs_fp; 281, identifier:coverage_thresh; 282, binary_operator:0.95 * 100; 283, assignment; 284, identifier:sim_thresh; 285, binary_operator:0.95 * 100; 286, assignment; 287, identifier:basename; 288, argument_list; 289, string_content:%s.sortmerna; 290, identifier:basename; 291, argument_list; 292, string; 293, tuple; 294, string_content:sortmerna; 295, string_content:--reads; 296, string_content:--ref; 297, string; 298, tuple; 299, string_content:--aligned; 300, string_content:--blast; 301, string_content:3; 302, string_content:--best; 303, string_content:1; 304, string_content:--print_all_reads; 305, string_content:-v; 306, string_content:-e; 307, string_content:100; 308, identifier:params; 309, attribute; 310, argument_list; 311, attribute; 312, argument_list; 313, attribute; 314, argument_list; 315, identifier:output_fp; 316, integer:0; 317, list; 318, string_content:%s.blast; 319, call; 320, as_pattern_target; 321, expression_statement; 322, comment:# if * means no match; 323, if_statement; 324, comment:# check if % identity[2] and coverage[13] are large enough; 325, if_statement; 326, identifier:x; 327, identifier:aligned_seq_ids; 328, comparison_operator:x in aligned_seq_ids; 329, identifier:output_fp; 330, string; 331, identifier:totalseqs; 332, integer:1; 333, identifier:label; 334, subscript; 335, identifier:op; 336, argument_list; 337, expression_statement; 338, expression_statement; 339, block; 340, string_content:total sequences %d, passing sequences %d,; 341, string_content:failing sequences %d; 342, string_content:file %s has size 0, continuing; 343, float:0.95; 344, integer:100; 345, identifier:coverage_thresh; 346, binary_operator:0.5 * 100; 347, float:0.95; 348, integer:100; 349, identifier:sim_thresh; 350, binary_operator:0.65 * 100; 351, identifier:seqs_fp; 352, identifier:seqs_fp; 353, string_content:running on ref_fp %s working dir %s refdb_fp %s seqs %s; 354, identifier:db; 355, identifier:working_dir; 356, subscript; 357, identifier:seqs_fp; 358, string_content:%s,%s; 359, identifier:db; 360, subscript; 361, identifier:logger; 362, identifier:error; 363, binary_operator:'sortmerna error on file %s' % seqs_fp; 364, identifier:logger; 365, identifier:error; 366, binary_operator:'stdout : %s' % sout; 367, identifier:logger; 368, identifier:error; 369, binary_operator:'stderr : %s' % serr; 370, identifier:open; 371, argument_list; 372, identifier:bfl; 373, assignment; 374, comparison_operator:line[1] == '*'; 375, block; 376, boolean_operator; 377, block; 378, identifier:x; 379, identifier:aligned_seq_ids; 380, string_content:w; 381, call; 382, integer:0; 383, identifier:label; 384, call; 385, augmented_assignment; 386, expression_statement; 387, float:0.5; 388, integer:100; 389, float:0.65; 390, integer:100; 391, identifier:ref_db_fp; 392, identifier:i; 393, identifier:ref_db_fp; 394, identifier:i; 395, string; 396, identifier:seqs_fp; 397, string; 398, identifier:sout; 399, string; 400, identifier:serr; 401, identifier:blast_output_filename; 402, string; 403, identifier:line; 404, call; 405, subscript; 406, string; 407, continue_statement; 408, boolean_operator; 409, line_continuation:\; 410, parenthesized_expression; 411, expression_statement; 412, attribute; 413, argument_list; 414, attribute; 415, argument_list; 416, identifier:okseqs; 417, integer:1; 418, augmented_assignment; 419, string_content:sortmerna error on file %s; 420, string_content:stdout : %s; 421, string_content:stderr : %s; 422, string_content:r; 423, attribute; 424, argument_list; 425, identifier:line; 426, integer:1; 427, string_content:*; 428, parenthesized_expression; 429, line_continuation:\; 430, parenthesized_expression; 431, comparison_operator:float(line[11]) >= bitscore_thresh * len(line[0]); 432, call; 433, identifier:label; 434, identifier:split; 435, identifier:out_f; 436, identifier:write; 437, binary_operator:">%s\n%s\n" % (label, seq); 438, identifier:badseqs; 439, integer:1; 440, call; 441, identifier:split; 442, string; 443, comparison_operator:float(line[2]) >= sim_thresh; 444, comparison_operator:float(line[13]) >= coverage_thresh; 445, call; 446, binary_operator:bitscore_thresh * len(line[0]); 447, attribute; 448, argument_list; 449, string:">%s\n%s\n"; 450, tuple; 451, attribute; 452, argument_list; 453, string_content; 454, call; 455, identifier:sim_thresh; 456, call; 457, identifier:coverage_thresh; 458, identifier:float; 459, argument_list; 460, identifier:bitscore_thresh; 461, call; 462, identifier:aligned_seq_ids; 463, identifier:add; 464, subscript; 465, identifier:label; 466, identifier:seq; 467, identifier:line; 468, identifier:strip; 469, escape_sequence:\t; 470, identifier:float; 471, argument_list; 472, identifier:float; 473, argument_list; 474, subscript; 475, identifier:len; 476, argument_list; 477, identifier:line; 478, integer:0; 479, subscript; 480, subscript; 481, identifier:line; 482, integer:11; 483, subscript; 484, identifier:line; 485, integer:2; 486, identifier:line; 487, integer:13; 488, identifier:line; 489, integer:0
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):
      """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]