#Use the output file as processed in 3.1.4, change FILEPATH to the folder where output_file2 is

#Run each command line by line unless otherwise stated

CombDat_QN = read.delim('FILEPATH/Fold Changes.txt',header=T)
CombDat_QN = as.matrix(CombDat_QN)
CombDat_QN = CombDat_QN[-c(1:4),grep('FC',colnames(CombDat_QN))]
mode(CombDat_QN) = 'numeric'

CombDat_QN = log2(CombDat_QN)

path = 'FILEPATH/'

#Please set the following variables prior to execute the code:
  ncols_multiplot = 5
  #a recommended starting value is 5, the distribution plots will be organized in 5 columns. Please adjust this number as see fit. 
  #Execute the following R code plot distributions of Log2Foldchange:
  
  multiplot_width = 300 * ncols_multiplot
  multiplot_height = 300 * ceiling((length(colnames(CombDat_QN))/ncols_multiplot))
  
############ HIGHLIGHT ALL BEFORE RUN /START ############
  
  #define multiplot function
  multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
    png(filename = paste0(path,"FC_distribPlot.png"),
        width = multiplot_width, height = multiplot_height, units = "px", pointsize = 12,
        bg = "white")
    print("plotting data distribution of each Sample", quote=FALSE)
    require(grid)
    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)
    
    numPlots = length(plots)
    
    # If layout is NULL, then use 'cols' to determine layout
    if (is.null(layout)) {
      # Make the panel
      # ncol: Number of columns of plots
      # nrow: Number of rows needed, calculated from # of cols
      layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                       ncol = cols, nrow = ceiling(numPlots/cols))
    }
    if (numPlots==1) {
      print(plots[[1]])
      
    } else {
      # Set up the page
      grid.newpage()
      pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
      
      # Make each plot, in the correct location
      for (i in 1:numPlots) {
        # Get the i,j matrix positions of the regions that contain this subplot
        matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
        
        print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                        layout.pos.col = matchidx$col))
      }
    }
    dev.off()
  }
  ############ HIGHLIGHT ALL BEFORE RUN /END ############
  
  start=1
  end=ncol(CombDat_QN)
  ncols=ncols_multiplot
  
  iterator_start=1
  iterator_end=end-start+1
  
  library(ggplot2)
  
  ############ HIGHLIGHT ALL BEFORE RUN /START ############
  plots <- list()  # new empty list
  for (i in iterator_start:iterator_end) {
    
    real_position=i+start-1
    label=colnames(CombDat_QN)[real_position]
    x=as.double(subset(CombDat_QN,select=label)) #get data vector to plot
    #x=as.double(CombDat_QN[,1])
    uplim=quantile(x,0.99,na.rm = TRUE)
    lowlim=quantile(x,0.02,na.rm = TRUE)
    by=(uplim-lowlim)/500
    mybreaks=seq(lowlim,uplim,by=by)
    df<-data.frame(label=factor(label),LogFoldChange=x) # build dataframe
    p1 = 	ggplot(df, aes(x=LogFoldChange))+
      geom_histogram(breaks=mybreaks)+
      ggtitle(label)+
      ylim(0,200)
    plots[[i]] <- p1  # add each plot into plot list
  }
  ############ HIGHLIGHT ALL BEFORE RUN /END ############
  
  multiplot(plotlist = plots, cols = ncols)