*** network aging simulation, net-aging-sim-2015Sep15.R

helen:ms02GINPPI hqin$ ll net-aging-sim-2014Sep*

-rw-r--r--  1 hqin  staff   4.2K Sep 14 18:37 net-aging-sim-2015Sep15.R

-rw-r--r--  1 hqin  staff   6.1K Sep 14 13:17 net-aging-sim-2015Sep7b.R

net-aging-sim-2015Sep15.R (aka net-aging-sim-20150915.R)

This script now takes argument for start and end ms02 networks, lambda, p and popSize

# R -f file --args start end lambda p popSize

#  R -f net-aging-sim-2014Sep15.R --args 1 1 0.002 0.95 10

# the network processing working inside of this work should be move outside of this script. 
# This script can iterate over ms02 networks with the same lambda, p, and popSize

#20150915Mon add argument for lambda, p, popSize

# 2014 Feb 17, change name "20131221.DIPandGIN.sim.aging_v2.R" to "net-aging-sim-2014Feb17.R"

# 2013 Dec 20, merge DIP PPI and Genetic Inxt Net -> Multi-net approach

rm(list=ls())

source("lifespan.r")

source("network.r")

#for debug

#  start =1; end = 1; lambda= 1/100; p = 0.95; popSize=5;

# R -f file --args start end lambda p popSize

#  R -f net-aging-sim-2014Sep15.R --args 1 1 0.002 0.95 10

options(echo=TRUE) # if you want see commands in output file

args <- commandArgs(trailingOnly = TRUE)

print(args)

# trailingOnly=TRUE means that only your arguments are returned, check:

# print(commandsArgs(trailingOnly=FALSE))

start = as.integer(args[1]); start; 

end = as.integer(args[2]); end; 

lambda = as.numeric(args[3]); lambda;

p = as.numeric(args[4]); p;

popSize = as.numeric(args[5]); popSize;

myhost = 'greenfield'  # 'byte' 'blacklight' 'mactower'

#myhost = 'byte'  # 'byte' 'blacklight' 'mactower'

myhost = 'helen';

mydir = "/crucible/mc48o9p/hqin2/mactower-network-failure-simulation-master/ms02GINPPI"

if (myhost == 'byte') {  mydir = "/Users/hqin/github/mactower-network-failure-simulation/ms02GINPPI"

} else if (myhost == 'helen') { mydir = "/Users/hqin/github/mactower-network-failure-simulation/ms02GINPPI";  

} 

mydir

setwd(mydir)

debug = 0; 

#essential gene info

essenTb = read.csv('SummaryRegressionHetHom2013Oct29.csv', colClasses=rep('character', 9))

for ( i in start:end) {

  path = paste('dipgin.ms02.output/', i, sep='')

  ms02file = paste('ms02_', i, ".tab", sep='')

  infile = paste( path, '/', ms02file, sep=""); print(infile)

  pairs = read.table( infile, header=T, sep="\t", colClass = c("character", "character", NA) )

  print(head(pairs))

  if(debug==9) {     pairs = pairs[1:1000,]  }

  pairs = pairs[ pairs$id1 != pairs$id2, ]  

  # How do the two data set overlap? DIP seems to contain some questionable orfs

  uniq.orf.from.pairs = unique(c(pairs$id1, pairs$id2)) #5496 ORF

  matches = uniq.orf.from.pairs %in% unique(essenTb$orf)

  table(matches)

  

  # label essential nodes, remove nonesse-nonessen pairs

  essentialORFs = essenTb$orf[essenTb$essenflag=='essential']

  pairs$essen1 = pairs$id1 %in% essentialORFs

  pairs$essen2 = pairs$id2 %in% essentialORFs

  #remove nonessen <-> nonessen intxn because they do not affect aging. 

  pairs$remove = ifelse( pairs$essen1==F & pairs$essen2==F, T, F  )

  pairs= pairs[! pairs$remove, ]  

  # 34052 for one ms02 network

  

  #how many essen <--> essen intxn? 

  pairs$inxnEE = pairs$essen1 & pairs$essen2

  table(pairs$inxnEE)

  #How many essen genes? 

  tmp = essentialORFs %in% unique(c(pairs$id1, pairs$id2)) 

  table(tmp)

  essentialORFsPPI = essentialORFs[tmp] #?????

  

  #get connectivities per node

  degreeTb = data.frame( table(c(pairs[,1], pairs[,2])))

  summary(degreeTb)

  degreeTb$ORF = as.character( degreeTb[,1])

  

  degreeCutoff = 4; #######!!!!!! degree=5 is the median

  tmp = essentialORFsPPI %in% degreeTb$ORF[degreeTb$Freq>degreeCutoff]

  GooddEssentialORFsPPI = essentialORFsPPI[tmp]

    

  if(debug >= 5){GooddEssentialORFsPPI = GooddEssentialORFsPPI[1:100]  }

  

  sim_names = c( "degreeCutoff","p", "lambda", "meanLS", "medianLS", "R","G", "GompAIC", "WeibAIC")

  sim       = t( c(NA,     NA,   NA,       NA,       NA,        NA,  NA,   NA,      NA))

  sim = data.frame(sim)

  names(sim) = sim_names

  full_age_dir = paste(path, '/', 'popages', sep='')

  system(paste('mkdir ', full_age_dir ))

  

      popAges = numeric(popSize)

      time1 = date()

      j=1; count = 0; 

      while ((j <= popSize) && ( count < popSize*30)) {

        count = count + 1;      

        print(paste("count=",count))

        currentNetworkAge = single_network_failure(lambda, p, pairs, GooddEssentialORFsPPI)

        if (currentNetworkAge > 0) {

          popAges[j] = currentNetworkAge      

          j = j+1

        } 

      }# end of j while-loop, population loop

            

      

      timestamp = format(Sys.time(), "%Y%b%d_%H%M%S")

      age.file.name=paste("cutoff", degreeCutoff, "p", p, "lambda", lambda, 'popsize',popSize, "time",timestamp, "txt", sep="." )

      full_age_file = paste( full_age_dir,'/', age.file.name, sep='')

      

      write.csv( popAges, full_age_file, row.names=F)

        

}