<file_sep># mouse_embryonic_miRNAs
The miRNA-mRNA analysis during mouse embryonic development
<file_sep>
# reading the expression data and the cluster info
miRNAclusters = read.delim("miRNA_clusters_16.txt",header = FALSE, stringsAsFactors = FALSE)
colnames(miRNAclusters) = c("Symbol","cluster.number")

mRNAclusters = read.delim("mRNA_clusters_30.txt",stringsAsFactors = FALSE,header = FALSE) # these clusters are with stringTie ids: MSTRG#
colnames(mRNAclusters) = c("Symbol","cluster.number")

mRNAtpms <- read.table("gene.tpm.modified.txt", header=T, row.names = "Symbol")

miRNAtpms <- read.csv("CPMs.TMMnormalized.modified.csv", header=T, row.names = "Symbol")
miRNAtpms <- miRNAtpms + 0.1
miRNAtpms[,colnames(miRNAtpms) %in% paste0("t",seq(1,48))] = NA

# Reading the design matrix and sample file
edesignfile <- read.csv("miRNA.maSigPro.design.csv", header=TRUE, stringsAsFactors = F)
samples <- read.csv("miRNA_samples.all.modified.csv", header=TRUE, stringsAsFactors = F)

colors = cbind(c("forebrain","#800000"),c("midbrain","#B22222"),c("hindbrain","#DC143C"),c("neural.tube","#FF8C00"),c("cranioface","gold"),c("liver", "#9370DB"),c("heart","#31a354"),c("kidney","#0000FF"),c("lung","green2"),c("intestine","#1E90FF"),c("stomach","#40E0D0"),c("limb","#EE82EE"))
colnames(colors) = colors[1,]

#averaging the expressions for each cluster over each sample point
miRNA.tpms = matrix(0, ncol = 16,nrow = 96)  
for(i in 1:length(unique(miRNAclusters$cluster.number))){
  for(j in 1:96){
    miRNA.tpms[j,i] = median(apply(miRNAtpms[miRNAclusters[which(miRNAclusters$cluster.number == i),1],(2*j-1):(2*j)],1,function(x) mean(as.numeric(x))),na.rm = TRUE)
  }
}

mRNA.tpms = matrix(0, ncol = 30,nrow = 96)
for(i in 1:30){
  for(j in 1:96){
    mRNA.tpms[j,i] = median(apply(mRNAtpms[mRNAclusters[which(mRNAclusters$cluster.number == i),1],(2*j-1):(2*j)],1,function(x) mean(as.numeric(x))),na.rm = TRUE)
  }
}

#calculating the standard deviations for each cluster across different stages                                   
tissues = unique(edesignfile$tissue)
stdevMat = as.data.frame(matrix(0, max(miRNAclusters$cluster.number),length(tissues)),colnames = tissues)
for(i in 1:length(unique(miRNAclusters$cluster.number))){
  for(j in 1:length(tissues)){
    s = (which(samples$tissue == tissues[j]))
    s = round((s[seq(1,length(s),2)]+1)/2)
    stdevMat[i,j] = sd(miRNA.tpms[s,i],na.rm = TRUE)
  }
}

# drawing the tissue specificity heatmap
library("gplots")
png("tissueSpecificity.png", width=5, height=5, units="in", res=300)
heatmap.2(apply(stdevMat,2,function(x) as.numeric(x)),scale = "row",Rowv = FALSE, Colv = FALSE,colsep = 1:12,rowsep = 1:max(miRNAclusters$cluster.number),ColSideColors = c(colors[2,as.character(unique(tissues))]), sepwidth = c(0.02,0.05), sepcolor = "purple",na.color = "grey",margin=c(8, 4),xlab = "Tissues", ylab = "Clusters",main = "Tissue specificity matrix",col = bluered,tracecol = NA,key = FALSE,labCol = as.character(unique(tissues)))
dev.off()



# defining tissue specificity matrix
scaled_stdev = scale(t(stdevMat),center = TRUE)
scaled_stdev[scaled_stdev>=0]<- 1
scaled_stdev[scaled_stdev<0] <- 0
tissue.specificity = c()
for(i in 1:12){
  tissue.specificity = c(tissue.specificity, rep(scaled_stdev[i,],8))
}
tissue.specificity = matrix(tissue.specificity,nrow =8*12,ncol = length(unique(miRNAclusters$cluster.number)),byrow = TRUE)

# calculating the partial correlations between miRNA clusters and mRNA clusters
tissue.specificity[tissue.specificity == 1] = TRUE
tissue.specificity[is.na(tissue.specificity)] = FALSE
corMat = matrix(0, ncol = 30,nrow = length(unique(miRNAclusters$cluster.number)))
for(i in 1:length(unique(miRNAclusters$cluster.number))){
  X = miRNA.tpms[as.logical(tissue.specificity[,i]),i]
  corMat[i,] = cor(X,mRNA.tpms[as.logical(tissue.specificity[,i]),],use = "pairwise.complete.obs")
}

# drawing the heatmap for the partial correlation matrix
png("SelectiveCorrelation.png",  width=9, height=6, units="in", res=300)
heatmap.2(apply(corMat,2,function(x) as.numeric(x)),margins = c(5,5),Rowv = FALSE, Colv = FALSE,trace="none",dendrogram = 'none',density.info = "none",colsep = 1:30,rowsep = 1:length(unique(miRNAclusters$cluster.number)), sepwidth = c(0.02,0.05), sepcolor = "black",na.color = "grey",xlab = "mRNA Clusters", ylab = "miRNA Clusters",col = colorpanel(5,"#f03b20","#bcbddc","#dadaeb"),tracecol = NA,main = "Partial correlation of miRNA-mRNA clusters",labCol = paste0("gC",1:30),labRow = paste0("miC",1:length(unique(miRNAclusters$cluster.number))),key = TRUE,keysize = 1,key.par=list(mar=c(3.5,0,3,0)),key.xtickfun=function() {
  breaks <- parent.frame()$breaks
  return(list(
    at=parent.frame()$scale01(c(breaks[1],
                                breaks[length(breaks)])),
    labels=c(as.character(breaks[1]),
             as.character(breaks[length(breaks)]))
  ))
},key.xlab = "Correlation Coefficient",key.title = "Pearson Correlation")
dev.off()

save(tissue.specificity,corMat,miRNA.tpms,mRNA.tpms,file = "correlations_16miRC.RData")
<file_sep># miRNAtap compiles the predicted targets of the miRNAs for 
library(miRNAtap)
library(gplots)

#reading the mRNA and miRNA clusters
miRNAclusters = read.delim("miRNA_clusters_16.txt",header = FALSE, stringsAsFactors = FALSE)
colnames(miRNAclusters) = c("Symbol","cluster.number")

mRNAclusters = read.delim("mRNA_clusters_30.txt",stringsAsFactors = FALSE,header = FALSE)
colnames(mRNAclusters) = c("Symbol","cluster.number")
# Reading the Ensembl to ID mapping file
entrezMap = read.delim("gencode.vM10.metadata.EntrezGene",stringsAsFactors = FALSE,header = FALSE)
colnames(entrezMap) = c("Symbol","Entrez.ID")

# geneMap has the Ensambl ids for the mRNAs
geneMap = read.delim("StringTie_geneNames_CommonNames.txt",stringsAsFactors = FALSE,header = FALSE)
colnames(geneMap) = c("Symbol","Ensambl","Name")


z = merge(mRNAclusters,geneMap, by = "Symbol")
mRNAclusters = merge(z,entrezMap, by.x = 3,by.y = 1)

# targets is a list with the size of miC# x gC which holds the target list for each miCxgC interaction
targets = list()
for (i in 1:length(unique(miRNAclusters$cluster.number))){
  targets[[i]] = list()
  for (l in 1:length(unique(mRNAclusters$cluster.number))){
    targets[[i]][[l]] = data.frame(matrix(ncol = 8))
    colnames(targets[[i]][[l]]) = c("Entrez.ID", "miR.Symbol", "rank_product","rank_final","Ensambl","Symbol","cluster.number","Name")
  }
  
  # going through the miRNAs in each cluster to get the predicted targets and compile them
  for (mirna in miRNAclusters[miRNAclusters$cluster.number == i,1]) {
    # miRNAtap uses 5 different resources for target prediction
    p <- getPredictedTargets(mirna, species = 'mmu',method = 'geom',min_src = 3,synonyms = TRUE)
    if(is.null(p)){ 
      print(paste("jumping for no prediction for:",mirna))
      next
    }
    p = as.data.frame(p)
    p$rank_final =  p$rank_final/max(p$rank_final)
    r = cbind(rownames(p),mirna,p[,c("rank_product","rank_final")])
    colnames(r) = c("Entrez.ID", "miR.Symbol","rank_product","rank_final")
    
    # distributing the targets of each miRNA between different mRNA clusters
    for (j in 1:length(unique(mRNAclusters$cluster.number))){
      mRNAcluster = mRNAclusters[mRNAclusters$cluster.number == j,]
      c = merge(r,mRNAcluster,by = "Entrez.ID")
      if(!all(isEmpty(c))){
        targets[[i]][[j]] = rbind(targets[[i]][[j]], c)
      }
      
    }
    
  }
  
}
# counting the unique number of target mRNAs in each interactions
mRNAt = data.frame(matrix(ncol = length(unique(mRNAclusters$cluster.number)),nrow = length(unique(miRNAclusters$cluster.number))))
mRNAn = c()
for (i in 1:length(unique(miRNAclusters$cluster.number))){
  for (j in 1:length(unique(mRNAclusters$cluster.number))){
    mRNAn[j] = length(mRNAclusters[mRNAclusters$cluster.number == j,1]) 
    mRNAt[i,j] = length(unique(targets[[i]][[j]]$Entrez.ID))
  }
}

# calculating the contingency table and the PValues for the Chi-square test

PVals = data.frame(matrix(ncol = length(unique(mRNAclusters$cluster.number)),nrow = length(unique(miRNAclusters$cluster.number))))
for (i in 1:length(unique(miRNAclusters$cluster.number))){
  for (j in 1:length(unique(mRNAclusters$cluster.number))){
    
    contigencyTable = rbind(c(mRNAt[i,j],sum(mRNAt[i,])-mRNAt[i,j]),c(mRNAn[j]-mRNAt[i,j],sum(mRNAn[-j])-sum(mRNAt[i,])+mRNAt[i,j]))
    test = chisq.test(contigencyTable)
    PVals[i,j] = test$p.value
    
  }
}


# plotting the -log10 of the P-value 

x = apply(-log10(PVals),2,function(x) as.numeric(as.character(x)))
breaks = c(0,4.2,6,8,10,12)
png("chisq-Pvals.png", width=9, height=6, units="in", res=300)
par(mar=c(4,4,1,1))
rownames(x) = paste0("miC",1:length(unique(miRNAclusters$cluster.number)))
colnames(x) = paste0("gC",1:length(unique(mRNAclusters$cluster.number)))
heatmap.2(x,Rowv = NA,Colv = NA,trace="none",dendrogram = 'none',density.info = "none",colsep = 0:length(unique(mRNAclusters$cluster.number)),rowsep = 0:length(unique(miRNAclusters$cluster.number)),sepcolor = "black", col = colorpanel(5,"#dadaeb","#9e9ac8","#54278f"),key = TRUE,xlab = "mRNA Clusters", ylab = "miRNA Clusters",key.ylab = "P-Values(-log10)",keysize = 1,key.par=list(mar=c(3.5,0,3,0)),key.xtickfun=function() {
  breaks <- parent.frame()$breaks
  return(list(
    at=parent.frame()$scale01(c(breaks[1],
                                breaks[length(breaks)])),
    labels=c(as.character(breaks[1]),
             as.character(breaks[length(breaks)]))
  ))
},key.xlab = "P-Values",main = "Target Enrichment (Chi-sq)")
dev.off() 

save(targets,mRNAt,mRNAt,PVals,file = "predictionTargets.RData")
<file_sep># Running maSigPro on miRNAs in hpc
library(maSigPro)
library(ggplot2)


edesignfile <- read.csv("mRNA.maSigPro.design.csv", header=TRUE, stringsAsFactors = F)


## make experimental design
time <- edesignfile$time
replicates <- edesignfile$reps
forebrain <- edesignfile$forebrain
midbrain <- edesignfile$midbrain
hindbrain <- edesignfile$hindbrain
cranioface <- edesignfile$cranioface
neural.tube <- edesignfile$neural.tube
liver <- edesignfile$liver
heart <- edesignfile$heart
limb <- edesignfile$limb
stomach <- edesignfile$stomach
intestine <- edesignfile$intestine
kidney <- edesignfile$kidney
lung <- edesignfile$lung

samples <- edesignfile$Sample
ss.edesign <- cbind(time,replicates,forebrain, midbrain,
                    hindbrain, cranioface, neural.tube,
                    liver, heart
                    , limb, stomach,
                    intestine, kidney, lung)
rownames(ss.edesign) <- paste(edesignfile$Sample)
design <- make.design.matrix(ss.edesign, degree=3)



## read the miRNA file and filter the miRNAs with zero variation across the experiments
data.counts = read.delim("gene.counts.txt", row.names="Symbol")
# remove the spike ins
data.counts = data.counts[,colnames(data.counts) %in% edesignfile$Sample]
data.filtered = data.counts[apply(data.counts,1,var,na.rm = TRUE) != 0,]
data.t <- data.filtered + 0.1


datap <- p.vector(data.t, design = design, counts = FALSE) # find the significant genes
datat <- T.fit(datap, alfa = 0.01) # find the significant differences

# saving the maSigpro regressed data
save(datat,file = "mRNA_maSigPro_TMM.RData")

get<-get.siggenes(datat,rsq=0.7, vars="all") # get the list of significant genes
x <- see.genes(get$sig.genes,k=30)$cut
y <- data.frame(x)
write.table(y,"mRNA_clusters_30.txt", sep="\t", quote=FALSE, col.names=FALSE)

colors = cbind(c("forebrain","#800000"),c("midbrain","#B22222"),c("hindbrain","#DC143C"),c("neural.tube","#FF8C00"),c("cranioface","gold"),c("liver", "#9370DB"),c("heart","#31a354"),c("kidney","#0000FF"),c("lung","green2"),c("intestine","#1E90FF"),c("stomach","#40E0D0"),c("limb","#EE82EE"))
colors = t(colors)
rownames(colors) = colors[,1]
samples <- read.csv("mRNA_samples.all.modified.csv", header=TRUE, stringsAsFactors = F)
tpms <- read.delim("gene.tpm.modified.txt", header=T, row.names = "Symbol")
tpms <- tpms + 0.1

tpms[,colnames(tpms) %in% paste0("t",seq(1,48))] = NA

# plotting the miRNA cluster profiles
for (clust in 1:30){
  cluster <- rownames(y)[y$x == clust]
  cluster.data <- tpms[cluster,]
  med <- apply(cluster.data, 2, median)
  
  tissues <- c()
  medians <- c()
  std.err <- c()
  
  for (tissue in unique(edesignfile$tissue)){
    tissue.medians <- med[samples$tissue == tissue]
    list.split <- split(tissue.medians, ceiling(seq_along(tissue.medians)/2))
    f.d <- c()
    se <- c()
    for (i in list.split){
      f.d <- c(f.d, mean(i))
      se <- c(se, sd(i)/sqrt(2))
    }
    tissues <- c(tissues , rep(tissue,8))
    medians <- c(medians, f.d)
    std.err <- c(std.err, se)
  }
  Tissues <- Tissue <- factor(tissues, 
                              levels=c("forebrain","midbrain","hindbrain","neural.tube",
                                       "cranioface","liver",
                                       "heart","kidney","lung","intestine"
                                       ,"stomach","limb"))
  df <- data.frame(tissue = Tissues, time=rep(c(10.5,11.5,12.5,13.5,14.5,15.5,16.5,18.5),12),
                   medians = medians, se = std.err)
  df <- na.omit(df)
  
  clust.name <- plot_num <- paste("cluster",clust,": ",length(cluster), " genes",sep = "")
  
  p <- ggplot(df, aes(x=time, y=medians, group=tissue, color = tissue)) + 
    geom_line() +
    geom_point()+
    geom_errorbar(aes(ymin=medians-se, ymax=medians+se), width=.8,
                  position=position_dodge(0.0))+ 
    #labs(y="",x="",title = clust.name) 
    labs(y="Median profile (TPM)",x="Developmental stage",title = clust.name)
  
  p <- p + theme(text = element_text(size = 20))
  #p <- p + geom_line(aes(size = 0))
  #p <- p + theme_bw()
  
  p <- p + theme_bw()
  p <- p + theme(panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 axis.line = element_line(colour = "black"))
  
  p <- p + scale_color_manual(values = as.character(colors[unique(df$tissue),2]), guide = F)
  
  p <- p + theme(axis.text=element_text(size=20),
                 axis.title=element_text(size=20, face = "bold"),plot.title = element_text(size=25) )
  p <- p + scale_x_discrete(limit = c(10.5,11.5, 12.5, 13.5, 14.5, 15.5, 16.5,18.5),
                            labels = c("e10.5","e11.5","e12.5", "e13.5", "e14.5", "e15.5", "e16.5","P0"), expand = c(0.05, 0.05))
  print(p)
  plot_num <- paste("cluster",clust,".pdf",sep="")
  assign(paste("fig",clust,sep=""),p)
}

require(gridExtra)
grid.arrange(fig1,fig2,fig3,fig4,fig5,fig6,fig7,fig8,fig9,fig10,fig11,fig12,
             fig13,fig14,fig15,fig16,fig17,fig18,fig19,fig20,fig21,fig22,fig23,
             fig24,fig25,fig26,fig27,fig28,fig29,fig30,
             ncol=4)
dev.print(pdf, paste0('miRNA_clusters_',16,'.pdf'),width = 50, height = 60)
<file_sep>
library(ggplot2)
require(gridExtra)

#reading the gene and miRNA clusters with parsing to match the targetfinding method
miRNAclusters = read.delim("miRNA_clusters_16.txt",header = FALSE, stringsAsFactors = FALSE)
colnames(miRNAclusters) = c("Symbol","cluster.number")

mRNAclusters = read.delim("mRNA_clusters_30.txt",stringsAsFactors = FALSE,header = FALSE)
colnames(mRNAclusters) = c("Symbol","cluster.number")

mRNAtpms <- read.table("gene.tpm.modified.txt", header=T, row.names = "Symbol")

miRNAtpms <- read.csv("CPMs.TMMnormalized.modified.csv", header=T, row.names = "Symbol")
miRNAtpms <- miRNAtpms + 0.1 # pseudo counts
miRNAtpms[,colnames(miRNAtpms) %in% paste0("t",seq(1,48))] = NA  # dealing with the missing data labeled as t1 to t48

# averaging the experessions for each cluster at each of the time-stage points
miRNA.tpms = matrix(0, ncol = 16,nrow = 96)  
for(i in 1:length(unique(miRNAclusters$cluster.number))){
  for(j in 1:96){
    miRNA.tpms[j,i] = median(apply(miRNAtpms[miRNAclusters[which(miRNAclusters$cluster.number == i),1],(2*j-1):(2*j)],1,function(x) mean(as.numeric(x))),na.rm = TRUE)
  }
}

mRNA.tpms = matrix(0, ncol = 30,nrow = 96)
for(i in 1:30){
  for(j in 1:96){
    mRNA.tpms[j,i] = median(apply(mRNAtpms[mRNAclusters[which(mRNAclusters$cluster.number == i),1],(2*j-1):(2*j)],1,function(x) mean(as.numeric(x))),na.rm = TRUE)
  }
}


# loading the correlation matrix (corMat), tissue specificity (tissue.specificity) and the cluster means of mRNA and miRNA expressions (mRNA.tpms and miRNA.tpms)
load("~/Documents/MortazaviLab/Integration/FinalCodes/correlations_16miRC.RData")

#loading the target prediction matrix (PVals and targets)
load("~/Documents/MortazaviLab/Integration/FinalCodes/predictionTargets.RData")

# selecting the interaction with significant enrichment
sigInteractions = arrayInd(which(PVals < 0.05/30/length(unique(miRNAclusters$cluster.number))),c(length(unique(miRNAclusters$cluster.number)),30))

samples <- read.csv("miRNA_samples.all.modified.csv", header=TRUE, stringsAsFactors = F)
Tissues <- factor(samples$tissue[seq(1,192,2)],levels=c("forebrain","midbrain","hindbrain","neural.tube",
                                                        "cranioface","liver",
                                                        "heart","kidney","lung","intestine"
                                                        ,"stomach","limb"))
Stages <- rep(c(10.5,11.5,12.5,13.5,14.5,15.5,16.5,18.5),12)
colors = cbind(c("forebrain","#800000"),c("midbrain","#B22222"),c("hindbrain","#DC143C"),c("neural.tube","#FF8C00"),c("cranioface","gold"),c("liver", "#9370DB"),c("heart","#31a354"),c("kidney","#0000FF"),c("lung","green2"),c("intestine","#1E90FF"),c("stomach","#40E0D0"),c("limb","#EE82EE"))
colnames(colors) = colors[1,]
colors = colors[,unique(Tissues)]

# iterating through the significant interactions with negative correlation to plot the miRNA-mRNA expression side-by-side                               
for(i in 1:length(unique(miRNAclusters$cluster.number))){
  for(j in 1:length(unique(mRNAclusters$cluster.number))){
    if(PVals[i,j]<0.05/30/length(unique(miRNAclusters$cluster.number)) && corMat[i,j]<0){
      df <- data.frame(tissue = Tissues, time= Stages,miRNA = miRNA.tpms[,i],mRNA = mRNA.tpms[,j])
      df <- na.omit(df)
      p <- ggplot(df, aes(x=time, y=miRNA, group=tissue, color=tissue)) + 
        geom_line() +
        geom_line(data = df[df$tissue %in% Tissues[tissue.specificity[,i]==1],],size = 1.5) +
        scale_y_log10()+
        labs(y="Median profile (logCPM)",x="Developmental stage",title = paste0("miRNA Cluster",i," (",length(unique(miRNAclusters[miRNAclusters$cluster.number == i,1])) ," miRNAs)"))
      
      p <- p + theme(text = element_text(size = 20))
      p <- p + theme_bw()
      p <- p + scale_color_manual(values = colors[2,], guide = F)
      p <- p + theme(axis.text=element_text(size=9),
                     axis.title=element_text(size=12),plot.title = element_text(size=12) )
      p <- p + scale_x_discrete(limit = c(10.5,11.5, 12.5, 13.5, 14.5, 15.5, 16.5,18.5),
                                labels = c("e10.5","e11.5","e12.5", "e13.5", "e14.5", "e15.5", "e16.5","P0"), expand = c(0.05, 0.05))
      print(p)
      assign("mifig1",p)
      
      
      
      p <- ggplot(df, aes(x=time, y=mRNA, group=tissue, color=tissue)) + 
        geom_line() +
        geom_line(data = df[df$tissue %in% Tissues[tissue.specificity[,i]==1],],size = 1.5)+
        scale_y_log10()+
        labs(y="Median profile (logCPM)",x="Developmental stage",title = paste0("mRNA Cluster",j," (",length(unique(mRNAclusters[mRNAclusters$cluster.number == j,1])) ," genes)"))
      p <- p + theme(text = element_text(size = 20))
      p <- p + theme_bw()
      p <- p + scale_color_manual(values = colors[2,], guide = F)
      p <- p + theme(axis.text=element_text(size=9),
                     axis.title=element_text(size=12),plot.title = element_text(size=12) )
      p <- p + scale_x_discrete(limit = c(10.5,11.5, 12.5, 13.5, 14.5, 15.5, 16.5,18.5),
                                labels = c("e10.5","e11.5","e12.5", "e13.5", "e14.5", "e15.5", "e16.5","P0"), expand = c(0.05, 0.05))
      print(p)
      assign("mfig1",p)
      grid.arrange(mifig1,mfig1, ncol=2)
      dev.print(pdf, paste0('miR',i,'-','gene',j,'_log.pdf'),width = 7.5, height = 3)
    }
  }
}

######################### GO analysis
geneMap = read.delim("StringTie_geneNames_CommonNames.txt",stringsAsFactors = FALSE,header = FALSE)
colnames(geneMap) = c("Symbol","Ensambl","Name")

# writting the input files for the metascape analysis at "http://metascape.org/gp/index.html#/main/step1"
for(n in 1:dim(sigInteractions)[1]){
  targetlist = targets[[sigInteractions[n,1]]][[sigInteractions[n,2]]][-1,]
  mRNAs = unique(targetlist$Symbol)
  write.table(unique(geneMap[which(geneMap$Symbol %in% mRNAs),3]),paste0("target",sigInteractions[n,1],"-",sigInteractions[n,2],".csv"),sep = ",",quote = FALSE,row.names = FALSE,col.names = FALSE)
}

# Once the GO analysis was run on the Metascape server the result table was saved with the following format and read for enrichment analysis
# "miC#-gC#_metascape_result.xlsx"                              
library(readxl)

# using the output of the Metascape analysis to plot the enriched terms 
for(n in 1:dim(sigInteractions)[1]){
  miR = sigInteractions[n,1]
  gene = sigInteractions[n,2]
  test = read_excel(paste0(miR,"-",gene,"_metascape_result.xlsx"), sheet = 2)
  test = test[!duplicated(test$Description),]
  test = test[order(test$LogP,decreasing = TRUE),]
  test$Description = factor(test$Description,levels = test$Description)
  test = test[order(test$LogP,decreasing = FALSE),]
  test$colors[test$LogP < -4] = "#dadaeb"
  test$colors[test$LogP < -6] = "#9e9ac8"
  test$colors[test$LogP < -12] = "#54278f"
  png(paste0(miR,"-",gene,"_metascape_GO1.png"), width=9, height=6, units="in", res=300)
  ggplot(test[1:15,], aes(x=Description, y=-LogP)) +
    geom_bar(stat='identity',fill = rev(test$colors[1:15]),color="black") +
    scale_y_continuous(limits=c(0,17))+
    coord_flip()+
    theme(axis.text = element_text(size = 10),axis.title=element_text(size=12))
  dev.off()
}



