lib/ims_indel/indel_detector.rb
def count_depth(all_pos, bam_path, target_chr, mapq)
      pos_max = all_pos.max 
      pos_min = all_pos.min
      all_pos_depth = Hash.new(0)

      SAMReader.new(@samtools, bam_path).each(chr: target_chr,
                                         filter_flag: SAMReader::DUPLICATE,
                                         output_flag: SAMReader::PROPER_PAIR) do |line|
        read_seq = line.read_seq
        next if target_chr != line.chr 
        next if line.map_score <= mapq # Condition1
        next if read_seq.count("N") > 0
        map_status = line.map_status

        # complete match
        if /^(\d+)M$/ =~ map_status # 161M
          match_size = $1.to_i  #  obj_SAM.map_status --> 70M or 101M or

          # B fragment (match + clip seq)
        elsif /^(\d+)M(?:\d+)S$/ =~ map_status # 80M20S
          len_M = $1.to_i
          match_size = len_M
          # B fragment (short INS + clip seq)
        elsif /^(\d+)M(?:\d+)I(\d+)M(?:\d+)S$/ =~ map_status # 70M2I8M20S
          len_M1 = $1.to_i #70
          len_M2 = $2.to_i #8
          match_size = len_M1 + len_M2
          # B fragment (short DEL + clip seq)
        elsif /^(\d+)M(\d+)D(\d+)M(?:\d+)S$/ =~ map_status # 70M2D10M20S
          len_M1 = $1.to_i #70
          len_D  = $2.to_i #2
          len_M2 = $3.to_i #8
          match_size = len_M1 + len_D + len_M2

          # F fragment (clip seq + match)
        elsif /^(?:\d+)S(\d+)M$/ =~ map_status # 20S80M
          len_M = $1.to_i
          match_size = len_M
          # F fragment (clip seq + short INS)
        elsif /^(?:\d+)S(\d+)M(?:\d+)I(\d+)M$/ =~ map_status # 20S70M2I8M
          len_M1 = $1.to_i #70
          len_M2 = $2.to_i #8
          match_size = len_M1 + len_M2
          # F fragment (clip seq + short DEL)
        elsif /^(?:\d+)S(\d+)M(\d+)D(\d+)M$/ =~ map_status # 20S70M2D10M
          len_M1 = $1.to_i #70
          len_D  = $2.to_i #2
          len_M2 = $3.to_i #8
          match_size = len_M1 + len_D + len_M2

          # -----------------------------------------
          # short INDEL
          # -----------------------------------------
          # 1 insertion 
        elsif /^(\d+)M(\d+)I(\d+)M$/ =~ map_status # 70M2I28M
          len_M1 = $1.to_i #70
          #len_I  = $2.to_i #2
          len_M2 = $3.to_i #28
          match_size = len_M1 + len_M2
          # 1 deletion
        elsif /^(\d+)M(\d+)D(\d+)M$/ =~ map_status # 70M5D30M
          len_M1 = $1.to_i #70
          len_D  = $2.to_i #5
          len_M2 = $3.to_i #30
          match_size = len_M1 + len_D + len_M2

          # 1 del + 1 del
        elsif /^(\d+)M(\d+)D(\d+)M(\d+)D(\d+)M$/ =~ map_status # 70M 2D 10M 5D 28M
          len_M1 = $1.to_i #70
          len_D1 = $2.to_i #2
          len_M2 = $3.to_i #10
          len_D2 = $4.to_i #5
          len_M3 = $5.to_i #28
          match_size  = len_M1 + len_D1 + len_M2 + len_D2 + len_M3
          # 1 del + 1 ins
        elsif /^(\d+)M(\d+)D(\d+)M(?:\d+)I(\d+)M$/ =~ map_status # 70M 2D 10M 5I 28M
          len_M1 = $1.to_i #70
          len_D1 = $2.to_i #2
          len_M2 = $3.to_i #10
          len_M3 = $4.to_i #28
          match_size  = len_M1 + len_D1 + len_M2 + len_M3
          # 1 ins + 1 del
        elsif /^(\d+)M(?:\d+)I(\d+)M(\d+)D(\d+)M$/ =~ map_status # 70M 2I 10M 5D 28M
          len_M1 = $1.to_i #70
          len_M2 = $2.to_i #10
          len_D1 = $3.to_i #5
          len_M3 = $4.to_i #28
          match_size  = len_M1 + len_M2 + len_D1 + len_M3
          # 1 ins + 1 ins
        elsif /^(\d+)M(?:\d+)I(\d+)M(?:\d+)I(\d+)M$/ =~ map_status # 70M 2I 10M 5I 28M
          len_M1 = $1.to_i #70
          len_M2 = $2.to_i #10
          len_M3 = $3.to_i #28
          match_size  = len_M1 + len_M2 + len_M3
        else
          next
        end

        sttpos = line.chrpos
        endpos = sttpos + match_size - 1
        next if endpos < pos_min
        break if pos_max < sttpos

        skip = 0
        all_pos.each_with_index do |pos, index|
          if pos < sttpos
            skip = index
            next
          elsif sttpos <= pos && pos <= endpos
            all_pos_depth[pos] += 1  
          elsif endpos < pos
            all_pos = all_pos[skip..-1] # all_pos update
            break 
          end
        end
      end
      return all_pos_depth
    end
# count refeerence allele depth
# count total depth
# ------------------------------------------
count refeerence allele depth
count total depth
count refeerence allele depth count total depth
def count_depth(all_pos, bam_path, target_chr, mapq)
      pos_max = all_pos.max 
      pos_min = all_pos.min
      all_pos_depth = Hash.new(0)
      SAMReader.new(@samtools, bam_path).each(chr: target_chr,
                                         filter_flag: SAMReader::DUPLICATE,
                                         output_flag: SAMReader::PROPER_PAIR) do |line|
        read_seq = line.read_seq
        next if target_chr != line.chr 
        next if line.map_score <= mapq 
        next if read_seq.count("N") > 0
        map_status = line.map_status
        if /^(\d+)M$/ =~ map_status 
          match_size = $1.to_i  
        elsif /^(\d+)M(?:\d+)S$/ =~ map_status 
          len_M = $1.to_i
          match_size = len_M
        elsif /^(\d+)M(?:\d+)I(\d+)M(?:\d+)S$/ =~ map_status 
          len_M1 = $1.to_i 
          len_M2 = $2.to_i 
          match_size = len_M1 + len_M2
        elsif /^(\d+)M(\d+)D(\d+)M(?:\d+)S$/ =~ map_status 
          len_M1 = $1.to_i 
          len_D  = $2.to_i 
          len_M2 = $3.to_i 
          match_size = len_M1 + len_D + len_M2
        elsif /^(?:\d+)S(\d+)M$/ =~ map_status 
          len_M = $1.to_i
          match_size = len_M
        elsif /^(?:\d+)S(\d+)M(?:\d+)I(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_M2 = $2.to_i 
          match_size = len_M1 + len_M2
        elsif /^(?:\d+)S(\d+)M(\d+)D(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_D  = $2.to_i 
          len_M2 = $3.to_i 
          match_size = len_M1 + len_D + len_M2
        elsif /^(\d+)M(\d+)I(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          #len_I  = $2.to_i 
          len_M2 = $3.to_i 8
          match_size = len_M1 + len_M2
        elsif /^(\d+)M(\d+)D(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_D  = $2.to_i 
          len_M2 = $3.to_i 
          match_size = len_M1 + len_D + len_M2
        elsif /^(\d+)M(\d+)D(\d+)M(\d+)D(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_D1 = $2.to_i 
          len_M2 = $3.to_i 
          len_D2 = $4.to_i 
          len_M3 = $5.to_i 8
          match_size  = len_M1 + len_D1 + len_M2 + len_D2 + len_M3
        elsif /^(\d+)M(\d+)D(\d+)M(?:\d+)I(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_D1 = $2.to_i 
          len_M2 = $3.to_i 
          len_M3 = $4.to_i 8
          match_size  = len_M1 + len_D1 + len_M2 + len_M3
        elsif /^(\d+)M(?:\d+)I(\d+)M(\d+)D(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_M2 = $2.to_i 
          len_D1 = $3.to_i 
          len_M3 = $4.to_i 8
          match_size  = len_M1 + len_M2 + len_D1 + len_M3
        elsif /^(\d+)M(?:\d+)I(\d+)M(?:\d+)I(\d+)M$/ =~ map_status 
          len_M1 = $1.to_i 
          len_M2 = $2.to_i 
          len_M3 = $3.to_i 8
          match_size  = len_M1 + len_M2 + len_M3
        else
          next
        end
        sttpos = line.chrpos
        endpos = sttpos + match_size - 1
        next if endpos < pos_min
        break if pos_max < sttpos
        skip = 0
        all_pos.each_with_index do |pos, index|
          if pos < sttpos
            skip = index
            next
          elsif sttpos <= pos && pos <= endpos
            all_pos_depth[pos] += 1  
          elsif endpos < pos
            all_pos = all_pos[skip..-1] 
            break 
          end
        end
      end
      return all_pos_depth
    end
def count_depth ( all_pos , bam_path , target_chr , mapq ) pos_max = all_pos . max pos_min = all_pos . min all_pos_depth = Hash . new ( 0 ) SAMReader . new ( @samtools , bam_path ) . each ( chr : target_chr , filter_flag : SAMReader :: DUPLICATE , output_flag : SAMReader :: PROPER_PAIR ) do | line | read_seq = line . read_seq next if target_chr != line . chr next if line . map_score <= mapq next if read_seq . count ( "N" ) > 0 map_status = line . map_status if / ^( \d +)M$ / =~ map_status match_size = $1 . to_i elsif / ^( \d +)M(?: \d +)S$ / =~ map_status len_M = $1 . to_i match_size = len_M elsif / ^( \d +)M(?: \d +)I( \d +)M(?: \d +)S$ / =~ map_status len_M1 = $1 . to_i len_M2 = $2 . to_i match_size = len_M1 + len_M2 elsif / ^( \d +)M( \d +)D( \d +)M(?: \d +)S$ / =~ map_status len_M1 = $1 . to_i len_D = $2 . to_i len_M2 = $3 . to_i match_size = len_M1 + len_D + len_M2 elsif / ^(?: \d +)S( \d +)M$ / =~ map_status len_M = $1 . to_i match_size = len_M elsif / ^(?: \d +)S( \d +)M(?: \d +)I( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_M2 = $2 . to_i match_size = len_M1 + len_M2 elsif / ^(?: \d +)S( \d +)M( \d +)D( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_D = $2 . to_i len_M2 = $3 . to_i match_size = len_M1 + len_D + len_M2 elsif / ^( \d +)M( \d +)I( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_M2 = $3 . to_i match_size = len_M1 + len_M2 elsif / ^( \d +)M( \d +)D( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_D = $2 . to_i len_M2 = $3 . to_i match_size = len_M1 + len_D + len_M2 elsif / ^( \d +)M( \d +)D( \d +)M( \d +)D( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_D1 = $2 . to_i len_M2 = $3 . to_i len_D2 = $4 . to_i len_M3 = $5 . to_i match_size = len_M1 + len_D1 + len_M2 + len_D2 + len_M3 elsif / ^( \d +)M( \d +)D( \d +)M(?: \d +)I( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_D1 = $2 . to_i len_M2 = $3 . to_i len_M3 = $4 . to_i match_size = len_M1 + len_D1 + len_M2 + len_M3 elsif / ^( \d +)M(?: \d +)I( \d +)M( \d +)D( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_M2 = $2 . to_i len_D1 = $3 . to_i len_M3 = $4 . to_i match_size = len_M1 + len_M2 + len_D1 + len_M3 elsif / ^( \d +)M(?: \d +)I( \d +)M(?: \d +)I( \d +)M$ / =~ map_status len_M1 = $1 . to_i len_M2 = $2 . to_i len_M3 = $3 . to_i match_size = len_M1 + len_M2 + len_M3 else next end sttpos = line . chrpos endpos = sttpos + match_size - 1 next if endpos < pos_min break if pos_max < sttpos skip = 0 all_pos . each_with_index do | pos , index | if pos < sttpos skip = index next elsif sttpos <= pos && pos <= endpos all_pos_depth [ pos ] += 1 elsif endpos < pos all_pos = all_pos [ skip .. - 1 ] break end end end return all_pos_depth end
count refeerence allele depth
count total depth
count refeerence allele depth count total depth
# Condition1 # complete match # 161M #  obj_SAM.map_status --> 70M or 101M or # B fragment (match + clip seq) # 80M20S # B fragment (short INS + clip seq) # 70M2I8M20S #70 #8 # B fragment (short DEL + clip seq) # 70M2D10M20S #70 #2 #8 # F fragment (clip seq + match) # 20S80M # F fragment (clip seq + short INS) # 20S70M2I8M #70 #8 # F fragment (clip seq + short DEL) # 20S70M2D10M #70 #2 #8 # ----------------------------------------- # short INDEL # ----------------------------------------- # 1 insertion  # 70M2I28M #70 #len_I  = $2.to_i #2 #28 # 1 deletion # 70M5D30M #70 #5 #30 # 1 del + 1 del # 70M 2D 10M 5D 28M #70 #2 #10 #5 #28 # 1 del + 1 ins # 70M 2D 10M 5I 28M #70 #2 #10 #28 # 1 ins + 1 del # 70M 2I 10M 5D 28M #70 #10 #5 #28 # 1 ins + 1 ins # 70M 2I 10M 5I 28M #70 #10 #28 # all_pos update
all_pos
bam_path
target_chr
mapq