R/hwe_geno.R

In agdamsbo/daDoctoR: Functions For Health Research

#' HWE calculation by genotypes
#'
#' For easy calculation.
#' @param mm First count of genotype.
#' @keywords hardy-weinberg-equllibrium
#' @export

hwe_geno<-function(mm,mn,nn,mo,no,oo,alleles=2)
{
  ## x is the number of common homozygote, y the heterozygote and z the rare homozygote.

  if (alleles==2){
    ## Biallelic has three degrees of freedom
    df=1

    a<-sum(mm,mn,nn)*2
    b<-sum(mm,mn,nn)
    p<-(2*mm+mn)/a
    q<-(2*nn+mn)/a

    al_dist<-round(rbind(cbind(p*a,q*a),cbind(p*100,q*100)),1)
    al_names<-c("m","n")

    p_p<-round(p^2*b,3)
    p_q<-round(2*p*q*b,3)
    q_q<-round(q^2*b,3)

    obs=rbind(mm,mn,nn)
    exp=rbind(p_p, p_q, q_q)

    hwe<-data.frame(obs,exp)

    hwe$dev<-hwe$obs-hwe$exp
    hwe$chi<-hwe$dev^2/hwe$exp

    snp_obs<-round(rbind(hwe$obs,hwe$obs/sum(hwe$obs)*100),1)

    snp_exp<-round(rbind(hwe$exp,hwe$exp/b*100),1)

    gen_names<-c("mm","mn","nn")

    chi<-sum(hwe$chi,na.rm = TRUE)

    p_v<-pchisq(chi, df=df, lower.tail=FALSE)
  }
  if(alleles==3){
    ## Triallelic has three degrees of freedom
    df=3

    a<-sum(mm,mn,nn,mo,no,oo)*2
    b<-sum(mm,mn,nn,mo,no,oo)
    p<-(2*mm+mn+mo)/a
    q<-(2*nn+mn+no)/a
    r<-(2*oo+no+mo)/a

    al_dist<-round(rbind(cbind(p*a,q*a,r*a),cbind(p*100,q*100,r*100)),1)
    al_names<-c("m","n","o")

    p_p<-round(p^2*b,3)
    p_q<-round(2*p*q*b,3)
    q_q<-round(q^2*b,3)
    p_r<-round(2*r*p*b,3)
    q_r<-round(2*q*r*b,3)
    r_r<-round(r^2*b,3)

    obs=rbind(mm,mn,nn,mo,no,oo)
    exp=rbind(p_p, p_q, q_q, p_r, q_r, r_r)

    hwe<-data.frame(obs,exp)

    hwe$dev<-hwe$obs-hwe$exp
    hwe$chi<-hwe$dev^2/hwe$exp

    snp_obs<-round(rbind(hwe$obs,hwe$obs/sum(hwe$obs)*100),1)

    snp_exp<-round(rbind(hwe$exp,hwe$exp/b*100),1)

    gen_names<-c("mm","mn","nn", "mo", "no", "oo")

    chi<-sum(hwe$chi,na.rm = TRUE)

    p_v<-pchisq(chi, df=df, lower.tail=FALSE)
  }


  colnames(al_dist)<-al_names
  colnames(snp_obs)<-gen_names
  colnames(snp_exp)<-gen_names

  rownames(al_dist)<-c("N","%")
  rownames(snp_obs)<-c("N","%")
  rownames(snp_exp)<-c("N","%")

  int<-ifelse(p_v<=0.05,"The null-hypothesis of difference from the HWE can be confirmed","The null-hypothesis of difference from the HWE can be rejected")

  t1<-"Chi-square test for Hardy-Weinberg equillibrium for a bi- or triallellic system. Theory: http://www.husdyr.kvl.dk/htm/kc/popgen/genetics/2/2.htm"

  list<-list(info=t1,n.total=b,allele.dist=al_dist,observed.dist=snp_obs,expected.dist=snp_exp,chi.value=chi,p.value=p_v,df=df,interpretation=int)

  return(list)
}