--- title: "A.get Microsatellite Analysis 2018" date: "Oct 9 2018" output: html_notebook --- ```{r setup, include=FALSE} #knitr::opts_chunk$set(echo = TRUE) ``` ```{r} #install_github("thibautjombart/adegenet") library(adegenet) library(hierfstat) library(pegas) squirrel.data <- read.csv("Aget_AlleleReport_All_0.csv") write.struct(squirrel.data, fname = "Aget_AlleleReport_All_0.str") squirrel.data <- read.structure("Aget_AlleleReport_All_0.str") #256 genotypes, 16 markers, 0, 1, enter, 0, y summary(squirrel.data) ``` Choose K ```{r} squirrel.aic <- snapclust.choose.k(max = 10, squirrel.data) plot(squirrel.aic) ``` Snapclust ```{r} squirrel.k5 <- snapclust(squirrel.data, 5) compoplot(squirrel.k5) ``` PCA ```{r} squirrel.pca.data <- scaleGen(squirrel.data, NA.method = "mean") squirrel.pca.plot <- dudi.pca(squirrel.pca.data, cent = FALSE, scale = FALSE) #10 axes s.class(squirrel.pca.plot$li, pop(squirrel.data), sub = "PCA - PC 1 and 2", xax = 1, yax = 2, csub = 2) s.class(squirrel.pca.plot$li, pop(squirrel.data), sub = "PCA - PC 1 and 3", xax = 1, yax = 3, csub = 2) s.class(squirrel.pca.plot$li, pop(squirrel.data), sub = "PCA - PC 2 and 3", xax = 2, yax = 3, csub = 2) #### # library(ggplot2) # library(ggfortify) # pca <- prcomp(data.frame(squirrel.pca.data)) # autoplot(pca, data = squirrel.pca.data) ``` Hardy-Weinberg ```{r} library(devtools) library(pegas) library(adegenet) library(gstudio) library(ggplot2) microdata <- read_population("Aget_AlleleReport_All_NA.csv", type = "column", locus.columns = 2:33, header = TRUE) #Allelic Diversity diversity <- genetic_diversity(microdata, mode = "A") diversity write.csv(diversity, file = "diversity.csv") freqs.strata <- frequencies(microdata) freqs.strata write.csv(freqs.strata, file = "allele freqs.csv") #Heterozygosity Ho <- Ho(microdata) He <- He(microdata) mean(Ho$Ho) mean(He$He) #Test for HW #microdata2 <- read.structure("Aget_AlleleReport_All_0.str", n.ind = 265, n.loc = 20, onerowperind = TRUE, col.lab = 1, col.pop = 0, row.marknames = 1) squirrel.dataHWE <- as.loci(squirrel.data) hw.test(squirrel.dataHWE) inbreeding <- inbreeding(squirrel.data, pop = NULL, truenames = TRUE) F <- sapply(inbreeding, mean) hist(F, main = "Average Inbreeding") which(F>0.4) which(F>0.9) max(F) #0.9028046 for ID 211 F <- inbreeding(squirrel.data, pop = NULL, truenames = TRUE, res.type = "function")[which(F>0.902)] #to get max F value instead of which() and max() F plot(F$'211', main = paste("Inbreeding of Individual", names(F)), xlab = "Inbreeding (F)", ylab = "Probability Density") #see http://adegenet.r-forge.r-project.org/files/tutorial-basics.pdf for calculating probability densities of inbred individuals ``` Missing Data ```{r} library(poppr) #find # loci containing missing values > 10% squirreldata.locina <- missingno(squirrel.data, type = "loci", cutoff = 0.1, quiet = TRUE, freq = FALSE) View(squirreldata.locina) #find # genotypes containing missing values > 10% squirreldata.genona <- missingno(squirrel.data, type - "geno", cutoff = 0.1, quiet = TRUE, freq = FALSE) View(squirreldata.genona) #replace all NA with 0 squirreldata.0 <- missingno(squirrel.data, type = "0", cutoff = 0.1, quiet= TRUE, freq = FALSE) View(squirreldata.0) ```