From 010b2a490f7ff14088d9a50bfe35ea1868829942 Mon Sep 17 00:00:00 2001 From: Kathryn Atherton <40177040+k-atherton@users.noreply.github.com> Date: Sat, 30 Jan 2021 19:46:43 -0500 Subject: [PATCH] Update Network_analysis_with_subset.r Subset by Site ID and loop through all sites. --- Network_analysis_with_subset.r | 115 +++++++++++++++++---------------- 1 file changed, 61 insertions(+), 54 deletions(-) diff --git a/Network_analysis_with_subset.r b/Network_analysis_with_subset.r index 736827d..6348190 100644 --- a/Network_analysis_with_subset.r +++ b/Network_analysis_with_subset.r @@ -28,68 +28,75 @@ minor.phylum.list <- rownames(minor.phylum) for (i in 1:length (minor.phylum.list)){ taxonomy$Phylum <- gsub(minor.phylum.list[i],"Others",taxonomy$Phylum)} -# Subset +# Add experimental design data for subsetting bind <- cbind (percent.t,DESIGN) -subset <- subset(bind, bind$Factor1=="AAA") # Change -subset <- subset(subset, subset$Factor2=="aaa") # Change -subset <- subset[,1:(ncol(subset)-ncol(DESIGN))] +# Subset by site +# If we get >4 samples per plot, we can change this to unique(bind$Plot40) +# However, the network package requires at least 4 samples to calculate the network, and we have 3-4 per plot +sites <- unique(bind$Site) -# Filter to pick up parts of ASVs -subset.t.filter <- subset[ ,colMeans(percent.t) >= 0.05] # To pick up >0.05% ASVs -subset.t.filter <- subset.t.filter[,colSums(subset.t.filter)>0] -print(c(ncol(subset),"versus",ncol(subset.t.filter))) +# Loop through all sites +for(j in 1:length(sites){ + subset <- subset(bind, bind$Site==sites[j]) + subset <- subset[,1:(ncol(subset)-ncol(DESIGN))] -# Calculate network -percent.cor <- rcorr(as.matrix(subset.t.filter), type="spearman") -percent.pval <- forceSymmetric(percent.cor$P) # Self-correlation as NA -#Select only the taxa for the filtered ASVs by using rownames of percent.pval -sel.tax <- taxonomy[rownames(percent.pval),,drop=FALSE] -#Sanity check --> should be "[1] TRUE" -all.equal(rownames(sel.tax), rownames(percent.pval)) + # Filter to pick up parts of ASVs + subset.t.filter <- subset[ ,colMeans(percent.t) >= 0.05] # To pick up >0.05% ASVs + subset.t.filter <- subset.t.filter[,colSums(subset.t.filter)>0] + print(c(ncol(subset),"versus",ncol(subset.t.filter))) -p.yes <- percent.cor$P<0.05 -r.yes <- percent.cor$r>0 -r.high <- percent.cor$r>0.7 -r.val <- percent.cor$r # select all the correlation values -p.r.yes = p.yes*r.yes*r.val*r.high -adjm<-p.r.yes + # Calculate network + percent.cor <- rcorr(as.matrix(subset.t.filter), type="spearman") + percent.pval <- forceSymmetric(percent.cor$P) # Self-correlation as NA + #Select only the taxa for the filtered ASVs by using rownames of percent.pval + sel.tax <- taxonomy[rownames(percent.pval),,drop=FALSE] + #Sanity check --> should be "[1] TRUE" + all.equal(rownames(sel.tax), rownames(percent.pval)) -net.grph=graph.adjacency(adjm,mode="undirected",weighted=TRUE,diag=FALSE) + p.yes <- percent.cor$P<0.05 + r.yes <- percent.cor$r>0 + r.high <- percent.cor$r>0.7 + r.val <- percent.cor$r # select all the correlation values + p.r.yes = p.yes*r.yes*r.val*r.high + adjm<-p.r.yes -# hs <- hub_score(net.grph, weights=NA)$vector　 -# as <- authority_score(net.grph, weights=NA)$vector -# pr <- page.rank(net.grph,directed=F)$vector　 -deg <- degree(net.grph, mode="all")　 + net.grph=graph.adjacency(adjm,mode="undirected",weighted=TRUE,diag=FALSE) -# Align taxonomy names -sel.tax$Phylum <- factor(sel.tax$Phylum) -others.n <- which(levels(sel.tax$Phylum)=="Others") -levels <- NULL -if (others.n == 1){ - for (i in 1:(length(levels(sel.tax$Phylum)))){ - levels <- c(levels, levels(sel.tax$Phylum)[i])} - levels <- c(levels, "Others")}else{ -for (i in 1:(others.n-1)){ - levels <- c(levels, levels(sel.tax$Phylum)[i])} -for (i in (others.n+1):(length(levels(sel.tax$Phylum)))){ - levels <- c(levels, levels(sel.tax$Phylum)[i])} -levels <- c(levels, "Others") -levels(sel.tax$Phylum) <- levels} + # hs <- hub_score(net.grph, weights=NA)$vector　 + # as <- authority_score(net.grph, weights=NA)$vector + # pr <- page.rank(net.grph,directed=F)$vector　 + deg <- degree(net.grph, mode="all")　 -# Illustrate -col=rainbow(length(levels(sel.tax$Phylum))) -plot.igraph(net.grph, vertex.size=deg*0.15,vertex.label=NA, vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.kamada.kawai) -# plot(net.grph, vertex.size=deg*0.15,vertex.label=NA,vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.random) -# plot(net.grph, vertex.size=deg*0.15,vertex.label=NA,vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.fruchterman.reingold) -legend(x = "bottomleft", legend = levels(sel.tax$Phylum), pch = 19, col = col, bty = "n", pt.cex=2, title = "Color legend") + # Align taxonomy names + sel.tax$Phylum <- factor(sel.tax$Phylum) + others.n <- which(levels(sel.tax$Phylum)=="Others") + levels <- NULL + if (others.n == 1){ + for (i in 1:(length(levels(sel.tax$Phylum)))){ + levels <- c(levels, levels(sel.tax$Phylum)[i])} + levels <- c(levels, "Others")}else{ + for (i in 1:(others.n-1)){ + levels <- c(levels, levels(sel.tax$Phylum)[i])} + for (i in (others.n+1):(length(levels(sel.tax$Phylum)))){ + levels <- c(levels, levels(sel.tax$Phylum)[i])} + levels <- c(levels, "Others") + levels(sel.tax$Phylum) <- levels} -gsize <- gsize(net.grph) -edge_density <- round(edge_density(net.grph),digit=5) -text(x=1,y=-1,paste("The number of edge = ", gsize)) -text(x=1,y=-1.1,paste("edge density = ", edge_density)) -title("XXXXXX") #Change + # Illustrate + col=rainbow(length(levels(sel.tax$Phylum))) + plot.igraph(net.grph, vertex.size=deg*0.15,vertex.label=NA, vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.kamada.kawai) + # plot(net.grph, vertex.size=deg*0.15,vertex.label=NA,vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.random) + # plot(net.grph, vertex.size=deg*0.15,vertex.label=NA,vertex.color=col[unclass(sel.tax$Phylum)],layout=layout.fruchterman.reingold) + legend(x = "bottomleft", legend = levels(sel.tax$Phylum), pch = 19, col = col, bty = "n", pt.cex=2, title = "Color legend") -# Save -dev.copy(pdf, file="~/R/Analysis/1_Test/ITS/Network.pdf", height=5, width=10) -dev.off() + gsize <- gsize(net.grph) + edge_density <- round(edge_density(net.grph),digit=5) + text(x=1,y=-1,paste("The number of edge = ", gsize)) + text(x=1,y=-1.1,paste("edge density = ", edge_density)) + title(paste("Site", sites[j])) + + # Save + dev.copy(pdf, file=paste0("~/R/Analysis/1_Test/ITS/Network_",sites[j],".pdf"), height=5, width=10) + dev.off() +}