class Bio::BFRTools::BFRContainer

Public Instance Methods

get_region(opts={}) click to toggle source
# File lib/bio/BFRTools.rb, line 386
def get_region(opts={})
  opts[:container] = self
  region = BFRRegion.new(opts)
end
init_counters() click to toggle source
# File lib/bio/BFRTools.rb, line 362
def init_counters
  @putative_snps      = 0
  @proccesed_regions  = 0
  @not_enogh_coverage = 0
  @total_avg_coverage_bulk_1 = 0.0
  @total_avg_coverage_bulk_2 = 0.0
  @total_snp_1kbp = 0.0
  @no_snps = 0
  @too_many_snps = 0

end
print_header(opts={}) click to toggle source
print_stats(opts={}) click to toggle source
process_region(opts={}) click to toggle source
# File lib/bio/BFRTools.rb, line 391
def process_region(opts={})        
  opts = { :min_cov=>20, :max_snp_1kbp => 10, :max_per=>0.20 }.merge!(opts)

  @proccesed_regions += 1
  output = opts[:output_file] ? opts[:output_file] : $stdout
  print_output = opts[:output_file] ? true : false
  opts[:container] = self

  region = BFRRegion.new(opts)

  #puts region.to_multi_fasta

  @total_snp_1kbp += region.snp_1kbp 
 # puts "SNPS: #{region.snp_1kbp}"
  if region.snp_count == 0
    @no_snps += 1 
    print_output = false 
  end

  if region.snp_1kbp  > opts[:max_snp_1kbp]
    @too_many_snps += 1
    print_output = false
  end



  @total_avg_coverage_bulk_2 += region.avg_cov_bulk_2
  @total_avg_coverage_bulk_1 += region.avg_cov_bulk_1

  if region.avg_cov_bulk_2 < opts[:min_cov] or region.avg_cov_bulk_1 < opts[:min_cov]
    @not_enogh_coverage += 1
    print_output = false
  end

  info = Array.new

  if print_output
    for i in (0..region.size-1)
      if region.coverages_1[i] > opts[:min_cov] and region.coverages_2[i] > opts[:min_cov]
        BASES.each do |base|

          info.clear
          if  Bio::NucleicAcid.is_valid( region.parental_1_sequence[i],  base.to_s  ) and 
            not  Bio::NucleicAcid.is_valid( region.parental_2_sequence[i],  base.to_s  )
            info << :first
          end

          if   Bio::NucleicAcid.is_valid( region.parental_2_sequence[i],  base.to_s  ) and 
            not Bio::NucleicAcid.is_valid( region.parental_1_sequence[i],  base.to_s  )
            info << :second
          end


          for informative in info
            line = region.get_bfr_line(i, base, informative)
            output.print  line , "\n"
          end
        end
      end
    end
  end


  @parental_1_sam.mpileup_clear_cache region
  @parental_2_sam.mpileup_clear_cache region
  @bulk_2_sam.mpileup_clear_cache region
  @bulk_1_sam.mpileup_clear_cache region
  return region
end