############# # Kinship ############# # read marker data M <- readRDS("M") dim(M) head(M[,1:14]) tail(M[,1:14]) source("https://bioconductor.org/biocLite.R") biocLite("impute") devtools::install_github(repo = 'italo-granato/snpReady', ref = 'dev') library(snpReady) # creates the K matrix - genomic relatioship matrix args(G.matrix) G <- G.matrix(M, method = "VanRaden", format = "wide", plot = TRUE) Ga <- G$Ga dim(Ga) Ga[1:7, 1:7] Gd <- G$Gd dim(Gd) Gd[1:7, 1:7] #################### # graphs analysis #################### # svd decomposition - by individuals svdGa <- svd(Ga, nu = 54, nv = 54) svdGd <- svd(Gd, nu = 54, nv = 54) par(mfrow = c(1,2)) plot(cumsum((svdGa$d[1:54])^2/sum(svdGa$d^2)), ylab = "proportion accumulated", xlab = "number of individuals", col = "red", main = "Ga") plot(cumsum((svdGd$d[1:54])^2/sum(svdGd$d^2)), ylab = "proportion accumulated", xlab = "number of individuals", col = "blue", main = "Gd") # obtainig the eigenvectors and eigenvalues for A and D pcsGa <- Ga %*% svdGa$v dim(pcsGa) pcsGa[1:14,1:6] pcsGd <- Gd %*% svdGd$v dim(pcsGd) pcsGd[1:14,1:6] # PCA graphs for A and D # proportion explained by the first componentes axispcsA <- paste((round(svdGa$d[1:54]^2/sum(svdGa$d^2)*100))[1:3], "%", sep = "") axispcsA axispcsD <- paste((round(svdGd$d[1:54]^2/sum(svdGd$d^2)*100))[1:3], "%", sep = "") axispcsD # 3D graph library(scatterplot3d) #for A par(mfrow = c(1,2)) scatterplot3d(pcsGa[,1], pcsGa[,2], pcsGa[,3], xlab = axispcsA[1], ylab = axispcsA[2], zlab = axispcsA[3], axis = TRUE, color = "red", highlight.3d = FALSE, box = TRUE, angle = 50, main = "Ga") #for D scatterplot3d(pcsGd[,1], pcsGd[,2], pcsGd[,3], xlab = axispcsD[1], ylab = axispcsD[2], zlab = axispcsD[3], axis = TRUE, color = "blue", highlight.3d = FALSE, box = TRUE, angle = 50, main = "Gd") # 2D graph par(mfrow = c(1,2)) # for A plot(x = pcsGa[,1], y = pcsGa[,2], xlab = axispcsA[1], ylab = axispcsA[2], col = "red", main = "Ga = PC 1 vs PC 2") plot(x = pcsGa[,1], y = pcsGa[,3], xlab = axispcsA[1], ylab = axispcsA[3], col = "red", main = "Ga = PC 1 vs PC 3") # for D plot(x = pcsGd[,1], y = pcsGd[,2], xlab = axispcsD[1], ylab = axispcsD[2], col = "blue", main = "Gd = PC 1 vs PC 2") plot(x = pcsGd[,1], y = pcsGd[,3], xlab = axispcsD[1], ylab = axispcsD[3], col = "red", main = "Gd = PC 1 vs PC 3") dev.off() ############## # heatmaps ############## install.packages("superheat") library(superheat) superheat(Ga, pretty.order.rows = T, pretty.order.cols = T, col.dendrogram = T, clustering.method = "kmeans", dist.method = "euclidean", bottom.label.text.size = 2, left.label.text.size = 2, legend.text.size = 5, title = "Ga") superheat(Gd, pretty.order.rows = T, pretty.order.cols = T, col.dendrogram = T, clustering.method = "kmeans", dist.method = "euclidean", bottom.label.text.size = 2, left.label.text.size = 2, legend.text.size = 5, title = "Gd") ################################## ## Endogamy and Ne ################################## # endogamy per individual - considering just the parents Fi <- round(diag(Ga)-1, 2)[1:14] Fi # Effective size per individual Ne.i <- 1/(2*Fi) Ne.i # Ne of population Ne.pop <- sum(Fi) Ne.pop # Population endogamy rate F.pop <- 1/(2*Ne.pop) F.pop