lib/ims_indel/consensus.rb
def detect_long_insertion(insert_pair, clip_Bs2, clip_Fs2, unmap_seq, alt_read_depth=5)
      paired_indel_list = []
      insert_pair.each do |consensus_b, consensus_f|
        total_depth = consensus_b.depth + consensus_f.depth
        clip_chrpos_all = [consensus_b.start_pos, consensus_b.end_pos, consensus_f.start_pos, consensus_f.end_pos]
        clip_start = clip_chrpos_all.min
        clip_end = clip_chrpos_all.max
        consensus, _ = mafft_consensus([consensus_b, consensus_f], 1.0) # %identity = 1.0
        consensus, trim_flag = trim_consensus(consensus)
        # ---------------------------------------------------------------------------

        if consensus.count("?") > 2 # SNP check, if the number of SNPs is more than 2 ( > 2), uncomplete Long Insertion
          #puts "uncomplete long insertion..."

          # re-alignement using unmapped read
          can_use_unmaps = unmap_seq.find_all do | unmap|
            (unmap.start_pos < clip_start && clip_start < unmap.end_pos) || (unmap.start_pos < clip_end && clip_end < unmap.end_pos)
          end
          if can_use_unmaps.empty? # unmap read is nothing
            clip_B_stt = consensus_b.start_pos # clip_B start pos
            clip_F_end = consensus_f.end_pos # clip_F end pos
            upper_LI  = clip_Bs2[clip_B_stt]
            bottom_LI = clip_Fs2[clip_F_end]
            consensus = "#{upper_LI}-----#{bottom_LI}"
            if total_depth >= alt_read_depth # end pos is equal to sttpos
              paired_indel_list << Read.new(type: :ULI, start_pos: clip_B_stt, end_pos: clip_B_stt, seq: consensus, depth: total_depth)
              @mafft_inputs[paired_indel_list.last] = [consensus_b, consensus_f] if @mafft_inputs
            end

          else # unmap reads exist
            new_group_reads = ([consensus_b, consensus_f] + can_use_unmaps)
            new_consensus, _ = mafft_consensus(new_group_reads, 1.0) # make a consensus seq with the unmapped reads
            new_consensus, trim_flag = trim_consensus(new_consensus)

            if new_consensus.empty? || new_consensus.count("?") > 2 # SNP check, if the number of SNPs is more than 2 ( >=3)
              clip_B_stt = consensus_b.start_pos # clip_B start pos
              clip_F_end = consensus_f.end_pos # clip_F end pos
              upper_LI  = clip_Bs2[clip_B_stt]
              bottom_LI = clip_Fs2[clip_F_end]
              consensus = "#{upper_LI}-----#{bottom_LI}"
              if total_depth >= alt_read_depth # end pos is equal to sttpos
                paired_indel_list << Read.new(type: :ULI, start_pos: clip_B_stt, end_pos: clip_B_stt, seq: consensus, depth: total_depth)
                @mafft_inputs[paired_indel_list.last] = new_group_reads if @mafft_inputs
              end
            else
              new_consensus, _ = mafft_consensus(new_group_reads, 0.5) # make a consensus seq with the unmapped reads
              new_consensus = trim_consensus_with_flag(new_consensus, trim_flag)

              new_total_depth = new_group_reads.inject(0){|res, read| res += read.depth}
              if new_total_depth >= alt_read_depth ####
                paired_indel_list << Read.new(type: :LI_wU, start_pos: clip_start, end_pos: clip_end, seq: new_consensus, depth: new_total_depth)
                @mafft_inputs[paired_indel_list.last] = new_group_reads if @mafft_inputs
              end
            end
          end

        else
          consensus, _ = mafft_consensus([consensus_b, consensus_f], 0.5) # %identity = 0.5, consensus update
          consensus = trim_consensus_with_flag(consensus, trim_flag)

          #puts "complete long insertion..."
          #puts align_reads_names
          #puts
          if total_depth >= alt_read_depth ####
            paired_indel_list << Read.new(type: :LI, start_pos: clip_start, end_pos: clip_end, seq: consensus, depth: total_depth)
            @mafft_inputs[paired_indel_list.last] = [consensus_b, consensus_f] if @mafft_inputs
          end
        end
      end
      return paired_indel_list
    end
# from insert_pairs, find complete insertion and uncomplete insertion
from insert_pairs, find complete insertion and uncomplete insertion
from insert_pairs find complete insertion and uncomplete insertion
def detect_long_insertion(insert_pair, clip_Bs2, clip_Fs2, unmap_seq, alt_read_depth=5)
      paired_indel_list = []
      insert_pair.each do |consensus_b, consensus_f|
        total_depth = consensus_b.depth + consensus_f.depth
        clip_chrpos_all = [consensus_b.start_pos, consensus_b.end_pos, consensus_f.start_pos, consensus_f.end_pos]
        clip_start = clip_chrpos_all.min
        clip_end = clip_chrpos_all.max
        consensus, _ = mafft_consensus([consensus_b, consensus_f], 1.0) 
        consensus, trim_flag = trim_consensus(consensus)
        if consensus.count("?") > 2 
          can_use_unmaps = unmap_seq.find_all do | unmap|
            (unmap.start_pos < clip_start && clip_start < unmap.end_pos) || (unmap.start_pos < clip_end && clip_end < unmap.end_pos)
          end
          if can_use_unmaps.empty? 
            clip_B_stt = consensus_b.start_pos 
            clip_F_end = consensus_f.end_pos 
            upper_LI  = clip_Bs2[clip_B_stt]
            bottom_LI = clip_Fs2[clip_F_end]
            consensus = "#{upper_LI}-----#{bottom_LI}"
            if total_depth >= alt_read_depth 
              paired_indel_list << Read.new(type: :ULI, start_pos: clip_B_stt, end_pos: clip_B_stt, seq: consensus, depth: total_depth)
              @mafft_inputs[paired_indel_list.last] = [consensus_b, consensus_f] if @mafft_inputs
            end
          else 
            new_group_reads = ([consensus_b, consensus_f] + can_use_unmaps)
            new_consensus, _ = mafft_consensus(new_group_reads, 1.0) 
            new_consensus, trim_flag = trim_consensus(new_consensus)
            if new_consensus.empty? || new_consensus.count("?") > 2 
              clip_B_stt = consensus_b.start_pos 
              clip_F_end = consensus_f.end_pos 
              upper_LI  = clip_Bs2[clip_B_stt]
              bottom_LI = clip_Fs2[clip_F_end]
              consensus = "#{upper_LI}-----#{bottom_LI}"
              if total_depth >= alt_read_depth 
                paired_indel_list << Read.new(type: :ULI, start_pos: clip_B_stt, end_pos: clip_B_stt, seq: consensus, depth: total_depth)
                @mafft_inputs[paired_indel_list.last] = new_group_reads if @mafft_inputs
              end
            else
              new_consensus, _ = mafft_consensus(new_group_reads, 0.5) 
              new_consensus = trim_consensus_with_flag(new_consensus, trim_flag)
              new_total_depth = new_group_reads.inject(0){|res, read| res += read.depth}
              if new_total_depth >= alt_read_depth 
                paired_indel_list << Read.new(type: :LI_wU, start_pos: clip_start, end_pos: clip_end, seq: new_consensus, depth: new_total_depth)
                @mafft_inputs[paired_indel_list.last] = new_group_reads if @mafft_inputs
              end
            end
          end
        else
          consensus, _ = mafft_consensus([consensus_b, consensus_f], 0.5) 
          consensus = trim_consensus_with_flag(consensus, trim_flag)
          if total_depth >= alt_read_depth 
            paired_indel_list << Read.new(type: :LI, start_pos: clip_start, end_pos: clip_end, seq: consensus, depth: total_depth)
            @mafft_inputs[paired_indel_list.last] = [consensus_b, consensus_f] if @mafft_inputs
          end
        end
      end
      return paired_indel_list
    end
def detect_long_insertion ( insert_pair , clip_Bs2 , clip_Fs2 , unmap_seq , alt_read_depth = 5 ) paired_indel_list = [ ] insert_pair . each do | consensus_b , consensus_f | total_depth = consensus_b . depth + consensus_f . depth clip_chrpos_all = [ consensus_b . start_pos , consensus_b . end_pos , consensus_f . start_pos , consensus_f . end_pos ] clip_start = clip_chrpos_all . min clip_end = clip_chrpos_all . max consensus , _ = mafft_consensus ( [ consensus_b , consensus_f ] , 1.0 ) consensus , trim_flag = trim_consensus ( consensus ) if consensus . count ( "?" ) > 2 can_use_unmaps = unmap_seq . find_all do | unmap | ( unmap . start_pos < clip_start && clip_start < unmap . end_pos ) || ( unmap . start_pos < clip_end && clip_end < unmap . end_pos ) end if can_use_unmaps . empty? clip_B_stt = consensus_b . start_pos clip_F_end = consensus_f . end_pos upper_LI = clip_Bs2 [ clip_B_stt ] bottom_LI = clip_Fs2 [ clip_F_end ] consensus = "#{upper_LI}-----#{bottom_LI}" if total_depth >= alt_read_depth paired_indel_list << Read . new ( type : :ULI , start_pos : clip_B_stt , end_pos : clip_B_stt , seq : consensus , depth : total_depth ) @mafft_inputs [ paired_indel_list . last ] = [ consensus_b , consensus_f ] if @mafft_inputs end else new_group_reads = ( [ consensus_b , consensus_f ] + can_use_unmaps ) new_consensus , _ = mafft_consensus ( new_group_reads , 1.0 ) new_consensus , trim_flag = trim_consensus ( new_consensus ) if new_consensus . empty? || new_consensus . count ( "?" ) > 2 clip_B_stt = consensus_b . start_pos clip_F_end = consensus_f . end_pos upper_LI = clip_Bs2 [ clip_B_stt ] bottom_LI = clip_Fs2 [ clip_F_end ] consensus = "#{upper_LI}-----#{bottom_LI}" if total_depth >= alt_read_depth paired_indel_list << Read . new ( type : :ULI , start_pos : clip_B_stt , end_pos : clip_B_stt , seq : consensus , depth : total_depth ) @mafft_inputs [ paired_indel_list . last ] = new_group_reads if @mafft_inputs end else new_consensus , _ = mafft_consensus ( new_group_reads , 0.5 ) new_consensus = trim_consensus_with_flag ( new_consensus , trim_flag ) new_total_depth = new_group_reads . inject ( 0 ) { | res , read | res += read . depth } if new_total_depth >= alt_read_depth paired_indel_list << Read . new ( type : :LI_wU , start_pos : clip_start , end_pos : clip_end , seq : new_consensus , depth : new_total_depth ) @mafft_inputs [ paired_indel_list . last ] = new_group_reads if @mafft_inputs end end end else consensus , _ = mafft_consensus ( [ consensus_b , consensus_f ] , 0.5 ) consensus = trim_consensus_with_flag ( consensus , trim_flag ) if total_depth >= alt_read_depth paired_indel_list << Read . new ( type : :LI , start_pos : clip_start , end_pos : clip_end , seq : consensus , depth : total_depth ) @mafft_inputs [ paired_indel_list . last ] = [ consensus_b , consensus_f ] if @mafft_inputs end end end return paired_indel_list end
from insert_pairs, find complete insertion and uncomplete insertion
from insert_pairs find complete insertion and uncomplete insertion
# %identity = 1.0 # --------------------------------------------------------------------------- # SNP check, if the number of SNPs is more than 2 ( > 2), uncomplete Long Insertion #puts "uncomplete long insertion..." # re-alignement using unmapped read # unmap read is nothing # clip_B start pos # clip_F end pos # end pos is equal to sttpos # unmap reads exist # make a consensus seq with the unmapped reads # SNP check, if the number of SNPs is more than 2 ( >=3) # clip_B start pos # clip_F end pos # end pos is equal to sttpos # make a consensus seq with the unmapped reads #### # %identity = 0.5, consensus update #puts "complete long insertion..." #puts align_reads_names #puts ####
insert_pair
clip_Bs2
clip_Fs2
unmap_seq
alt_read_depth