Report about GRMs file

1) read giv file and get giv matrix

GI.giv <- read.table("/Users/xiaohanj/Documents/rproject/genetics/Genotypes/GI.giv", header = FALSE, sep = "")


# GRM_inverse and giv comparision 
ids <- unique(c(GI.giv$V1, GI.giv$V2))
GRM_giv <- matrix(0, nrow = length(ids), ncol = length(ids))

rownames(GRM_giv) <- ids
colnames(GRM_giv ) <- ids

for(i in 1:nrow(GI.giv)) {
  id1 <- GI.giv[i, "V1"]
  id2 <- GI.giv[i, "V2"]
  value <- GI.giv[i, "V3"]
  GRM_giv[as.character(id1), as.character(id2)] <- value
  GRM_giv[as.character(id2), as.character(id1)] <- value  
}

# get GRMs by inverting giv
GRMs<-solve(GRM_giv)

2) distribution of diagnoal values of GRMs

hist(diag(GRMs), 
     main = "Distribution of diagnoal values of GRMs", 
     xlab = "Diagnoal values of GRMs", 
     ylab = "Frequency", 
     col = "lightblue", 
     border = "black")

summarize diagnoal values of GRMs

summary(diag(GRMs))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8857  0.9701  1.0000  1.0000  1.0269  1.1245

3) distribution of relationship (related to upper triangle)

hist(GRMs[upper.tri(GRMs)], 
     main = "Distribution of Relationship of GRMs", 
     xlab = "Relationship Coefficient", 
     ylab = "Frequency", 
     col = "lightblue", 
     border = "black")

summarize relationship of GRMs

summary(GRMs[upper.tri(GRMs)])

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.18246 -0.04121 -0.01457  0.00000  0.01945  0.61727

4) compare GRMs with A matrix

load("Am_reduced.RData")
# Identify common IDs
common_ids <- intersect(rownames(Am_reduced), rownames(GRMs))

# Subset both matrices by common IDs
Am_reduced_match <- Am_reduced[common_ids, common_ids]
GRMs_match <- GRMs[common_ids, common_ids]

# obtain upper triangle values
Am_values_upper <- Am_reduced_match[upper.tri(Am_reduced_match, diag = TRUE)]
GRMs_values_upper <- GRMs_match[upper.tri(GRMs_match, diag = TRUE)]
Am_GRMs_upper <- data.frame(Am_value = Am_values_upper, GRMs_value = GRMs_values_upper)


# make violin graph
library(ggplot2)
ggplot(Am_GRMs_upper, aes(x = factor(Am_value), y = GRMs_value, fill = factor(Am_value))) +
  geom_violin(trim = FALSE) +
  stat_summary(fun = median, geom = "point", size = 2, color = "black") +
  labs(x = "A matrix values", 
       y = "GRMs values", 
       fill = "A value",
       title = "Violin Plot of GRMs and A matrix") +
  theme_minimal()