--- Goal: "Proportion of assortative mating at MC1R" article: "Frequency of colour-related MC1R genotypes do not follow Hardy-Weinberg expectation in the barn owl (Tyto alba)" authors: "Valerie DUCRET, Arnaud GAIGHER, Celine SIMON, Jerome GOUDET and Alexandre ROULIN" date: "2015" --- ###################################################################################################### #####Script to measure the global proportion of assortative pairing and the power to detect it ######## ############################################################################################################ #Explanation : In order to estimate the potential effect of non-random pairing on the overall deviation from HW proportion, we calculated the proportion of couples that paired according to MC1R-related plumage colour, while the reminder assort randomly. Because II and VI individuals are globally similarly coloured, and hence strongly differ from VV individuals (San-Jose et al., 2015), we pooled II and VI individuals into the category “Rufous” (R) and V individuals are referred to as “White” (W). # Observation # Over the investigated period, we observed 72 Red X Red couples, 127+135 Red X White couples and 358 White X White couples. data<-matrix(c(72,127,135,358),nrow=2,byrow=T) dimnames(data)[[1]]<-dimnames(data)[[2]]<-c("R","W") data # This brings the following questions: #1. What is the proportion of assortative mating in this population? #2. Can this proportion of assortative mating explains the observed heterozygote deficiency? # Results # Proportion of assortative mating We call P the proportion of individuals with phenotype R and Q = 1 – P the proportion with phenotype W. If the population is pairing at random, we expect to find a proportion P2 of R*R couples, 2PQ of R*W and W*R couples and Q2 of W*W couples. With a fraction “x” of assortative pairing, the proportions of R*R and W*W couples are each increased by a fraction PQx, while the proportion of R*W and W*R couples decreases from 2PQ to 2PQ(1-x). In order to estimate the proportion of couples that assort according to MC1R-related plumage colour, we seek the value that minimizes the squared difference between the expected and observed number of the different types of breeding pairs. n<-sum(data) P<-(data[1,2]+data[2,1]+2*data[2,2])/2/n Q<-1-P x<-(0:1000)/1000 A<-((Q^2+P*Q*x)*n-data[1,1])^2 C<-((P^2+P*Q*x)*n-data[2,2])^2 B<-(2*P*Q*n*(1-x)-data[2,1]-data[1,2])^2 res<-A+B+C plot(x,res,type="l",xlab="Prop assortative mating",ylab="Distance from obs",col="blue",log="y") prop.assort<-x[which.min(res)] abline(v=prop.assort,col="red") # Can this proportion explain the observed heterozygote deficiency? # For this, we need to estimate the frequency of the `I` and `V` alleles. The table of adult genotypes is given below (from table S4): ad.genos<-matrix(c(28,10,171,197,493,485),ncol=2,byrow=T) dimnames(ad.genos)[[1]]<-c("II","IV","VV") dimnames(ad.genos)[[2]]<-c("males","females") ad.genos #From this, we obtain the following frequencies for `I` and `V` in males and females: ninds<-apply(ad.genos,2,sum) round(pI<-(2*ad.genos[1,]+ad.genos[2,])/ninds/2,digits=4) round(pV<-1-pI,digits=4) exp.het<-pI[1]*pV[2]+pI[2]*pV[1] # is the expected proportion of heterozygotes in the offspring if there was random mating #With the observed proportion of assortative mating, however, we expect the following proportion fo each genotypes, using data from table 1: couples<-c(2,5+7,58,21+1,106+134,358) names(couples)<-c("II/II","II/IV","IV/IV","II/VV","IV/VV","VV/VV") offsp<-matrix(c(1,0,0,.5,.5,0,.25,.5,.25,0,1,0,0,.5,.5,0,0,1),ncol=3,byrow=T) dimnames(offsp)[[1]]<-names(couples) dimnames(offsp)[[2]]<-c("II","IV","VV") cbind(couples,offsp) round(exp.offs.genos<-couples%*%offsp/sum(couples),digits=4) FIS<-1-exp.offs.genos[2]/exp.het #If the heterozygote deficiency was entirely due to assortative pairing, we should oberved an FIS of 0.05 which is not far from the observed FIS in males (0.075) but less than half the observed FIS in females (0.129). It thus looks like other factors than assortative mating influences the heterozygote deficiency, at least in females. #In order to obtain a confidence interval for the expected FIS given a proportion of assortative mating of 8.7%, we generate random draws from the expected couples distributions, with and without assortative mating: #proportion of adult genotypes prop.ad<-sweep(ad.genos,2,apply(ad.genos,2,sum),FUN="/") #expected proportion of couples under random mating rm.ad<-outer(prop.ad[,1],prop.ad[,2]) #expected proportion of couples under pure assortative mating dum1<-sweep(ad.genos[-3,],2,apply(ad.genos[-3,],2,sum),FUN="/") dum2<-outer(dum1[,1],dum1[,2]) ass.ad<-cbind(rbind(dum2,c(0,0)),c(0,0,0)) ass.ad[3,3]<-prop.ad[3,1]*prop.ad[3,2] ass.ad[-3,-3]<-ass.ad[-3,-3]*(1-ass.ad[3,3]) #frequency of the different types of couples, given prop.assort assortative matings prop.couples.ass<-(1-prop.assort)*c(rm.ad[1,1],rm.ad[1,2]+rm.ad[2,1],rm.ad[2,2],rm.ad[1,3]+rm.ad[3,1],rm.ad[2,3]+rm.ad[3,2],rm.ad[3,3])+prop.assort*c(ass.ad[1,1],ass.ad[1,2]+ass.ad[2,1],ass.ad[2,2],ass.ad[1,3]+ass.ad[3,1],ass.ad[2,3]+ass.ad[3,2],ass.ad[3,3]) #frequencies of the different types of couples, with pure random mating prop.couples.rm<-c(rm.ad[1,1],rm.ad[1,2]+rm.ad[2,1],rm.ad[2,2],rm.ad[1,3]+rm.ad[3,1],rm.ad[2,3]+rm.ad[3,2],rm.ad[3,3]) nperm<-10000 dist.fis<-numeric(nperm) nc<-sum(couples) #draws of couples assuming assortative mating dum<-rmultinom(nperm,nc,prop.couples.ass) #draws of couples assuming random mating dum1<-rmultinom(nperm,nc,prop.couples.rm) for (i in 1:nperm){ #dist.fis[i]<-1-(dum[,i]%*%offsp/nc)[2]/exp.het #no variation of couples for RM dist.fis[i]<-1-(dum[,i]%*%offsp)[2]/(dum1[,i]%*%offsp)[2] #couples are drawn both in AM and RM } hist(dist.fis,breaks=seq(-0.3,0.3,by=0.01),xlab=expression(F[IS]),ylab="",main=paste("Distribution of expected FIS \n with 8.7% of assortative mating",sep=""));abline(v=c(FIS,0.075,0.129),col=c("blue","red","red"),lwd=2); quantile(dist.fis,c(0.025,0.975)) sum(dist.fis< 0.129)/10000 sum(dist.fis> 0.129)/10000 sum(dist.fis< 0.075)/10000 sum(dist.fis> 0.075)/10000 # We see from this graph that the distribution of Fis with this proportion of assortative mating is only slightly altered. The vertical blue line of the figure corresponds to the expected Fis, given the observed couples, while the 2 vertical red lines correspond to the observed FIS in males and females. Observed male FIS belongs to the 95% CI of the distribution, while observed female FIS is larger. This confirms that other factors than assortative mating must be invoked to explain the observed heterozygote deficiency in females. # How likely are we to detect this proportion of assortative mating anyway? # The statistical power to detect assortative pairing was assessed by sampling 10,000 times from a multinomial distribution the breeding pairs according to their observed frequencies and to the proportion of assortative pairing. The 10,000 3x3 matrices of pairs of genotypes obtained this way were compared to what should be expected under random pairing by means of a chi-square test. The statistical power was estimated as the proportion of the 10,000 chi-square tests giving a p-value less than 0.05. x<-rmultinom(nperm,nc,(1-prop.assort)*rm.ad+prop.assort*ass.ad) y<-apply(x,2,function(y) chisq.test(matrix(y,nrow=3))$p.value) dum<-which(is.na(y)) pow<-sum(y<=0.05,na.rm=TRUE)/(nperm-length(dum)) # Despite this low proportion of assortative mating, we would have detected it in 70% of the cases. ###################################################################################################### #####Script to measure the proportion of assortative pairing each year ######## ############################################################################################################ #1998 as an example #Observation data<-matrix(c(6,8,8,23),nrow=2,byrow=T) dimnames(data)[[1]]<-dimnames(data)[[2]]<-c("R","W") data #Results ##Proportion of assortative mating n<-sum(data) P<-(data[1,2]+data[2,1]+2*data[2,2])/2/n Q<-1-P x<-(-1000:1000)/1000 A<-((Q^2+P*Q*x)*n-data[1,1])^2 C<-((P^2+P*Q*x)*n-data[2,2])^2 B<-(2*P*Q*n*(1-x)-data[2,1]-data[1,2])^2 res<-A+B+C plot(x,res,type="l",xlab="Prop assortative mating",ylab="Distance from obs",col="blue",log="y") prop.assort<-x[which.min(res)] abline(v=prop.assort,col="red")