lib/ims_indel/indel_detector.rb
def detect_indels(indel_list, target_chr, max_indel_size, reffa)
      fin_indels = Hash.new {|h,k| h[k] = [] } # {31973646_31973663_tgttatttataataatga" => ["INS_SID_18(depth)_indelsttpos_seq", ..
      @is_checked_pos = {} # same Indel check (temprary use)
      all_pos    = []
      pairwise_consensus = PairwiseConsensus.new(@glsearch_bin, @glsearch_mat, @temp_path)

      #ttl_indel = indel_list.size
      indel_list.each_with_index do |read, index| # read = [":LI(or :ULI or :LD or :B or :F or :SID), clip_chrpos_all.min, clip_chrpos_all.max, consensus, ttl_depth]
        flag = read.type
        sttpos = read.start_pos
        endpos = read.end_pos
        seq = read.seq
        depth = read.depth

        if flag == :ULI # uncomplete long indel
          key = "#{sttpos}_#{endpos}"
          if @is_checked_pos[key]
            #puts "#{index+1}/#{ttl_indel}\t#{flag}_#{sttpos}-#{endpos}\t#{seq}\t#{seq.size-5}+\tdepth: #{depth}: INS\tSAME\n"
            fin_indels[key] << IndelResult.new("INS", flag, depth, sttpos, seq, read)
            all_pos << sttpos
          else
            @is_checked_pos[key] = true
            #puts "#{index+1}/#{ttl_indel}\t#{flag}_#{sttpos}-#{endpos}\t#{seq}\t#{seq.size-5}+\tdepth: #{depth}: INS\n"
            fin_indels[key] << IndelResult.new("INS", flag, depth, sttpos, seq, read)
            all_pos << sttpos
          end
          next
        end

        # cut a part of reference seq
        if sttpos - (max_indel_size/2) < 1
          ref_stt = 1
        else
          ref_stt = sttpos - (max_indel_size/2)
        end
        ref_end = endpos + (max_indel_size/2)
        ref_seq = `samtools faidx #{reffa} #{target_chr}:#{ref_stt}-#{ref_end}`.chomp.split("\n")[1..-1].join("").upcase
        ref  = [:R, ref_stt, ref_end, ref_seq]
        ipt_read = [read.type, read.start_pos, read.end_pos, read.seq]
        group_reads = [ipt_read, ref]
        indel_stt = ipt_read[1] # indel start position before comparison of reference genome, #######################
        consensus, align_reads_names, trim_flag = pairwise_consensus.make_consensus(group_reads, 1.0, nil)
        next if consensus.nil?

        if consensus.empty? || consensus.count("?") > 2 # ?2個はO.K.
          #align_reads_names.each_with_index do |line, subindex|
            #if subindex == 2
              #puts "#{index+1}/#{ttl_indel}\t#{line}\tNOT!!"
            #elsif /^R_/ =~ line.chomp.split("\t")[-1]
              #puts "Ref\t#{line}\n"
            #else
              #puts "Seq\t#{line}\n"
            #end
          #end
          next
        end

        consensus, align_reads_names = pairwise_consensus.make_consensus(group_reads, 0.5, trim_flag)
        next if consensus.nil?

        read_seq, _ = align_reads_names[0].split("\t")
        ref_seq, ref_name   = align_reads_names[1].split("\t")

        # special for LI, 150525
        # ref_seq = cagaaaaa--aaaaa---------ccagtcccgttctaaccagtcccgttctaa, :first --はmiss 
        #puts "ref\n#{ref_seq}\nread\n#{read_seq}\n"

        if (flag == :LI || flag == :LI_wU) && /(^\w+)(-+)(\w+)(-+)(\w+)$/ =~ ref_seq
          #puts "LI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n#{flag}\n"
          fir_str = $1.size
          fir_gap = $2.size
          sec_str = $3.size
          sec_gap = $4.size
          next unless /(^\-+)(\w+)(\-+)$/ =~ read_seq
          indel_info = process_long_insersion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
          fin_indels[indel_info.key] << IndelResult.new("INS", flag, depth, indel_stt, indel_info.str, read)
          all_pos << indel_stt
          #puts_indel_info("INS", ref_seq, read_seq, index, ttl_indel, align_reads_names, flag, depth, indel_info)

        elsif /(^\w+)(-+)(\w+)$/ =~ ref_seq # ref = cagaaaaa-----ccagtcccgttctaaccagtcccgttctaa
          # insertion: 右寄せのパターン or 真ん中が
          # 1. mono polymer
          # read = "TTTaaaaGGGGGaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
          # ref  = "---aaaaGGG--aaaaaaaaaaaaaaaaaaaaaa-------"
          # 2. multi polymer (右寄せ)
          # ref  = "TTTaaaa<a><b> -  - TTTT"
          # read = "---aaaa<a><b><a><b>TTTT"
          fir_str = $1.size
          fir_gap = $2.size
          sec_str = $3.size
          next unless /(^\-+)(\w+)(\-+)$/ =~ read_seq

          if no_overlap?(ref_seq, read_seq, fir_str, fir_gap)
            #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}"
            #puts "#{index+1}/#{ttl_indel}\tNO Overlap(Ins)!!\n"
            next
          end
          if few_over_lap_insertion?(read_seq, fir_str, fir_gap)
            #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}"
            #puts "#{index+1}/#{ttl_indel}\tFew Overlap(Ins)!!"
            next
          end
          indel_info = process_right_insersion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
          fin_indels[indel_info.key] << IndelResult.new("INS", flag, depth, indel_stt, indel_info.str, read)
          all_pos << indel_stt
          #puts_indel_info("INS", ref_seq, read_seq, index, ttl_indel, align_reads_names, flag, depth, indel_info)

        elsif ref_seq.count("-") == 0
          # deletion: 右寄せのパターンが多い
          # 1. mono polymer
          # ref  = "TTTaaaaGGGGGaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
          # read = "---aaaaGGG--aaaaaaaaaaaaaaaaaaaaaa-------"
          # 2. multi polymer
          # ref  = "TTTaaaaATGAGTAGTaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
          # read = "---aaaaAGTAGT---aaaaaaaaaaaaaaaaaaaaaa-------"
          #
          # first str >= 10 and sec str >= 10
          if /(^\-+)(\w+)(\-+)(\w+)(?:\-+)$/ =~ read_seq
            fir_gap  = $1.size
            fir_str  = $2.size
            sec_gap  = $3.size
            sec_str  = $4.size
            if few_over_lap_deletion?(fir_str, sec_str)
              #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}"
              #puts "#{index+1}/#{ttl_indel}\tFew Overlap(Del)!!\n"
              next
            end
            indel_info = process_deletion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
            fin_indels[indel_info.key] << IndelResult.new("DEL", flag, depth, indel_stt, indel_info.str, read)
            all_pos << indel_stt
          else
            #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}"
            #puts "#{index+1}/#{ttl_indel}\t#{align_reads_names[-1]}"
            #puts "#{index+1}/#{ttl_indel}\t#{flag}\tNOT_INDEL\n"
          end
        end
      end
      @is_checked_pos = nil
      all_pos.sort!
      all_pos.uniq!
      return all_pos, fin_indels
    end
# compared with reference seq, detect indels precisely
compared with reference seq, detect indels precisely
compared with reference seq detect indels precisely
def detect_indels(indel_list, target_chr, max_indel_size, reffa)
      fin_indels = Hash.new {|h,k| h[k] = [] } 
      @is_checked_pos = {} 
      all_pos    = []
      pairwise_consensus = PairwiseConsensus.new(@glsearch_bin, @glsearch_mat, @temp_path)
      indel_list.each_with_index do |read, index| 
        flag = read.type
        sttpos = read.start_pos
        endpos = read.end_pos
        seq = read.seq
        depth = read.depth
        if flag == :ULI 
          key = "{sttpos}_{endpos}"
          if @is_checked_pos[key]
            fin_indels[key] << IndelResult.new("INS", flag, depth, sttpos, seq, read)
            all_pos << sttpos
          else
            @is_checked_pos[key] = true
            fin_indels[key] << IndelResult.new("INS", flag, depth, sttpos, seq, read)
            all_pos << sttpos
          end
          next
        end
        if sttpos - (max_indel_size/2) < 1
          ref_stt = 1
        else
          ref_stt = sttpos - (max_indel_size/2)
        end
        ref_end = endpos + (max_indel_size/2)
        ref_seq = `samtools faidx {reffa} {target_chr}:{ref_stt}-{ref_end}`.chomp.split("\n")[1..-1].join("").upcase
        ref  = [:R, ref_stt, ref_end, ref_seq]
        ipt_read = [read.type, read.start_pos, read.end_pos, read.seq]
        group_reads = [ipt_read, ref]
        indel_stt = ipt_read[1] 
        consensus, align_reads_names, trim_flag = pairwise_consensus.make_consensus(group_reads, 1.0, nil)
        next if consensus.nil?
        if consensus.empty? || consensus.count("?") > 2 
          next
        end
        consensus, align_reads_names = pairwise_consensus.make_consensus(group_reads, 0.5, trim_flag)
        next if consensus.nil?
        read_seq, _ = align_reads_names[0].split("\t")
        ref_seq, ref_name   = align_reads_names[1].split("\t")
        if (flag == :LI || flag == :LI_wU) && /(^\w+)(-+)(\w+)(-+)(\w+)$/ =~ ref_seq
          fir_str = $1.size
          fir_gap = $2.size
          sec_str = $3.size
          sec_gap = $4.size
          next unless /(^\-+)(\w+)(\-+)$/ =~ read_seq
          indel_info = process_long_insersion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
          fin_indels[indel_info.key] << IndelResult.new("INS", flag, depth, indel_stt, indel_info.str, read)
          all_pos << indel_stt
        elsif /(^\w+)(-+)(\w+)$/ =~ ref_seq 
          fir_str = $1.size
          fir_gap = $2.size
          sec_str = $3.size
          next unless /(^\-+)(\w+)(\-+)$/ =~ read_seq
          if no_overlap?(ref_seq, read_seq, fir_str, fir_gap)
            next
          end
          if few_over_lap_insertion?(read_seq, fir_str, fir_gap)
            next
          end
          indel_info = process_right_insersion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
          fin_indels[indel_info.key] << IndelResult.new("INS", flag, depth, indel_stt, indel_info.str, read)
          all_pos << indel_stt
        elsif ref_seq.count("-") == 0
           first str >= 10 and sec str >= 10
          if /(^\-+)(\w+)(\-+)(\w+)(?:\-+)$/ =~ read_seq
            fir_gap  = $1.size
            fir_str  = $2.size
            sec_gap  = $3.size
            sec_str  = $4.size
            if few_over_lap_deletion?(fir_str, sec_str)
              puts "{index+1}/{ttl_indel}\tFew Overlap(Del)!!\n"
              next
            end
            indel_info = process_deletion(fir_gap, fir_str, sec_gap, sec_str, read_seq, ref_name, ref_seq)
            fin_indels[indel_info.key] << IndelResult.new("DEL", flag, depth, indel_stt, indel_info.str, read)
            all_pos << indel_stt
          else
            puts "{index+1}/{ttl_indel}\t{align_reads_names[-1]}"
            puts "{index+1}/{ttl_indel}\t{flag}\tNOT_INDEL\n"
          end
        end
      end
      @is_checked_pos = nil
      all_pos.sort!
      all_pos.uniq!
      return all_pos, fin_indels
    end
def detect_indels ( indel_list , target_chr , max_indel_size , reffa ) fin_indels = Hash . new { | h , k | h [ k ] = [ ] } @is_checked_pos = { } all_pos = [ ] pairwise_consensus = PairwiseConsensus . new ( @glsearch_bin , @glsearch_mat , @temp_path ) indel_list . each_with_index do | read , index | flag = read . type sttpos = read . start_pos endpos = read . end_pos seq = read . seq depth = read . depth if flag == :ULI key = "#{sttpos}_#{endpos}" if @is_checked_pos [ key ] fin_indels [ key ] << IndelResult . new ( "INS" , flag , depth , sttpos , seq , read ) all_pos << sttpos else @is_checked_pos [ key ] = true fin_indels [ key ] << IndelResult . new ( "INS" , flag , depth , sttpos , seq , read ) all_pos << sttpos end next end if sttpos - ( max_indel_size / 2 ) < 1 ref_stt = 1 else ref_stt = sttpos - ( max_indel_size / 2 ) end ref_end = endpos + ( max_indel_size / 2 ) ref_seq = ` samtools faidx  #{ reffa } #{ target_chr } : #{ ref_stt } - #{ ref_end } ` . chomp . split ( "\n" ) [ 1 .. - 1 ] . join ( "" ) . upcase ref = [ :R , ref_stt , ref_end , ref_seq ] ipt_read = [ read . type , read . start_pos , read . end_pos , read . seq ] group_reads = [ ipt_read , ref ] indel_stt = ipt_read [ 1 ] consensus , align_reads_names , trim_flag = pairwise_consensus . make_consensus ( group_reads , 1.0 , nil ) next if consensus . nil? if consensus . empty? || consensus . count ( "?" ) > 2 next end consensus , align_reads_names = pairwise_consensus . make_consensus ( group_reads , 0.5 , trim_flag ) next if consensus . nil? read_seq , _ = align_reads_names [ 0 ] . split ( "\t" ) ref_seq , ref_name = align_reads_names [ 1 ] . split ( "\t" ) if ( flag == :LI || flag == :LI_wU ) && / (^ \w +)(-+)( \w +)(-+)( \w +)$ / =~ ref_seq fir_str = $1 . size fir_gap = $2 . size sec_str = $3 . size sec_gap = $4 . size next unless / (^ \- +)( \w +)( \- +)$ / =~ read_seq indel_info = process_long_insersion ( fir_gap , fir_str , sec_gap , sec_str , read_seq , ref_name , ref_seq ) fin_indels [ indel_info . key ] << IndelResult . new ( "INS" , flag , depth , indel_stt , indel_info . str , read ) all_pos << indel_stt elsif / (^ \w +)(-+)( \w +)$ / =~ ref_seq fir_str = $1 . size fir_gap = $2 . size sec_str = $3 . size next unless / (^ \- +)( \w +)( \- +)$ / =~ read_seq if no_overlap? ( ref_seq , read_seq , fir_str , fir_gap ) next end if few_over_lap_insertion? ( read_seq , fir_str , fir_gap ) next end indel_info = process_right_insersion ( fir_gap , fir_str , sec_gap , sec_str , read_seq , ref_name , ref_seq ) fin_indels [ indel_info . key ] << IndelResult . new ( "INS" , flag , depth , indel_stt , indel_info . str , read ) all_pos << indel_stt elsif ref_seq . count ( "-" ) == 0 if / (^ \- +)( \w +)( \- +)( \w +)(?: \- +)$ / =~ read_seq fir_gap = $1 . size fir_str = $2 . size sec_gap = $3 . size sec_str = $4 . size if few_over_lap_deletion? ( fir_str , sec_str ) next end indel_info = process_deletion ( fir_gap , fir_str , sec_gap , sec_str , read_seq , ref_name , ref_seq ) fin_indels [ indel_info . key ] << IndelResult . new ( "DEL" , flag , depth , indel_stt , indel_info . str , read ) all_pos << indel_stt else end end end @is_checked_pos = nil all_pos . sort! all_pos . uniq! return all_pos , fin_indels end
compared with reference seq, detect indels precisely
compared with reference seq detect indels precisely
# {31973646_31973663_tgttatttataataatga" => ["INS_SID_18(depth)_indelsttpos_seq", .. # same Indel check (temprary use) #ttl_indel = indel_list.size # read = [":LI(or :ULI or :LD or :B or :F or :SID), clip_chrpos_all.min, clip_chrpos_all.max, consensus, ttl_depth] # uncomplete long indel #puts "#{index+1}/#{ttl_indel}\t#{flag}_#{sttpos}-#{endpos}\t#{seq}\t#{seq.size-5}+\tdepth: #{depth}: INS\tSAME\n" #puts "#{index+1}/#{ttl_indel}\t#{flag}_#{sttpos}-#{endpos}\t#{seq}\t#{seq.size-5}+\tdepth: #{depth}: INS\n" # cut a part of reference seq # indel start position before comparison of reference genome, ####################### # ?2個はO.K. #align_reads_names.each_with_index do |line, subindex| #if subindex == 2 #puts "#{index+1}/#{ttl_indel}\t#{line}\tNOT!!" #elsif /^R_/ =~ line.chomp.split("\t")[-1] #puts "Ref\t#{line}\n" #else #puts "Seq\t#{line}\n" #end #end # special for LI, 150525 # ref_seq = cagaaaaa--aaaaa---------ccagtcccgttctaaccagtcccgttctaa, :first --はmiss  #puts "ref\n#{ref_seq}\nread\n#{read_seq}\n" #puts "LI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n#{flag}\n" #puts_indel_info("INS", ref_seq, read_seq, index, ttl_indel, align_reads_names, flag, depth, indel_info) # ref = cagaaaaa-----ccagtcccgttctaaccagtcccgttctaa # insertion: 右寄せのパターン or 真ん中が # 1. mono polymer # read = "TTTaaaaGGGGGaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" # ref  = "---aaaaGGG--aaaaaaaaaaaaaaaaaaaaaa-------" # 2. multi polymer (右寄せ) # ref  = "TTTaaaa<a><b> -  - TTTT" # read = "---aaaa<a><b><a><b>TTTT" #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}" #puts "#{index+1}/#{ttl_indel}\tNO Overlap(Ins)!!\n" #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}" #puts "#{index+1}/#{ttl_indel}\tFew Overlap(Ins)!!" #puts_indel_info("INS", ref_seq, read_seq, index, ttl_indel, align_reads_names, flag, depth, indel_info) # deletion: 右寄せのパターンが多い # 1. mono polymer # ref  = "TTTaaaaGGGGGaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" # read = "---aaaaGGG--aaaaaaaaaaaaaaaaaaaaaa-------" # 2. multi polymer # ref  = "TTTaaaaATGAGTAGTaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" # read = "---aaaaAGTAGT---aaaaaaaaaaaaaaaaaaaaaa-------" # # first str >= 10 and sec str >= 10 #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}" #puts "#{index+1}/#{ttl_indel}\tFew Overlap(Del)!!\n" #puts "Ref\t#{ref_seq}\nSeq\t#{read_seq}" #puts "#{index+1}/#{ttl_indel}\t#{align_reads_names[-1]}" #puts "#{index+1}/#{ttl_indel}\t#{flag}\tNOT_INDEL\n"
indel_list
target_chr
max_indel_size
reffa