############# # 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) # additive Ga <- G$Ga dim(Ga) Ga[1:7, 1:7] # dominance Gd <- G$Gd dim(Gd) Gd[1:7, 1:7] #################### # graphs analysis #################### # svd decomposition - by individuals svdGa <- svd(Ga, nu = ncol(Ga), nv = nrow(Ga)) svdGd <- svd(Gd, nu = ncol(Gd), nv = nrow(Gd)) par(mfrow = c(1,2)) plot(cumsum((svdGa$d[1: ncol(Gd)])^2/sum(svdGa$d^2)), ylab = "proportion accumulated", xlab = "number of individuals", col = "red", main = "Ga") plot(cumsum((svdGd$d[1: ncol(Gd)])^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:5] pcsGd <- Gd %*% svdGd$v dim(pcsGd) pcsGd[1:14,1:5] # PCA graphs for A and D # proportion explained by the first componentes axispcsA <- paste((round(svdGa$d[1: ncol(Ga)]^2/sum(svdGa$d^2)*100))[1:3], "%", sep = "") axispcsA axispcsD <- paste((round(svdGd$d[1: ncol(Gd)]^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 = "blue", main = "Ga = PC 1 vs PC 3") # for D plot(x = pcsGd[,1], y = pcsGd[,2], xlab = axispcsD[1], ylab = axispcsD[2], col = "red", main = "Gd = PC 1 vs PC 2") plot(x = pcsGd[,1], y = pcsGd[,3], xlab = axispcsD[1], ylab = axispcsD[3], col = "blue", 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]) # Effective size per individual (Ne.i <- 1/(2*Fi)) # Ne of population (Ne.pop <- sum(Ne.i)) # Population endogamy rate (Fst.pop <- 1/(2*Ne.pop)) #############