# the input data: ASV table contain  time series of 25 vagianl communities. 
# the aim of this script is to construct correlation network for each vaginal community, 
# and extract the weight of edges from each network to a vector.
#the ouput data: feature matrix consisted by network feature vector and corresponding edge names. 
#note: the ASVs in ASV table have been ranked in descending order by their sum of counts  across 25 subjects.

setwd("../..")  #set the path to the file folder of ASV table
library(SpiecEasi)
#divide the subjects into three types
SBV=c(53,3,127,128,130,135,59,5,112,116,35,55,15,17,40) 
ABV=c(77,82,12,23,27,60)
H=c(96,7,49,52)

subjects=union(union(SBV,ABV),H)
df=read.csv("ASV.csv",head=T,sep=",")

#select the ASVs with sum of count in top 60 to construct network
Num_ASV=60

#generate edge names for network
ASV_names=names(df)[5:(Num_ASV+5)]
ASV_pair_name=data.frame()
for ( i in  1:(Num_features-1)){
    for ( j in  (i+1):Num_features){
                ASV_pair_name=rbind(ASV_pair_name,paste(ASV_names[i],"_",ASV_names[j],sep="") )
}
}

#generate weight for each edge  using  sparcc function
features=data.frame()
for (i in subjects){
      temp=subset(df,subject==i)
      otu=temp[,-c(1,2,3,4)]
      sparcc.bv <-sparcc(otu[,1:Num_features])
      cc= sparcc.bv $Cor
      features=rbind(features,cc[upper.tri(cc,diag = F)])
}    


edges=paste("edge_",1:dim(features)[2],sep="")    

#delete the zero-weight edges in 25 subjects.
flags=colSums(abs(features))>0
edges=edges[flags]
features=features[,flags]

# output the feature matrix and corresponding edge names
names(features)=edges
write.csv( ASV_pair_name[flags,],"ASV_pair_name.csv",row.names=F)
write.csv(cbind(subjects,features),"ASV_features.csv",row.names=F)