mirror of
https://github.com/agdamsbo/daDoctoR.git
synced 2024-11-21 19:30:22 +01:00
102 lines
2.5 KiB
R
102 lines
2.5 KiB
R
#' 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)
|
|
}
|