############################################################################################## ############################################################################################## ## Cryptic Biological Invasions: a General Model of Hybridization #© 2017, Claudio S. Quilodrán, Frédéric Austerlitz, Mathias Currat and Juan I. Montoya-Burgos. #Department of Genetics and Evolution. University of Geneva. Switzerland #corresponding author:juan.montoya@unige.ch ############################################################################################## ############################################################################################## #Our model describes the genetic consequences of the interbreeding between a species (N1) and a closely related taxon (N0). One of them usually represents newly arrived organisms. The community is composed of diploid parental taxa, initially differentiated by multiple independent biallelic loci and assuming normal recombination between homologous chromosomes, with case dependent interbreeding rates. ##################Breeding pairs function################## #We used an approximate hypergeometric model to decompose the possible genotype space in a discrete number of genetic classes, representing the proportions of loci coming from each parental species. This model is included in equations 4, 5 and 6 of the main text. require(combinat) ################################## ## Equation 4 cross <- function(nl, k, l, i) { sumb = 0.0; a=2*nl sumb=sapply( 0:round(nl*k), function(g){ sapply( 0:round(nl*l), function(h){ sumb + prob(nl,k, g/a)*prob(nl, l, h/a)*bC(k,l,g/a,h/a,i,nl) }) }) return(sum(sumb)) } ################################## ## Equation 5 ss = 0 prob <- function(nl,k,g) { u = k - 2*g ss = 0 ss = sapply( 0:((k-2*g)*2*nl), function(t){ s = t ss + nCm(round(nl - 2*nl * g), round( s ))*nCm(round(nl - 2 * nl * g - s), round(2*nl*k - 4 * nl * g - s)) }) ss=sum(ss) return(ss * nCm (nl, round(2* nl *g)) / nCm(2* nl, round(2* nl *k))) } ################################## ## Equation 6 bC <- function (k,l,g,h,i,nl) { uu = round(2*nl*k-4*nl*g+2*nl*l-4*nl*h) aa = round(2*nl*i-2*nl*g-2*nl*h) if (uu<0 || aa<0 || uu 1){sum(sapply(regmatches(loci[1:(position-1)], gregexpr("[a-z]", loci[1:(position-1)], perl=TRUE)), length)) } else {0} Npheno<-sapply(regmatches(genotype, gregexpr("[a-z]", genotype, perl=TRUE)), length) Pall<-if(na==0){ AA((n-position+1), (Npheno-arest)/(2*(n-position+1))) } else if(na==1){ Aa((n-position+1), (Npheno-arest)/(2*(n-position+1))) } else{ aa((n-position+1), (Npheno-arest)/(2*(n-position+1))) } return(Pall) } #The final probability is presented in equation 11 of the main text. Genotype<-function(genotype){ loci<-substring(genotype, seq(1,nchar(genotype),2), seq(2,nchar(genotype),2)) n=length(loci) Npheno<-sapply(regmatches(genotype, gregexpr("[a-z]", genotype, perl=TRUE)), length) mat<-sapply(1:n, function(i){ Ploc(genotype,i) }) return(cbind(ngen=n, class=Npheno/(2*n), freq=prod(mat) )) } #Example: the genotype “AABbCcddEe” has 5 genes. It is part of the genetic class 0.5, in which it is present at a frequency of 0.03 genotype="AABbCcddEe" Genotype(genotype)