spiec_test_Dal_Acne.R

#libraries
library(readr)
library(devtools)
library(SpiecEasi)
library(phyloseq)
library(tibble)
library(igraph)
library(stringr)
library(tidyr)
library(ggplot2)
library(Matrix)
library(orca)
library(ggnetwork)
#install_github('wilkox/forcats')
library(wilkoxmisc)
library(forcats)
library(intergraph)
library(network)
library(plyr)
library(dplyr)
library(tidyverse)

#Filter OTUs from table where prevalence is below 25% of dataset (bacteria)
OTUTableDCA <- read_tsv("OTU_Table_Bacterial_Dalian_Acne.txt")
# cast OTU table
OTUTableDCA <- dcast(OTUTableDCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableDCA[1] <- NULL
# select most prevalent OTUs by occurrence frequency
OTUTableDCA <- 1*(OTUTableDCA>0)
freqDCA <- apply(OTUTableDCA, 2, mean)
freqDCA <- freqDCA[freqDCA >=0.25]
freqDCA <- as.matrix(freqDCA)
abundOTUsDCA <- rownames(freqDCA)
prevalenceDCA <- freqDCA[,1]


#calculate the abundance of prevalent OTUs (bacteria)
OTUTableDCA <- read_tsv("OTU_Table_Bacterial_Dalian_Acne.txt")
# cast OTU table
OTUTableDCA <- dcast(OTUTableDCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableDCA[1] <- NULL
N <- apply(OTUTableDCA,1,sum)
OTUTableDCA <- read_tsv("OTU_Table_Bacterial_Dalian_Acne.txt") %>% filter(OTU %in% abundOTUsDCA)
OTUTableDCA <- dcast(OTUTableDCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableDCA[1] <- NULL
abundDCA <- OTUTableDCA/N
abundDCA <- apply(abundDCA,2,mean)
abundanceDCA <- abundDCA

#Store taxonomic information of prevalent OTUs (bacteria)
OTUTableDCA <- read_tsv("OTU_Table_Bacterial_Dalian_Acne.txt") %>%
filter(OTU %in% abundOTUsDCA) %>%
select(OTU, Kingdom, Phylum, Class, Order, Family, Genus) %>%
unique()
GenusDCA <- OTUTableDCA$Genus
ClassDCA <- OTUTableDCA$Class
PhylumDCA <- OTUTableDCA$Phylum
KingdomDCA <- OTUTableDCA$Kingdom

#Prepare phyloseq OTU table with prevalent OTUs (bacteria)
OTUTableDCA <- "OTU_Table_Bacterial_Dalian_Acne.txt" %>%
read_tsv %>%
filter(OTU %in% abundOTUsDCA) %>%
select(OTU, Sample, Count) %>%
mutate(OTU = word(OTU, 1, sep = ";")) %>%
# Select only most abundant OTUs
# Delete the marked lines below to include all OTUs
#group_by(OTU) %>%                   # <- Delete
#mutate(TotalCount = sum(Count)) %>% # <- Delete
#ungroup %>%                         # <- Delete
#filter(TotalCount >= 500) %>%      # <- Delete
#select(-TotalCount) %>%             # <- Delete
spread(OTU, Count, fill = 0) %>%
as.data.frame %>%
remove_rownames %>%
column_to_rownames("Sample") %>%
as.matrix %>%
otu_table(taxa_are_rows = FALSE)

#Prepare phyloseq taxa table (bacteria)
TaxaTableDCA <- "OTU_Table_Bacterial_Dalian_Acne.txt" %>%
read_tsv %>%
filter(OTU %in% abundOTUsDCA) %>%
select(OTU, Domain = Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
mutate(OTU = word(OTU, 1, sep = ";")) %>%
unique %>%
as.data.frame %>%
remove_rownames %>%
column_to_rownames("OTU") %>%
as.matrix %>%
tax_table

PhyseqDCA <- phyloseq(OTUTableDCA, TaxaTableDCA)


#Step 1 of SPIEC-EASI: perform normalization, pseudo-count, and CLR-transformation of count data (just bacteria data)
data_clr_dca <- t(clr(OTUTableDCA+1, 1))



##Optional: subselect first 100 columns (OTUs) to test script
#data_clr_bac_small <- data_clr_bac[,1:100]


################################################
################################################
################################################
################################################
################################################


#Filter OTUs from table where prevalence is below 25% of dataset (fungi)
OTUTableFCA <- read_tsv("OTU_Table_Fungal_Dalian_Acne.txt")
# cast OTU table
OTUTableFCA <- dcast(OTUTableFCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableFCA[1] <- NULL
# select most prevalent OTUs by occurrence frequency
OTUTableFCA <- 1*(OTUTableFCA>0)
freqFCA <- apply(OTUTableFCA, 2, mean)
freqFCA <- freqFCA[freqFCA >=0.25]
freqFCA <- as.matrix(freqFCA)
abundOTUsFCA <- rownames(freqFCA)
prevalenceFCA <- freqFCA[,1]


#calculate the abundance of prevalent OTUs (fungi)
OTUTableFCA <- read_tsv("OTU_Table_Fungal_Dalian_Acne.txt")
# cast OTU table
OTUTableFCA <- dcast(OTUTableFCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableFCA[1] <- NULL
N <- apply(OTUTableFCA,1,sum)
OTUTableFCA <- read_tsv("OTU_Table_Fungal_Dalian_Acne.txt") %>% filter(OTU %in% abundOTUsFCA)
OTUTableFCA <- dcast(OTUTableFCA, Sample ~ OTU, value.var = "Count", fill = 0)
OTUTableFCA[1] <- NULL
abundFCA <- OTUTableFCA/N
abundFCA <- apply(abundFCA,2,mean)
abundanceFCA <- abundFCA


#Store taxonomic information (fungi)
OTUTableFCA <- read_tsv("OTU_Table_Fungal_Dalian_Acne.txt") %>%
filter(OTU %in% abundOTUsFCA) %>%
select(OTU, Kingdom, Phylum, Class, Order, Family, Genus) %>%
unique()
GenusFCA <- OTUTableFCA$Genus
ClassFCA <- OTUTableFCA$Class
PhylumFCA <- OTUTableFCA$Phylum
KingdomFCA <- OTUTableFCA$Kingdom

#Prepare phyloseq OTU table (fungi)
OTUTableFCA <- "OTU_Table_Fungal_Dalian_Acne.txt" %>%
read_tsv %>%
filter(OTU %in% abundOTUsFCA) %>%
select(OTU, Sample, Count) %>%
mutate(OTU = word(OTU, 1, sep = ";")) %>%
# Select only most abundant OTUs
# Delete the marked lines below to include all OTUs
#group_by(OTU) %>%                   # <- Delete
#mutate(TotalCount = sum(Count)) %>% # <- Delete
#ungroup %>%                         # <- Delete
#filter(TotalCount >= 500) %>%      # <- Delete
#select(-TotalCount) %>%             # <- Delete
spread(OTU, Count, fill = 0) %>%
as.data.frame %>%
remove_rownames %>%
column_to_rownames("Sample") %>%
as.matrix %>%
otu_table(taxa_are_rows = FALSE)

#Prepare phyloseq taxa table (fungi)
TaxaTableFCA <- "OTU_Table_Fungal_Dalian_Acne.txt" %>%
read_tsv %>%
filter(OTU %in% abundOTUsFCA) %>%
select(OTU, Domain = Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
mutate(OTU = word(OTU, 1, sep = ";")) %>%
unique %>%
as.data.frame %>%
remove_rownames %>%
column_to_rownames("OTU") %>%
as.matrix %>%
tax_table

PhyseqFCA <- phyloseq(OTUTableFCA, TaxaTableFCA)

#Step 1 of SPIEC-EASI: perform normalization, pseudo-count, and CLR-transformation of count data (fungal)
data_clr_fca <- t(clr(OTUTableFCA+1, 1))


#Combine bacterial and fungal CLR-transformed data. This combined file will be used for steps 2-3 of SPIEC-EASI pipeline
data_clr_combined_dca <- cbind(data_clr_dca,data_clr_fca)

#Combine bacterial and fungal abundance and prevalence data by combining tables from separate domains. Want to retain abundance data as percentage of community within each domain, not total relative abundance of the two domains combined
freqDCA <- c(freqDCA,freqFCA)
abundDCA <- c(abundDCA,abundFCA)
prevalenceDCA <- freqDCA
abundanceDCA <- abundDCA
GenusCombinedDCA <- c(GenusDCA,GenusFCA)
ClassCombinedDCA <- c(ClassDCA,ClassFCA)
PhylumCombinedDCA <- c(PhylumDCA,PhylumFCA)
KingdomCombinedDCA <- c(KingdomDCA,KingdomFCA)

#Step 2 of SPIEC-EASI: Inverse covariance estimation (from SPIEC-EASI github)
sparseiCov <- function(data, method, npn=FALSE, verbose=FALSE, cov.output = TRUE, ...) {
    
    if (npn) data <- huge::huge.npn(data, verbose=verbose)
    
    args <- list(...)
    
    method <- switch(method, glasso = "glasso", mb = "mb", stop("Method not supported"))
    
    if (is.null(args$lambda.min.ratio)) args$lambda.min.ratio <- 1e-3
    
    if (method %in% c("glasso")) {
        do.call(huge::huge, c(args, list(x=data, method=method, verbose=verbose,
        cov.output = cov.output)))
        
    } else if (method %in% c('mb')) {
        est <- do.call(huge::huge.mb, c(args, list(x=data, verbose=verbose)))
        est$method <- 'mb'
        est$data <- data
        est$sym  <- ifelse(!is.null(args$sym), args$sym, 'or')
        return(est)
    }
}

sparseiCov_dca <- sparseiCov(data_clr_combined_dca, method="mb", lambda.min.ratio=1e-3, verbose=TRUE, nlambda=20)


#Step 3 of SPIEC-EASI: Model selection to pick right lambda (from SPIEC-EASI github)

#' Model selection for picking the right \code{lambda} penalty.
#' This is identical to huge::huge.stars except that the subsampling loop is replaced with an mclapply function to add parallelization capabilities.
#'
#' @param est an estimate/model as produced by the sparseiCov function
#' @param criterion character string specifying criterion/method for model selection accepts 'stars' [default], 'ric', 'ebic'
#' @param stars.thresh variability threshold for stars selection
#' @param ebic.gamma tuning parameter for ebic
#' @param stars.subsample.ratio The default value 'is 10*sqrt(n)/n' when 'n>144' and '0.8' when 'n<=144', where 'n' is the sample size.
#' @param rep.num number of subsamplings when \code{criterion} = stars.
#' @param ncores number of cores to use. Need multiple processers if \code{ncores > 1}
#' @param normfun normalize internally if data should be renormalized
#' @importFrom parallel mclapply
#' @export
icov.select <- function(est, criterion = 'stars', stars.thresh = 0.05, ebic.gamma = 0.5,
stars.subsample.ratio = NULL, rep.num = 20, ncores=1, normfun=function(x) x, verbose=FALSE) {
    gcinfo(FALSE)
    if (est$cov.input) {
        message("Model selection is not available when using the covariance matrix as input.")
        class(est) = "select"
        return(est)
    }
    if (!est$cov.input) {
        if (est$method == "mb" && is.null(criterion))
        criterion = "stars"
        if (est$method == "ct" && is.null(criterion))
        criterion = "ebic"
        n = nrow(est$data)
        d = ncol(est$data)
        nlambda = length(est$lambda)
        if (criterion == "ric") {
            if (verbose) {
                message("Conducting rotation information criterion (ric) selection....")
                #        flush.console()
            }
            if (n > rep.num) {
                nr = rep.num
                r = sample(n, rep.num)
            }
            if (n <= rep.num) {
                nr = n
                r = 1:n
            }
            out = .C("RIC", X = as.double(est$data), dd = as.integer(d),
            nn = as.integer(n), r = as.integer(r), nr = as.integer(nr),
            lambda_opt = as.double(0), PACKAGE = "huge")
            est$opt.lambda = out$lambda_opt/n
            rm(out)
            gc()
            if (verbose) {
                message("done\n")
                #        flush.console()
            }
            if (verbose) {
                message("Computing the optimal graph....")
                #        flush.console()
            }
            if (est$opt.lambda > max(cor(est$data)))
            est$refit = Matrix(0, d, d)
            else {
                if (est$method == "mb")
                est$refit = huge::huge.mb(est$data, lambda = est$opt.lambda,
                sym = est$sym, idx.mat = est$idx.mat, verbose = FALSE)$path[[1]]
                if (est$method == "glasso") {
                    if (!is.null(est$cov)) {
                        tmp = huge::huge.glasso(est$data, lambda = est$opt.lambda,
                        scr = est$scr, cov.output = TRUE, verbose = FALSE)
                        est$opt.cov = tmp$cov[[1]]
                    }
                    if (is.null(est$cov))
                    tmp = huge::huge.glasso(est$data, lambda = est$opt.lambda,
                    verbose = FALSE)
                    est$refit = tmp$path[[1]]
                    est$opt.icov = tmp$icov[[1]]
                    rm(tmp)
                    gc()
                }
                if (est$method == "ct")
                est$refit = huge::huge.ct(est$data, lambda = est$opt.lambda,
                verbose = FALSE)$path[[1]]
            }
            est$opt.sparsity = sum(est$refit)/d/(d - 1)
            if (verbose) {
                cat("done\n")
                #        flush.console()
            }
        }
        if (criterion == "ebic" && est$method == "glasso") {
            if (verbose) {
                cat("Conducting extended Bayesian information criterion (ebic) selection....")
                #        flush.console()
            }
            est$ebic.score = -n * est$loglik + log(n) * est$df + 4 * ebic.gamma * log(d) * est$df
            est$opt.index = which.min(est$ebic.score)
            est$refit = est$path[[est$opt.index]]
            est$opt.icov = est$icov[[est$opt.index]]
            if (est$cov.output)
            est$opt.cov = est$cov[[est$opt.index]]
            est$opt.lambda = est$lambda[est$opt.index]
            est$opt.sparsity = est$sparsity[est$opt.index]
            if (verbose) {
                message("done\n")
                #        flush.console()
            }
        }
        if (criterion == "stars") {
            if (is.null(stars.subsample.ratio)) {
                if (n > 144)
                stars.subsample.ratio = 10 * sqrt(n)/n
                if (n <= 144)
                stars.subsample.ratio = 0.8
            }
            
            #            for (i in 1:nlambda) merge[[i]] <- Matrix(0, d, d)
            
            if (verbose) {
                mes = "Conducting Subsampling....."
                message(mes, appendLF = FALSE)
                #        cat("\n")
                #        flush.console()
            }
            #    for (i in 1:rep.num) {
            premerge <- parallel::mclapply(1:rep.num, function(i) {
                #                if (verbose) {
                #                  mes <- paste(c("Conducting Subsampling....in progress:",
                #                    floor(100 * i/rep.num), "%"), collapse = "")
                #                  cat(mes, "\r")
                #                  flush.console()
                #                }
                #                merge <- replicate(nlambda, Matrix(0, d,d))
                ind.sample = sample(c(1:n), floor(n * stars.subsample.ratio),
                replace = FALSE)
                if (est$method == "mb")
                tmp = huge::huge.mb(normfun(est$data[ind.sample, ]), lambda = est$lambda,
                scr = est$scr, idx.mat = est$idx.mat, sym = est$sym,
                verbose = FALSE)$path
                if (est$method == "ct")
                tmp = huge::huge.ct(normfun(est$data[ind.sample, ]), lambda = est$lambda,
                verbose = FALSE)$path
                if (est$method == "glasso")
                tmp = huge::huge.glasso(normfun(est$data[ind.sample, ]), lambda = est$lambda,
                scr = est$scr, verbose = FALSE)$path
                #                for (j in 1:nlambda) merge[[j]] <- merge[[j]] + tmp[[j]]
                
                rm(ind.sample)
                gc()
                return(tmp)
            }, mc.cores=ncores)
            #  }
            # merge <- lapply(merge, as.matrix)
            #      merge <- lapply(merge, simplify2array)
            #      est$merge <- lapply(1:dim(merge)[3], function(i) merge[,,i]/rep.num)
            
            merge <- Reduce(function(l1, l2) lapply(1:length(l1),
            function(i) l1[[i]] + l2[[i]]), premerge, accumulate=FALSE)
            
            if (verbose) {
                message("done")
                #        cat("\n")
                #        flush.console()
            }
            est$variability = rep(0, nlambda)
            est$merge <- vector('list', nlambda)
            for (i in 1:nlambda) {
                est$merge[[i]]  <- merge[[i]]/rep.num
                est$variability[i] <- 4 * sum(est$merge[[i]] * (1 - est$merge[[i]]))/(d * (d - 1))
            }
            est$opt.index = max(which.max(est$variability >=
            stars.thresh)[1] - 1, 1)
            est$refit = est$path[[est$opt.index]]
            est$opt.lambda = est$lambda[est$opt.index]
            est$opt.sparsity = est$sparsity[est$opt.index]
            if (est$method == "glasso") {
                est$opt.icov = est$icov[[est$opt.index]]
                if (!is.null(est$cov))
                est$opt.cov = est$cov[[est$opt.index]]
            }
        }
        est$criterion = criterion
        class(est) = "select"
        return(est)
    }
}

spiec.easi.combined_dca <- icov.select(sparseiCov_dca, rep.num = 20, ncores = 16, verbose=TRUE)

#Merge bacterial and fungal phyloseq
PhyseqCombined_DC_Acne <- merge_phyloseq(PhyseqDCA,PhyseqFCA)

# 6 Convert to igraph and network objects
ig.combined_dca <- spiec.easi.combined_dca %>%
.$refit %>%
adj2igraph(vertex.attr = list(name = taxa_names(PhyseqCombined_DC_Acne)))

vertex_attr(ig.combined_dca)
edge_attr(ig.combined_dca)

betweenness_dca <- betweenness(ig.combined_dca)
ig.combined_dca <- set_vertex_attr(ig.combined_dca,"betweenness",value=betweenness_dca)
degree_dca <- degree(ig.combined_dca)
ig.combined_dca <- set_vertex_attr(ig.combined_dca,"degree", value=degree_dca)


# set OTUs vertex attributes
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "prevalence", value=prevalenceDCA)
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "abundance", value=abundanceDCA)
GenusCombinedDCA[is.na(GenusCombinedDCA)] <- "unclassified_Genus"
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "Genus", value=GenusCombinedDCA)
ClassCombinedDCA[is.na(ClassCombinedDCA)] <- "unclassified_Class"
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "Class", value=ClassCombinedDCA)
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "Phylum", value=PhylumCombinedDCA)
ig.combined_dca <- set_vertex_attr(ig.combined_dca, "Kingdom", value=KingdomCombinedDCA)

#Extrac OTU names for later use
OTUNamescombined_DC_Acne <- PhyseqCombined_DC_Acne %>% otu_table %>% colnames

# 7 Determine correlation weights and direction (positive or negative
# correlation), and save
elist.weighted.combined_dca <- spiec.easi.combined_dca %>%
getOptBeta %>%
symBeta(mode = 'maxabs') %>%
summary %>%
# Assuming OTU numbers in i and j correlate to the order of OTUs in
# se.summer3zw, restore OTU names
left_join(data.frame(OTU1 = OTUNamescombined_DC_Acne, i = 1:length(OTUNamescombined_DC_Acne))) %>%
left_join(data.frame(OTU2 = OTUNamescombined_DC_Acne, j = 1:length(OTUNamescombined_DC_Acne)))

elist.weighted.combined_dca <- elist.weighted.combined_dca %>%
mutate(CorSign = ifelse(x > 0, "Positive", "Negative")) %>%
mutate(CorMagnitude = abs(x))

#Save elist document
#write.tidy(elist.weighted.combined_dca,"elist.weighted.combined_dca.txt")

ig.combined_dca <- set_edge_attr(ig.combined_dca, "CorrelationSign", value=elist.weighted.combined_dca$CorSign)
ig.combined_dca <- set_edge_attr(ig.combined_dca, "CorrelationMagnitude", value=elist.weighted.combined_dca$CorMagnitude)


# remove weak edges (< 0.15)
Weak = which(E(ig.combined_dca)$CorrelationMagnitude<0.15)
ig.combined_dca = delete_edges(ig.combined_dca, Weak) 

# remove vertices without strong edges (independent vertices)
Isolated = which(degree(ig.combined_dca)==0)
ig.combined_dca = delete_vertices(ig.combined_dca, Isolated)

# plot the network
library(RColorBrewer)
Cor <- E(ig.combined_dca)$CorrelationMagnitude


par(mfrow=c(1,2))
coords <- layout.fruchterman.reingold(ig.combined_dca)
E(ig.combined_dca)$color[E(ig.combined_dca)$CorrelationSign == "Negative"] <- 'red'
E(ig.combined_dca)$color[E(ig.combined_dca)$CorrelationSign == "Positive"] <- 'blue'


pal <- brewer.pal(length(unique(ig.combined_dca$Kingdom)), "Set3")

pal2 <- c("#ff7f00", "#c41185",  "#320956", "#b15928","black")

fine = 5000 # this will adjust the resolving power.
plot(ig.combined_dca,
vertex.label=NA,
vertex.size= sqrt(V(ig.combined_dca)$abundance*500)+2,
vertex.color= pal2[as.numeric(as.factor(vertex_attr(ig.combined_dca, "Kingdom")))],
edge.width=Cor*4,
#layout_nicely(ig.summer),
layout=coords)
legend("topleft", legend=levels(as.factor(V(ig.combined_dca)$Kingdom)), bg= "#e1e5de",col = pal2, bty="o",pch=20, pt.cex = 1.3, cex = 0.6, text.col=pal2 , horiz = FALSE, inset = c(0.3, 0.87))
title("Dalian Cheek - Acne", line=-1)

# Perform fast-greedy community detection on network graph
kc = fastgreedy.community(ig.combined_dca)

# Perform edge-betweenness community detection on network graph
gc = edge.betweenness.community(ig.combined_dca)

# Plot community networks determined by fast-greedy and edge-betweenness methods side-by-side
# color the vertex by class rank
#pal <- brewer.pal(length(unique(V(ig.summer)$Class)), "Paired")
pal <- c("#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6",  "#ffff99", "#b15928","grey", "black")
my_color =pal[as.numeric(as.factor(V(ig.combined)$ClassCombinedDCA))]
pal2 <- c("red", "grey")
my_color2 =pal2[as.numeric(as.factor(E(ig.combined)$CorrelationSignDCA))]


par(mfrow = c(1, 2))
plot(kc, ig.combined_dca,
vertex.label = NA,
vertex.size = 5,
mark.groups=kc,
edge.color=my_color2,
col=my_color,
#membership=NULL,
#vertex.color= my_color,
main = "Fast-greedy \nmodularity=0.6 groups=8",
edge.width=Cor*5,
layout=coords)
plot(gc, ig.combined_dca,
vertex.label =NA,
vertex.size = 5,
edge.color=my_color2,
col=my_color,
#vertex.color= pal[as.numeric(as.factor(vertex_attr(ig.summer, "Class")))],
main = "Edge-betweenness \nmodularity=0.59 groups=9",
edge.width=Cor*5,
#layout_nicely(ig.bac.small),
layout=coords)




##############draw degree distribution

dd.BCH <- degree_distribution(ig.combined_bch, cumulative=FALSE)
sum(seq_along(dd.BCH)*dd.BCH)-1

dd.BCA <- degree_distribution(ig.combined_bca, cumulative=FALSE)
sum(seq_along(dd.BCA)*dd.BCA)-1

dd.DCH <- degree_distribution(ig.combined_dch, cumulative=FALSE)
sum(seq_along(dd.DCH)*dd.DCH)-1

dd.DCA <- degree_distribution(ig.combined_dca, cumulative=FALSE)
sum(seq_along(dd.DCA)*dd.DCA)-1

pdf("degree_distribution_cheek.pdf",width=5,height=5)
par(mar=c(4,4,4,4))
plot(seq_along(dd.BCH)-1, dd.BCH, type='b', xlim=c(0,30), ylim=c(0,0.2),
ylab="Frequency", xlab="Degree", col='blue')
points(seq_along(dd.BCA)-1, dd.BCA, type='b', col='green')
points(seq_along(dd.DCH)-1, dd.DCH, type='b', col='red')
points(seq_along(dd.DCA)-1, dd.DCA, type='b', col='black')
legend("topright", c("Baoding Healthy","Baoding Acne","Dalian Healthy", "Dalian Acne"), col=c("blue","green","red","black"),pch=1, lty=1, cex=0.6)
title(main="Cheek")
dev.off()
      
#Network stability estimation by natural connectivity
natcon <- function(ig) {
  N <- vcount(ig)
  adj <- get.adjacency(ig)
  evals <- eigen(adj)$value
  nc <- log(mean(exp(evals)))
  nc / (N - log(N))
}
nc.attack <- function(ig) {
  hubord <- -order(rank(betweenness(ig)), rank(degree(ig)), decreasing=TRUE)
  sapply(1:round(vcount(ig)*.95), function(i) {
    ind <- hubord[1:i]
    tmp <- delete_vertices(ig, V(ig)$name[ind])
    natcon(tmp)
  })
}


natcon <- function(ig) {
  N <- vcount(ig)
  adj <- get.adjacency(ig)
  evals <- eigen(adj)$value
  nc <- log(mean(exp(evals)))
  nc / (N - log(N))
}
nc.attack <- function(ig) {
  hubord <- order(rank(degree(ig)), decreasing=TRUE)
  sapply(1:round(vcount(ig)*.95), function(i) {
    ind <- hubord[1:i]
    tmp <- delete_vertices(ig, V(ig)$name[ind])
    natcon(tmp)
  })
}


nc.DCA <- nc.attack(ig.combined_dca)

##Plot Natural connectivity

png("connectivity_cheek.png",width=4,height=4,units="in",res=1200)
plot(seq(0,0.8,len=length(nc.Winter)),nc.Winter, type='l', col='blue', ylim=c(0,max(nc.Winter)), xlab="Proportion of Removed Nodes", ylab="Natural Connectivity",font.lab=2)
points(seq(0,0.8,len=length(nc.Spring)),nc.Spring,type='l', col='green')
points(seq(0,0.8,len=length(nc.Summer)),nc.Summer,type='l', col='red')
points(seq(0,0.8,len=length(nc.DCA)),nc.DCA,type='l')
dev.off()