Update Network_analysis_with_subset.r

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()    
+}