graphlet_correlation_distance.R

library(SpiecEasi)

data(amgut1.filt)
depths <- rowSums(amgut1.filt)
amgut1.filt.n  <- t(apply(amgut1.filt, 1, norm_to_total))
amgut1.filt.cs <- round(amgut1.filt.n * min(depths))

d <- ncol(amgut1.filt.cs)
n <- nrow(amgut1.filt.cs)

## Simulate data with random OTU size
make_data <- function(type="cluster") {
  p <- rpois(1, d)
  graph <- make_graph(type, p, p)
  Prec  <- graph2prec(graph)
  Cor   <- cov2cor(prec2cov(Prec))
  x <- amgut1.filt.cs[,sample(1:d, p, replace=TRUE)]
  synth_comm_from_counts(x, mar=2, distr='zinegbin', Sigma=Cor, n=n)
}

## two types of networks
set.seed(10010)
Xli <- c(replicate(5, make_data('cluster'),  simplify=FALSE), 
         replicate(5, make_data('band'),  simplify=FALSE))
type <- as.factor(rep(c('cluster', 'band'), each=5))

## Run SPIEC-EASI mb
se <- lapply(Xli, spiec.easi, method='mb', lambda.min.ratio=1e-1, pulsar.params=list(ncores=5))
refits <- lapply(se, function(x) x$refit$stars)


## Compute graphlet correlation vectors
library(orca) # not imported
gcmat <- t(sapply(refits, pulsar::gcvec))

## embedding and plot
mds <- cmdscale(dist(gcmat))
plot(mds, col=type)