nodes: 0, module; 1, function_definition; 2, function_name:get_region_nt_counts; 3, parameters; 4, block; 5, identifier:region; 6, identifier:bam; 7, default_parameter; 8, expression_statement; 9, comment:# TODO: I should figure out what the different possible values are that; 10, comment:# pysam could give me back (so far I only have ATCGN). Can I get deletions; 11, comment:# and insertions?; 12, comment:# TODO: This could probably be parallelized.; 13, if_statement; 14, if_statement; 15, expression_statement; 16, if_statement; 17, expression_statement; 18, for_statement; 19, return_statement; 20, identifier:stranded; 21, False; 22, comment:"""
      Get counts of each nucleotide from a bam file for a given region. If R1 and
      R2 reads both overlap a position, only one count will be added. If the R1
      and R2 reads disagree at a position they both overlap, that read pair is not
      used for that position.  Can optionally output strand-specific counts.
  
      Parameters
      ----------
      region : str or list
          Region of type chrom:start-end, chrom:start-end:strand, or [chrom,
          start, end]. The strand is ignored for chrom:start-end:strand. For
          chrom:start-end, the coordinates are one-based inclusive. For example,
          the query chr1:10-11 will give you the counts for the 10th and 11th
          bases of chr1. For [chrom, start, end], the coordinates are zero-based
          and end exclusive (like a bed file). The query [chr1, 9, 11] will give
          you the coverage of the 10th and 11th bases of chr1. The region value is
          passed directly to pysam's pileup function. 
  
      bam : pysam.calignmentfile.AlignmentFile or str
          Bam file opened with pysam or path to bam file (must be sorted and
          indexed).
  
      stranded : boolean
          Boolean indicating whether read data is stranded and stranded nucleotide
          counts should be returned. Assumes R1 read on reverse strand implies +
          strand coverage etc.
  
      Returns
      -------
      counts : pandas.DataFrame
          Data frame with the counts for each base in the region. The index of
          this data frame is one-based for compatibility with VCF files.
  
      """; 23, comparison_operator:type(bam) == str; 24, block; 25, comparison_operator:type(region) is str; 26, block; 27, elif_clause; 28, assignment; 29, identifier:stranded; 30, block; 31, assignment; 32, identifier:pc; 33, identifier:pp; 34, comment:# Most of this code deals with R1 and R2 reads that overlap so that we; 35, comment:# don't get two counts from one fragment.; 36, block; 37, identifier:counts; 38, call; 39, identifier:str; 40, expression_statement; 41, call; 42, identifier:str; 43, expression_statement; 44, if_statement; 45, expression_statement; 46, expression_statement; 47, expression_statement; 48, expression_statement; 49, comparison_operator:type(region) is (list or tuple); 50, block; 51, identifier:cols; 52, list; 53, expression_statement; 54, identifier:counts; 55, call; 56, expression_statement; 57, expression_statement; 58, expression_statement; 59, expression_statement; 60, expression_statement; 61, for_statement; 62, expression_statement; 63, expression_statement; 64, expression_statement; 65, expression_statement; 66, expression_statement; 67, expression_statement; 68, expression_statement; 69, expression_statement; 70, expression_statement; 71, expression_statement; 72, for_statement; 73, identifier:type; 74, argument_list; 75, assignment; 76, identifier:type; 77, argument_list; 78, assignment; 79, comparison_operator:len(r) == 3; 80, block; 81, elif_clause; 82, assignment; 83, assignment; 84, assignment; 85, assignment; 86, call; 87, parenthesized_expression; 88, expression_statement; 89, expression_statement; 90, expression_statement; 91, string; 92, string; 93, string; 94, string; 95, string; 96, assignment; 97, attribute; 98, argument_list; 99, assignment; 100, assignment; 101, assignment; 102, assignment; 103, assignment; 104, identifier:pr; 105, attribute; 106, block; 107, assignment; 108, assignment; 109, assignment; 110, assignment; 111, assignment; 112, assignment; 113, call; 114, call; 115, assignment; 116, call; 117, identifier:vc; 118, identifier:vcs; 119, block; 120, identifier:bam; 121, identifier:bam; 122, call; 123, identifier:region; 124, identifier:r; 125, call; 126, call; 127, integer:3; 128, expression_statement; 129, comparison_operator:len(r) == 4; 130, block; 131, identifier:start; 132, call; 133, identifier:end; 134, call; 135, identifier:ind; 136, list_comprehension; 137, identifier:pp; 138, call; 139, identifier:type; 140, argument_list; 141, boolean_operator; 142, assignment; 143, assignment; 144, assignment; 145, string_content:A; 146, string_content:T; 147, string_content:C; 148, string_content:G; 149, string_content:N; 150, identifier:cols; 151, binary_operator:['{}+'.format(x) for x in cols] + ['{}-'.format(x) for x in cols]; 152, identifier:pd; 153, identifier:DataFrame; 154, integer:0; 155, keyword_argument; 156, keyword_argument; 157, identifier:pos; 158, binary_operator:pc.reference_pos + 1; 159, identifier:r1_qnames; 160, list; 161, identifier:r1_nts; 162, list; 163, identifier:r2_qnames; 164, list; 165, identifier:r2_nts; 166, list; 167, identifier:pc; 168, identifier:pileups; 169, expression_statement; 170, expression_statement; 171, expression_statement; 172, if_statement; 173, identifier:r1; 174, call; 175, identifier:r2; 176, call; 177, identifier:df; 178, attribute; 179, identifier:singles; 180, subscript; 181, identifier:doubles; 182, call; 183, identifier:vcs; 184, list; 185, attribute; 186, argument_list; 187, attribute; 188, argument_list; 189, identifier:doubles; 190, subscript; 191, attribute; 192, argument_list; 193, expression_statement; 194, attribute; 195, argument_list; 196, identifier:parse_region; 197, argument_list; 198, identifier:len; 199, argument_list; 200, assignment; 201, call; 202, integer:4; 203, expression_statement; 204, identifier:int; 205, argument_list; 206, identifier:int; 207, argument_list; 208, call; 209, for_in_clause; 210, attribute; 211, argument_list; 212, identifier:region; 213, identifier:list; 214, identifier:tuple; 215, pattern_list; 216, identifier:region; 217, identifier:ind; 218, list_comprehension; 219, identifier:pp; 220, call; 221, list_comprehension; 222, list_comprehension; 223, identifier:index; 224, identifier:ind; 225, identifier:columns; 226, identifier:cols; 227, attribute; 228, integer:1; 229, assignment; 230, assignment; 231, assignment; 232, identifier:nt; 233, block; 234, attribute; 235, argument_list; 236, attribute; 237, argument_list; 238, call; 239, identifier:T; 240, identifier:df; 241, comparison_operator:df.isnull().sum(axis=1) == 1; 242, attribute; 243, argument_list; 244, identifier:vcs; 245, identifier:append; 246, call; 247, identifier:vcs; 248, identifier:append; 249, call; 250, identifier:doubles; 251, comparison_operator:doubles.R1 == doubles.R2; 252, identifier:vcs; 253, identifier:append; 254, call; 255, augmented_assignment; 256, identifier:pysam; 257, identifier:AlignmentFile; 258, identifier:bam; 259, string; 260, identifier:region; 261, identifier:r; 262, pattern_list; 263, identifier:r; 264, identifier:len; 265, argument_list; 266, assignment; 267, identifier:start; 268, identifier:end; 269, attribute; 270, argument_list; 271, identifier:x; 272, call; 273, identifier:bam; 274, identifier:pileup; 275, keyword_argument; 276, keyword_argument; 277, identifier:chrom; 278, identifier:start; 279, identifier:end; 280, call; 281, for_in_clause; 282, attribute; 283, argument_list; 284, call; 285, for_in_clause; 286, call; 287, for_in_clause; 288, identifier:pc; 289, identifier:reference_pos; 290, identifier:qnames; 291, subscript; 292, identifier:nts; 293, subscript; 294, identifier:nt; 295, call; 296, expression_statement; 297, expression_statement; 298, identifier:pd; 299, identifier:Series; 300, identifier:r1_nts; 301, keyword_argument; 302, identifier:pd; 303, identifier:Series; 304, identifier:r2_nts; 305, keyword_argument; 306, attribute; 307, argument_list; 308, call; 309, integer:1; 310, identifier:df; 311, identifier:dropna; 312, attribute; 313, argument_list; 314, attribute; 315, argument_list; 316, attribute; 317, attribute; 318, attribute; 319, argument_list; 320, subscript; 321, identifier:vc; 322, string_content:rb; 323, identifier:chrom; 324, identifier:start; 325, identifier:end; 326, identifier:r; 327, pattern_list; 328, identifier:r; 329, string; 330, identifier:format; 331, identifier:chrom; 332, identifier:x; 333, identifier:range; 334, argument_list; 335, identifier:region; 336, identifier:region; 337, identifier:truncate; 338, True; 339, attribute; 340, argument_list; 341, identifier:x; 342, call; 343, identifier:bam; 344, identifier:pileup; 345, identifier:chrom; 346, identifier:start; 347, identifier:end; 348, keyword_argument; 349, attribute; 350, argument_list; 351, identifier:x; 352, identifier:cols; 353, attribute; 354, argument_list; 355, identifier:x; 356, identifier:cols; 357, list; 358, attribute; 359, list; 360, attribute; 361, identifier:_pos_nt; 362, argument_list; 363, call; 364, call; 365, identifier:index; 366, identifier:r1_qnames; 367, identifier:index; 368, identifier:r2_qnames; 369, identifier:pd; 370, identifier:DataFrame; 371, list; 372, keyword_argument; 373, attribute; 374, argument_list; 375, subscript; 376, identifier:value_counts; 377, subscript; 378, identifier:value_counts; 379, identifier:doubles; 380, identifier:R1; 381, identifier:doubles; 382, identifier:R2; 383, attribute; 384, identifier:value_counts; 385, attribute; 386, call; 387, attribute; 388, identifier:chrom; 389, identifier:start; 390, identifier:end; 391, identifier:strand; 392, string_content:{}:{}; 393, identifier:start; 394, binary_operator:end + 1; 395, string; 396, identifier:format; 397, identifier:chrom; 398, identifier:x; 399, identifier:range; 400, argument_list; 401, identifier:truncate; 402, True; 403, string; 404, identifier:format; 405, identifier:x; 406, string; 407, identifier:format; 408, identifier:x; 409, identifier:r1_qnames; 410, identifier:r2_qnames; 411, attribute; 412, identifier:is_read2; 413, identifier:r1_nts; 414, identifier:r2_nts; 415, attribute; 416, identifier:is_read2; 417, identifier:pr; 418, attribute; 419, identifier:stranded; 420, attribute; 421, argument_list; 422, attribute; 423, argument_list; 424, identifier:r1; 425, identifier:r2; 426, identifier:index; 427, list; 428, call; 429, identifier:sum; 430, keyword_argument; 431, identifier:singles; 432, string; 433, identifier:singles; 434, string; 435, identifier:doubles; 436, identifier:R1; 437, identifier:counts; 438, identifier:ix; 439, attribute; 440, argument_list; 441, identifier:vc; 442, identifier:index; 443, identifier:end; 444, integer:1; 445, string_content:{}:{}; 446, binary_operator:int(start) + 1; 447, binary_operator:int(end) + 1; 448, string_content:{}+; 449, string_content:{}-; 450, identifier:pr; 451, identifier:alignment; 452, identifier:pr; 453, identifier:alignment; 454, identifier:pc; 455, identifier:reference_pos; 456, identifier:qnames; 457, identifier:append; 458, attribute; 459, identifier:nts; 460, identifier:append; 461, identifier:nt; 462, string; 463, string; 464, attribute; 465, argument_list; 466, identifier:axis; 467, integer:1; 468, string_content:R1; 469, string_content:R2; 470, string; 471, identifier:format; 472, identifier:chrom; 473, identifier:pos; 474, call; 475, integer:1; 476, call; 477, integer:1; 478, attribute; 479, identifier:qname; 480, string_content:R1; 481, string_content:R2; 482, identifier:df; 483, identifier:isnull; 484, string_content:{}:{}; 485, identifier:int; 486, argument_list; 487, identifier:int; 488, argument_list; 489, identifier:pr; 490, identifier:alignment; 491, identifier:start; 492, identifier:end
code: def get_region_nt_counts(region, bam, stranded=False):
      """
      Get counts of each nucleotide from a bam file for a given region. If R1 and
      R2 reads both overlap a position, only one count will be added. If the R1
      and R2 reads disagree at a position they both overlap, that read pair is not
      used for that position.  Can optionally output strand-specific counts.
  
      Parameters
      ----------
      region : str or list
          Region of type chrom:start-end, chrom:start-end:strand, or [chrom,
          start, end]. The strand is ignored for chrom:start-end:strand. For
          chrom:start-end, the coordinates are one-based inclusive. For example,
          the query chr1:10-11 will give you the counts for the 10th and 11th
          bases of chr1. For [chrom, start, end], the coordinates are zero-based
          and end exclusive (like a bed file). The query [chr1, 9, 11] will give
          you the coverage of the 10th and 11th bases of chr1. The region value is
          passed directly to pysam's pileup function. 
  
      bam : pysam.calignmentfile.AlignmentFile or str
          Bam file opened with pysam or path to bam file (must be sorted and
          indexed).
  
      stranded : boolean
          Boolean indicating whether read data is stranded and stranded nucleotide
          counts should be returned. Assumes R1 read on reverse strand implies +
          strand coverage etc.
  
      Returns
      -------
      counts : pandas.DataFrame
          Data frame with the counts for each base in the region. The index of
          this data frame is one-based for compatibility with VCF files.
  
      """
      # TODO: I should figure out what the different possible values are that
      # pysam could give me back (so far I only have ATCGN). Can I get deletions
      # and insertions? 
      # TODO: This could probably be parallelized.
      if type(bam) == str:
          bam = pysam.AlignmentFile(bam, 'rb')
  
      if type(region) is str:
          r = parse_region(region)
          if len(r) == 3:
              chrom, start, end = r
          elif len(r) == 4:
              chrom, start, end, strand = r
          start = int(start)
          end = int(end)
          ind = ['{}:{}'.format(chrom, x) for 
                 x in range(start, end + 1)]
          pp = bam.pileup(region=region, truncate=True)
      
      elif type(region) is (list or tuple):
          chrom, start, end = region
          ind = ['{}:{}'.format(chrom, x) for 
                 x in range(int(start) + 1, int(end) + 1)]
          pp = bam.pileup(chrom, start, end, truncate=True)
     
      cols = ['A', 'T', 'C', 'G', 'N']
      if stranded:
          cols = ['{}+'.format(x) for x in cols] + ['{}-'.format(x) for x in cols]
      counts = pd.DataFrame(0, index=ind, columns=cols)
      
      for pc in pp: 
          # Most of this code deals with R1 and R2 reads that overlap so that we
          # don't get two counts from one fragment.
          pos = pc.reference_pos + 1
          r1_qnames = []
          r1_nts = []
          r2_qnames = []
          r2_nts = []
          for pr in pc.pileups:
              qnames = [r1_qnames, r2_qnames][pr.alignment.is_read2]
              nts = [r1_nts, r2_nts][pr.alignment.is_read2]
              nt = _pos_nt(pr, pc.reference_pos, stranded)
              if nt:
                  qnames.append(pr.alignment.qname)
                  nts.append(nt)
          r1 = pd.Series(r1_nts, index=r1_qnames)
          r2 = pd.Series(r2_nts, index=r2_qnames)
          df = pd.DataFrame([r1, r2], index=['R1', 'R2']).T
          singles = df[df.isnull().sum(axis=1) == 1]
          doubles = df.dropna()
          vcs = []
          vcs.append(singles['R1'].value_counts())
          vcs.append(singles['R2'].value_counts())
          doubles = doubles[doubles.R1 == doubles.R2]
          vcs.append(doubles.R1.value_counts())
          for vc in vcs:
              counts.ix['{}:{}'.format(chrom, pos), vc.index] += vc
      return counts