-
Notifications
You must be signed in to change notification settings - Fork 6
Development of analysis functions
#library("strataG.devel") #development version from GitHub
library("strataGdevel") #development version from GitHub
library("swfscMisc") #development version from GitHub
library("adegenet") #2.0
library("apex")
library("pegas")
library("hierfstat") #development version from GitHub
#GENIND
##LOCUS
#####each calculated per locus and per population x locus
####Target stats: num.samples, num.alleles, allelic.richness, obs.He, exp.He, theta, percent.unique.alleles
smryLoci <- summarizeLoci(msats)
smryPop <- summarizeLoci(msats, by.strata = TRUE)
####Target stats: HWE.df, HWE.p.val
#using genind object
hw.test(msats)
####Target stats: num.private.alleles
alleleFreqs <- alleleFreqs(msats, by.strata = TRUE)
by.loc <- sapply(alleleFreqs, function(loc) {
mat <- loc[, "freq", ]
rowSums(apply(mat, 1, function(r) {
result <- rep(FALSE, length(r))
if(sum(r > 0) == 1) result[r > 0] <- TRUE
result
}))
})
rownames(by.loc) <- strataNames(msats)
perLocus <- colSums(by.loc) #this has the number of alleles that are private per locus
t(by.loc) #the rows will be have the private alleles for each population by locus
####Target stats: mratio
#Sean's code using genind object
####Target stats: mratio.p.val
#To implement
####Target stats: FIS, FIT, FST
#convert from genind to loci, which is used by pegas
nancycatsloci<-genind2loci(nancycats)
#for loci
FSTloci<-Fst(nancycatsloci)
Fisloci <- FSTloci[ , c("Fis")]
Fitloci <- FSTloci[ , c("Fit")]
Fstloci <- FSTloci[ , c("Fst")]
#for pops
FSTpop<-Fst(nancycatsloci, pop = #column with pop information#)
Fispop <- FSTpop[ , c("Fis")]
Fitpop <- FSTpop[ , c("Fit")]
Fstpop <- FSTpop[ , c("Fst")]
####Target stats: Ne
#To implement
##GLOBAL ####Target stats: Chi2.p.val; and D,Fst,F'st,Gst,G'st,G"st + p.vals ovl <- overallTest(msats, nrep = 5, stat.list = statList("chi2"), quietly = TRUE)
##PAIRWISE ####Target stats: Chi2.p.val, D, Fst, F'st, Gst, G'st, G"st + p.vals pws <- pairwiseTest(msats, nrep = 5,stat.list = list(statGst), quietly = TRUE) pws[1]$result ####Target stats: shared.alleles.locus #####calculated per locus and per population pair sharedAlleles(msats) ####Target stats: Chord.dist # convert genind to hierfstat loci format. This is currently a work-around, it will change nancycatshier <- hierfstat:::.genind2hierfstat(nancycats) #genind2hierfstat is a hidden function in hierfstat chord.dist <- genet.dist(nancycatshier,diploid=TRUE,method="Dch")
#MULTIDNA
##LOCUS #####each calculated per locus and per population x locus ####Target stats: nucleotide.diversity nucleotideDiversity(g) ####Target stats: Fu.F, Fu.F.p.val FusFs(g) ####Target stats: Taj.D, Taj.D.p.val tajimasD(g) ####Target stats: num.samples, num.haps, hap.richness, pct.unique.haps, Heterozygosity smryLoci <- summary(g) smryPop <- summary(g, by.strata = TRUE) ####Target stats: nucleotide.divergence, mean.pct.within dA<-nucleotideDivergence(dloop) dA1$within ####Target stats: num.private.haps hapFreqs <- alleleFreqs(dloop.g, by.strata = TRUE) by.loc <- sapply(hapFreqs, function(loc) { mat <- loc[, "freq", ] rowSums(apply(mat, 1, function(r) { result <- rep(FALSE, length(r)) if(sum(r > 0) == 1) result[r > 0] <- TRUE result })) }) rownames(by.loc) <- strataNames(dloop.g) perLocus <- colSums(by.loc)#this has the number of haplotypes that are private per locus t(by.loc)# the rows will be have the private haplotypes for each population by locus ####Target stats: Ne #To implement
##GLOBAL ####Target stats: Chi2.p.val; and Fst, PHist + p.vals ovl <- overallTest(g, nrep = 5,stat.list = statList("chi2"), quietly = TRUE)
##PAIRWISE ####Target stats: dA, mean.pct.between dA<-nucleotideDivergence(dloop) dA1$between ####Target stats: Chi2.p.val; and Fst, PHist + p.vals pws <- pairwiseTest(g, nrep = 5, stat.list = list(statGst), quietly = TRUE) pws[1]$result ####Target stats: shared.haps.locus sharedAlleles(g)