R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(readr, verbose = F)
library(tidyverse, verbose=F)
mat<-read_tsv('/Users/Cyrus/Desktop/828O/Assignment1/Data File S2. Raw genetic interaction datasets_ Matrix format/SGA_NxN_clustered.cdt')

amat <- mat %>%
  # extract data rows and columns
  slice(-(1:5)) %>%
  select(starts_with("dma")) %>%
  
  # convert to numeric matrix
  type_convert() %>%
  as.matrix()
  
# turn into unweighted, undirected adjacency matrix  
amat <- 1 * (abs(amat) > 0.2)
  
# extract data about arrays (columns)
coldata <- mat %>%
  slice(1:5) %>%
  select(GID, starts_with("dma")) %>%
  slice(2) %>% select(-1) %>%
  gather(dma, orf, starts_with("dma"))
 
# extract ORF ids for queries (rows)
row_orf <- mat %>%
  slice(-(1:5)) %>%
  pull(ORF)

# match row and column ORF ids
m <- match(row_orf, coldata$orf)
rows_to_use <- !is.na(m)
cols_to_use <- m[rows_to_use]

# subset matrix into ORFs found in both rows and columns
amat <- amat[rows_to_use, cols_to_use]

# set diagonal and missing entries in matrix to 0
diag(amat) <- 0
amat[is.na(amat)] <- 0

# make the adjacency matrix diagonal
amat <- ceiling(0.5 * (amat + t(amat)))

#Removing repeated rows
Repeated_Rows_id= which(duplicated(row_orf[rows_to_use]))
Repeated_Rows_Name = row_orf[rows_to_use][which(duplicated(row_orf[rows_to_use]))]
Repeated_col_id= match(Repeated_Rows_Name, coldata$orf[cols_to_use])

#Matrix without any repeated rows or column
amat=amat[-Repeated_Rows_id, -Repeated_col_id]

# of Vertices, V
V<- ncol(amat)
# of Egdes, E 
E<- sum(amat)
#Average Degree, avDeg
avDeg<- (E/V)/2 
#Density
D= (2*avDeg)/(V-1)
print(c(Vertices=V, Edges=E, AvgDegree=avDeg, Density=D))

    Vertices        Edges    AvgDegree      Density 
2.801000e+03 6.676100e+04 1.191735e+01 8.512394e-03 

#Degree Distibution
Deg_Dist=apply(amat, 1, sum)
hist(Deg_Dist, freq=F, breaks=200)

#Log-Log plot
Deg_Dist=Deg_Dist[Deg_Dist>0]
Deg_Prob_logscale=log(table(Deg_Dist)/V)
Deg_logscale=log(as.numeric(rownames(Deg_Prob_logscale)))
plot(Deg_logscale, Deg_Prob_logscale)

#BFS is reinventing the wheel
library(spa)
#Shortest_PathNxN=floyd(amat)

#Average_Distance
mean(Shortest_PathNxN[!is.infinite(Shortest_PathNxN)], na.rm=T)

[1] 0

#Diameter of the Network
max(Shortest_PathNxN[!is.infinite(Shortest_PathNxN)], na.rm=T)

[1] 0

library(ggplot2)
Temp_Data=data.frame()
hist(Shortest_PathNxN)

#Calculating Clustering Coefficient
cluster_coefficient<-function(i){
  tryCatch({
  Neighbors_id=which(amat[i,]==1)
  All_Poss_comb= length(Neighbors_id)*(length(Neighbors_id)-1)/2
  Ngbr_Conn_count=sum(sapply(Neighbors_id, function(x) sum(amat[x,Neighbors_id])))
  data.frame(Cf=log(All_Poss_comb/Ngbr_Conn_count), Vertices=log(length(Neighbors_id)))
  },  error= function(err){ data.frame(Cf=0,Vertices=0) })
}
#Clustering Coefficient:: CC
CC=as.data.frame(t(sapply(1:ncol(amat), function(i) cluster_coefficient(i))))
CC=data.frame(Vertices=unlist(CC$Vertices), Cf=unlist(CC$Cf))
CC=do.call(data.frame,lapply(CC, function(x) replace(x, is.infinite(x),NA)))
CC=na.omit(CC)
#Plotting CC vs Vertices log-log scale
library(ggplot2)
TT=aggregate(CC, list(CC$Vertices), mean)
plot(x=TT$Vertices, y=TT$Cf, main="Vertices vs Cluster Coefficient LogLog", xlab = "Log(Vertices)", ylab="Log(Cluster_Coef")

END