In ravinpoudel/myFunctions: Collection of custom scripts

This funtion allows to extract a graph object and create a table with nodes, based on node degree. Table also include taxonomic information as well as network attributes such as node degree (nd) or size (frequence of observation). Output table is ordered from highest to lowest node degree.

vtable.ON.nodeDegree <- function(graph_object = net.grph, selected_taxonomy = sel.tax){
  aa <- data.frame(node = V(net.grph)$name, nd = degree(net.grph))
  bb <- aa[order(aa$nd, decreasing = T), ]
  test <- sel.tax[rownames(sel.tax) %in% rownames(bb),]
  kkk <- test[rownames(bb),,drop=FALSE]
  ddd <- cbind(bb,kkk)
  return(ddd)
}

Example

library(devtools)
install_github("ravinpoudel/myFunctions")
library(myFunctions)
library(igraph)
library(magrittr)
library(tidyverse)
library(Hmisc) 
library(Matrix)
library(knitr)

# import count data with read.mothur.shared function
otu_data <- read.mothur.shared(Sys.glob("*.shared"))
otu_data [1:4, 1:3]

# upload taxanomy file using read.mothur.taxonomy function
tax <- read.mothur.taxonomy(Sys.glob("*.taxonomy"))
head(tax)

# if everything went smoothly then rownames of these two files should get you a true output
all.equal(colnames(otu_data), rownames(tax))

# Keep the OTUs with more than 100 counts
otu.table.filter <- otu_data[ , colSums(otu_data) > 100]

# Check for the decrease in the number of OTUs
dim(otu.table.filter)


# Calculate pairwise correlation between OTUs
otu.cor <- rcorr(as.matrix(otu.table.filter), type="spearman")

# Get p-value matrix
otu.pval <- forceSymmetric(otu.cor$P) # Self-correlation as NA

# Select only the taxa for the filtered OTUs by using rownames of otu.pval
sel.tax <-tax[rownames(otu.pval),,drop=FALSE]

# Sanity check 
all.equal(rownames(sel.tax), rownames(otu.pval))


# Filter the association based on p-values and level of correlations
p.yes <- otu.pval < 0.001

# Select the r values for p.yes
r.val <- otu.cor$r # select all the correlation values 
p.yes.r <- r.val * p.yes # only select correlation values based on p-value criterion 

# Select OTUs by level of correlation
p.yes.r <- abs(p.yes.r) > 0.7 # output is logical vector
p.yes.rr <- p.yes.r * r.val # use logical vector for subscripting.


# Create an adjacency matrix
adjm <- as.matrix(p.yes.rr)

# Add taxonomic information associated with adjacency matrix
colnames(adjm) <- as.vector(sel.tax$Family)
rownames(adjm) <- as.vector(sel.tax$Family)

# Create an adjacency matrix in igraph format                               
net.grph <- graph.adjacency(adjm,mode="undirected", weighted=TRUE,diag=FALSE)


# Calculate edge weight == level of correlation
edgew <- E(net.grph)$weight

# Identify isolated nodes
bad.vs <- V(net.grph)[degree(net.grph) == 0] 

# Remove isolated nodes
net.grph <- delete.vertices(net.grph, bad.vs)


# Plot the graph object
plot(net.grph,
     vertex.size=4,
     vertex.frame.color="black",
     edge.curved=F,
     edge.width=1.5,
     layout=layout.fruchterman.reingold,
     edge.color=ifelse(edgew<0,"red","blue"),
     vertex.label=NA,
     vertex.label.color="black",
     vertex.label.family="Times New Roman",
     vertex.label.font=2)

# Plot the graph object with label.
plot(net.grph,
     vertex.size=4,
     vertex.frame.color="black",
     edge.curved=F,
     edge.width=1.5,
     layout=layout.fruchterman.reingold,
     edge.color=ifelse(edgew<0,"red","blue"),
     vertex.label.color="black",
     vertex.label.family="Times New Roman",
     vertex.label.font=3,
     vertex.label.cex= 0.5)

nodes_info_basedon_nodedegree <- vtable.ON.nodeDegree(net.grph, sel.tax)
kable(nodes_info_basedon_nodedegree)