SEM path diagram (using igraph)

(Hao slack, Jan 7th 2021)

Saved on server: /space/chen-syn01/1/data/cinliu/plots/path_diagrams_igraph

original data from:
- /space/chen-syn01/1/data/haowang/glus_gwas/genomicSEM/*.Rdata
- /space/chen-syn01/1/data/haowang/glus_gwas/genomicSEM/results_CommonFactor_UKB_raw/*.Rdata

Goal plots:

The goal is to make a total of 9 plots using igraph.

7 plots plotting 1 factor only.
2 plots plotting 2 factors.

Sample of 1 factor plots

ABCD_EUR_DWLs.tiff:

List of the 7 plots (1 factor):

2 Factor plots: UKB_raw_area_CFA.tiff and UKB_thickness_CFA.tiff

Original Instructions

*instructions given through zoom call

Hao made SEM calculations using GenomicSEM package (https://github.com/MichelNivard/GenomicSEM/wiki/3.-Models-without-Individual-SNP-effects), and now wishes to use semPlot to make path diagram
After exploring I was unable to find way to use semPlot with our results (I think GenomicsSEM is not a package semPlots/ lavvan supports just yet; so I've made the plots in igraph for now
- (note: semPlots require the raw data and the model that was used to make the calculations. I was not able to find the formulation for the "Common Factor Model" that GenomicsSEM uses (common factor model mentioned in step 3: https://github.com/MichelNivard/GenomicSEM/wiki/3.-Models-without-Individual-SNP-effects)

Files given (data)

Files are given in .Rdata from. There were 2 sets of data:
- /space/chen-syn01/1/data/haowang/glus_gwas/genomicSEM/*.Rdata
  - 7 .Rdata files. We will produce 7 plots from this directory.
- /space/chen-syn01/1/data/haowang/glus_gwas/genomicSEM/results_CommonFactor_UKB_raw/*.Rdata
  - 4 .Rdata files. We will use it to produce 4 graphs total

Code (Rscripts)

path_diagrams_7plots.R

rm(list = ls())# igraph setwd("~/Desktop/2020lab/Hao")
library(igraph) ####VERY IMPORTANT NOTE ABOUT igraph package#####the size of the self loop can only be changed by changing the source code #https://stackoverflow.com/questions/40558019/reduce-size-of-self-loops-in-igraph
#1) use Rstudio #2) use the line of code below: #trace("plot.igraph", edit=TRUE) #3) a pop up window should occur and you should be able to change the source code#4) please look for the line of code : (it is under the loop function section )# cp <- matrix(c(x0, y0, x0 + 0.4, y0 + 0.2, x0 + 0.4, y0 - 0.2, x0, y0), ncol = 2, byrow = TRUE)##and change it to:# cp <- matrix(c(x0, y0, x0 + 3 y0 + 3 x0 + 3, y0 - 3, x0, y0), ncol = 2, byrow = TRUE)#5)save and the editing window should close #### instructions end####
library(DiagrammeR)library("networkD3")library("dplyr")library("stringr")
#file path: /space/chen-syn01/1/data/haowang/glus_gwas/genomicSEMfile_list= c("ABCD_ALL_DWLs","ABCD_EUR_DWLs", "ABCD_NonEUR_DWLs","ABCD_singletons_ALL_DWLs", "ABCD_singletons_EUR_DWLs","ABCD_singletons_NonEUR_DWLs","UKB_DWLs")
for (j in 1:length(file_list)) {file_name = file_list[j]load(paste0("~/Desktop/2020lab/Hao/genomicSEM/",file_name,".Rdata"))
r <-CommonFactor_DWLS$results
#rename labels to using abbreviation Chenpheno_abbrev <- read.delim("~/Downloads/Chenpheno_abbrev.txt") #need to double check, I made the list with process of elimination index <- match(r$lhs, Chenpheno_abbrev$Shorttened_name, nomatch = 0)r$old_lhs <- r$lhs #for safety double check r$lhs[index != 0] <- Chenpheno_abbrev$Abbreviation[index]index <- match(r$rhs, Chenpheno_abbrev$Shorttened_name, nomatch = 0)r$old_rhs <- r$rhs #for safety double check r$rhs[index != 0] <- Chenpheno_abbrev$Abbreviation[index]
#making nodesnames <- unique(r$lhs)unames <- paste0("u",names[2:25]) names <- c(names,unames)nodes <- data.frame(names)nodes_1 <- create_node_df( nodes = nodes$names, label = nodes$names, type = "lower", style = "filled", color = "aqua", #incase you want to try using "networkD3" package shape = "circle", n = 49 )
#making edgesedge <- data.frame(r$lhs, r$rhs)colnames(edge) <- c("edge_from", "edge_to")val <- format(round(r$Standardized_Est,2), nsmall = 2)val2 <- format(round(r$Standardized_SE,2), nsmall =2)#make edge label edge$label <- paste0(val,"(", round(r$Standardized_SE,2) ,")" )edge$label edge$edge_from[25:48] <- paste0("u",edge$edge_to[25:48])
##matching edge with edge label head(edge)edge_1 <- create_edge_df( from = nodes_1$id[match(edge$edge_from, nodes_1$label)], to = nodes_1$id[match(edge$edge_to, nodes_1$label)], rel = edge$label, label = edge$label )nodes_1edge_1#nodes x position nodes_1$x <- NAnodes_1$x[1:25] <- 1:25nodes_1$x[26:49] <- 2:25#nodes y position nodes_1$y <- 4nodes_1$y[26:49] <- 1#nodes position for F1nodes_1$y[1] <- 17nodes_1$x[1] <-13.5#making edge n <- nodes_1[,c(3,8:9)]e <- edge[,1:2]colnames(e) <- c("to", "from")
to <- e$to[25:48]from <- e$to[25:48]new <- data.frame(to, from)

e<- rbind(e,new)
n$x <- n$x*2a<- graph_from_data_frame(vertices = n, d= e, directed = T)
#nodes_1$inodes_1$i <- nodes_1$ynodes_1$i[2:49] <- nodes_1$y[2:49]+1.5nodes_1$i[26:49] <- nodes_1$y[26:49]-3
nodes_1$i <- 5.5nodes_1$i[2:25][c(TRUE, FALSE)] <- 6.6nodes_1$i[26:49]<- -1
#manual label position adjustments y_cord <- nodes_1$i[2:49]y_cord <- y_cord +0.5y_cord[25:48] <- -.6#y_cord[21] <- 6
x_cord <- c( 8, 8.5, 11.6, 11.8, 14.9 ,15.5 ,17.8 ,19, 21, 22, 24, 25.3, 27, 29, 30.1, 32.5, 33.6, 36.3, 36.6, 39.4 ,40.2, 43, 43.2,46.8, 4 , 6, 8, 10, 12, 14 ,16 ,18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46,48, 50)x_cord[1:24] <-x_cord[1:24] -0.5x_cord[25:49] <-x_cord[25:49] -1
#####plot making####tiff(paste0(file_name,".tiff"), units="in", width=15, height=8,res=300) par(mar=c(0, 0, 0, 0), xaxs='i', yaxs='i')plot.igraph(a, vertex.size=150, vertex.label.cex = 0.65, vertex.color="white", vertex.label.color= "black", vertex.label.family="Helvetica", rescale = FALSE, ylim=c(6,10),xlim=c(1.5,51.5), asp = 1, edge.arrow.size=0.5, #the next 4 lines of comment will make the edge labels but the edge lines will run over the numbers #edge.label=edge$label, #edge.label.y = nodes_1$i[2:49], #edge.label.x =c( 8, 8.5, 11.3, 12, 15 ,15.7 ,17.6 ,19, 21, 22, 24, 26, 28, 29.7, 31.2, 33.3, 34.8, 36.5, 37.6, 40 ,40.5, 43.5, 43.5,46.8, # 4 , 6, 8, 10, 12, 14 ,16 ,18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46,48, 50), #edge.label.cex =0.6, #edge.label.family ="Helvetica", box.col = "lightblue", edge.label.color = "dodgerblue4", edge.width =1, edge.color = "gray89", edge.loop.angle=1.5)
text(40,13, paste(file_name, "\n\n Model fitChi-square:",round(CommonFactor_DWLS$modelfit$chisq, 2), "\n", "df:", round(CommonFactor_DWLS$modelfit$df,2), "\n", "AIC:", round(CommonFactor_DWLS$modelfit$AIC,2),"\n", "CFI:", round(CommonFactor_DWLS$modelfit$CFI,2),"\n", "SRMR:", round(CommonFactor_DWLS$modelfit$SRMR,2)), pos = 4)#make edge label with background for (i in 1:49) { par(font=2) legend(x_cord[i], y_cord[i], paste(edge$label[i]), box.col = "white", bg = "white", adj = 0.3,cex=0.6, text.width = strwidth(edge$label)[1]*0.3, #adjust the size of the box behind the text text.col = "dodgerblue4")}dev.off()}

path_diagram_area_thickness_2plots.R

rm(list = ls())setwd("~/Desktop/2020lab/Hao/genomicSEM/results_CommonFactor_UKB_raw/")#/space/chen-syn01/1/data/haowang/glus_gwas/genomicSEM/results_CommonFactor_UKB_raw
library(igraph) ####VERY IMPORTANT NOTE ABOUT igraph package#####the size of the self loop can only be changed by changing the source code #https://stackoverflow.com/questions/40558019/reduce-size-of-self-loops-in-igraph
#1) use Rstudio #2) use the line of code below: #trace("plot.igraph", edit=TRUE) #3) a pop up window should occur and you should be able to change the source code#4) please look for the line of code : (it is under the loop function section )# cp <- matrix(c(x0, y0, x0 + 0.4, y0 + 0.2, x0 + 0.4, y0 - 0.2, x0, y0), ncol = 2, byrow = TRUE)##and change it to:# cp <- matrix(c(x0, y0, x0 + 3 y0 + 3, x0 + 3, y0 - 3, x0, y0), ncol = 2, byrow = TRUE)#5)save and the editing window should close #### instructions end####
library(DiagrammeR)library("networkD3")library("dplyr")library("stringr")
file_list= c("UKB_raw_area_DWLs","UKB_raw_thickness_DWLs")
for (j in 1:length(file_list)) { file_name = file_list[j]load(paste0("~/Desktop/2020lab/Hao/genomicSEM/results_CommonFactor_UKB_raw/",file_name,".Rdata"))
r <-CommonFactor_DWLS$results
#relabel Chenpheno_abbrev <- read.delim("~/Downloads/Chenpheno_abbrev.txt") #need to double check, I made the list with process of elimination #relabel lhsindex <- match(r$lhs, Chenpheno_abbrev$Shorttened_name, nomatch = 0)r$old_lhs <- r$lhs #for safety double check r$lhs[index != 0] <- Chenpheno_abbrev$Abbreviation[index]#relabel rhsindex <- match(r$rhs, Chenpheno_abbrev$Shorttened_name, nomatch = 0)r$old_rhs <- r$rhs #for safety double check r$rhs[index != 0] <- Chenpheno_abbrev$Abbreviation[index]

names <- unique(r$lhs)unames <- paste0("u",names[2:13])
names <- c(names,unames)nodes <- data.frame(names)nodes_1 <- create_node_df( nodes = nodes$names, label = nodes$names, type = "lower", style = "filled", color = "aqua", shape = "circle", n = 25 )
edge <- data.frame(r$lhs, r$rhs)colnames(edge) <- c("edge_from", "edge_to")#edge$style <- "solid"val <- format(round(r$Standardized_Est,2), nsmall = 2)val2 <- format(round(r$Standardized_SE,2), nsmall =2)
edge$label <- paste0(val,"(", round(r$Standardized_SE,2) ,")" )edge$label edge$edge_from[13:24] <- paste0("u",edge$edge_to[13:24])
##matchinghead(edge)edge_1 <- create_edge_df( from = nodes_1$id[match(edge$edge_from, nodes_1$label)], to = nodes_1$id[match(edge$edge_to, nodes_1$label)], rel = edge$label, label = edge$label )nodes_1$x <- NAnodes_1$x[1:13] <- 1:13nodes_1$x[14:25] <- 2:13
nodes_1$y <- 4nodes_1$y[14:25] <- 1
nodes_1$y[1] <- 16nodes_1$x[1] <-7
n <- nodes_1[,c(3,8:9)]e <- edge[,1:2]colnames(e) <- c("to", "from")
to <- e$to[13:24]from <- e$to[13:24]new <- data.frame(to, from)
e<- rbind(e,new)
n$x <- n$x*2a<- graph_from_data_frame(vertices = n, d= e, directed = T)#edge label position nodes_1$i <- nodes_1$ynodes_1$i[2:25] <- nodes_1$y[2:25]+1.5nodes_1$i[14:25] <- nodes_1$y[14:25]-3
nodes_1$i <- 5.5nodes_1$i[2:13][c(TRUE, FALSE)] <- 6.6nodes_1$i[14:25] <- -1
#manual label position adjustments for text with white backgroudy_cord <- nodes_1$i[2:25]y_cord[1:12] <- y_cord[1:12] +1
x_cord <- c(6 , 7, 9, 10, 12.2, 13.8 ,15.3 ,17, 18, 20, 21, 24, 4 , 6, 8, 10, 12, 14 ,16 ,18, 20, 22.1, 24.2, 26.3)x_cord[1:12] <-x_cord[1:12] -0.5x_cord[13:24] <-x_cord[13:24] -1
#plottiff(paste0(file_name,".tiff"), units="in", width=7, height=6, res=300) plot.igraph(a, #edge.label=edge$label, vertex.size=180, vertex.label.cex = 0.6, vertex.color="white", vertex.label.color= "black", vertex.label.family="Helvetica", rescale = FALSE, ylim=c(-1.5,16.5),xlim=c(3.9,26.2), asp = 0, edge.arrow.size=0.5, #edge.label.y = nodes_1$i[2:25], #edge.label.x =c( 6 , 7, 8.8, 10, 12.2, 14 ,15.8 ,18, 18.7, 21, 22.1, 25, # 4 , 6, 8, 10, 12, 14 ,16 ,18, 20, 22.1, 24.2, 26.3), #edge.label.cex =0.6, #edge.label.family ="Helvetica", #edge.label.color = "dodgerblue4", edge.width =1, edge.color = "gray89", edge.loop.angle=1.5)
text(19,13, paste(file_name, "\n\n Model fitChi-square:",round(CommonFactor_DWLS$modelfit$chisq, 2), "\n", "df:", round(CommonFactor_DWLS$modelfit$df,3), "\n", "AIC:", round(CommonFactor_DWLS$modelfit$AIC,2),"\n", "CFI:", round(CommonFactor_DWLS$modelfit$CFI,2),"\n", "SRMR:", round(CommonFactor_DWLS$modelfit$SRMR,2)), pos = 4, cex = 0.9)for (i in 1:49) { legend(x_cord[i], y_cord[i], paste(edge$label[i]), box.col = "white", bg = "white", adj = 0.3,cex=0.56, text.width = strwidth(edge$label)[1]*0.3, #adjust the size of the box behind the text text.col = "dodgerblue4")}dev.off()}

Page updated

Report abuse