shyamrallapalli/cheripic

# encoding: utf-8
require 'bio'
require 'forwardable'

module Cheripic

# Custom error handling for Contig class
  class ContigError < CheripicError; end

  # A contig object from assembly that stores positions of
  # homozygous, heterozygous and hemi-variants
  #
  # @!attribute [rw] hm_pos
  #   @return [Hash] a hash of homozygous variant positions as keys and allele frequency as values
  # @!attribute [rw] ht_pos
  #   @return [Hash] a hash of heterozygous variant positions as keys and allele frequency as values
  # @!attribute [rw] hemi_pos
  #   @return [Hash] a hash of hemi-variant positions as keys and allele frequency as values
  # @!attribute [r] id
  #   @return [String] id of the contig in assembly taken from fasta file
  # @!attribute [r] length
  #   @return [Integer] length of contig in bases
  class Contig

    attr_accessor :hm_pos, :ht_pos, :hemi_pos, :mean_depth, :sd_depth
    attr_reader :id, :length

    # creates a Contig object using fasta entry
    # @param fasta [Bio::FastaFormat] an individual fasta entry from input assembly file
    def initialize (fasta)
      @id = fasta.entry_id
      @length = fasta.length
      @hm_pos = {}
      @ht_pos = {}
      @hemi_pos = {}
      @mean_depth = nil
      @sd_depth = nil
    end

    # Number of homozygous variants identified in the contig
    # @return [Integer]
    def hm_num
      self.hm_pos.length
    end

    # Number of heterozygous variants identified in the contig
    # @return [Integer]
    def ht_num
      self.ht_pos.length
    end

    # Homozygosity enrichment score calculated using
    # hm_num and ht_num of the contig object
    # @return [Float]
    def hme_score
      hmes_adjust = Options.hmes_adjust
      if self.hm_num == 0 and self.ht_num == 0
        0.0
      else
        (self.hm_num + hmes_adjust) / (self.ht_num + hmes_adjust)
      end
    end

    # Number of hemi-variants identified in the contig
    # @return [Integer]
    def hemi_num
      self.hemi_pos.length
    end

    # Mean of bulk frequency ratios (bfr) calculated using
    # bfr values all hemi_pos of the contig
    # @return [Float]
    def bfr_score
      if self.hemi_pos.values.empty?
        0.0
      else
        geom_mean(self.hemi_pos.values)
      end
    end

    # Calculates mean of an array of numbers
    # @param array [Array] an array of bfr values from hemi_snp
    # @return [Float] mean value as float
    def geom_mean(array)
      return array[0].to_f if array.length == 1
      array.reduce(:+) / array.size.to_f
      # sum = 0.0
      # array.each{ |v| sum += Math.log(v.to_f) }
      # sum /= array.size
      # Math.exp sum
    end

  end # Contig

end # Cheripic