Thursday, December 19, 2013
Tuesday, December 17, 2013
power law network note,
In special case that $t_k = t_0 *(k-1)$, we have $ G = 1/t_0 $. This case also means $\lambda = \lambda_0 / (k-1)$.
My current solutions for the network $R$ and $G$:
From binomial to exponential, Taylor expansion
The approximation of binomial term using an exponential term is a critical step in the development of reliability model of aging. This approximation is perfect in linear term, then have errors in 2nd, 3rd order terms.
So, G= (n-1)/z, where $z$ is $t_0$.
PPI aging simulation, mactower
I logmein into mactower, run the simulation, and then use Dropbox to send back the results.
20131216.ppi.sim.aging_v2_helen.R #for cutoff=4 for degree
20131216.ppi.sim.aging_v2_helen_cutoff5.R
# todo,
# permutation effect on aging?
# lambda ~ 1/connectivity of nodes
rm(list=ls())
require('flexsurv')
#source("/Users/hongqin/lib/R/lifespan.r")
source("lifespan.r")
#setwd("~/projects/0.network.aging.prj/0.ppi.reliability.simulation")
list.files(path='data', )
debug = 0;
#yeast PPI
pairs = read.csv('data/pairs.csv', colClasses=c('character','character'))
#this yeast ppi dataset is consistent with Taiwan group's report.
#essential gene info
essenTb = read.csv('data/SummaryRegressionHetHom2013Oct29.csv', colClasses=rep('character', 9))
#######################
# How do the two data set overlap? DIP seems to contain some questionable orfs
uniq.orf.from.pairs = unique(c(pairs$ORF1, pairs$ORF2))
matches = uniq.orf.from.pairs %in% unique(essenTb$orf)
table(matches)
#FALSE TRUE
# 261 4217
# YDL026W dubious, YAR009C transposable element
dubiousORF = uniq.orf.from.pairs[! matches]
matches = uniq.orf.from.pairs %in% unique(essenTb$orf[essenTb$essenflag=='essential'])
table(matches)
#FALSE TRUE
#3506 972 #this is amazingly consistent with the Taiwan group's report.
matches = uniq.orf.from.pairs %in% unique(essenTb$orf[essenTb$essenflag=='nonessential'])
table(matches)
# FALSE TRUE
# 1287 3191 #this is amazingly consistent with Taiwan group's report.
#remove dubious orfs from PPI
pairs$Removeflag = ifelse( pairs$ORF1 %in%dubiousORF | pairs$ORF2 %in%dubiousORF, T,F )
table(pairs$Removeflag)
#FALSE TRUE
#12180 1487
pairs = pairs[! pairs$Removeflag, ]
table(pairs$Removeflag)
pairs = pairs[,1:2] ##This set of pairs is read for analysis
###############################
# label essential nodes, remove nonesse-nonessen pairs
essentialORFs = essenTb$orf[essenTb$essenflag=='essential']
pairs$essen1 = pairs$ORF1 %in% essentialORFs
pairs$essen2 = pairs$ORF2 %in% essentialORFs
head(pairs)
#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,1:4 ] #only 6279 intxn left
#remove self-intxn, just to make sure
pairs = pairs[ pairs$ORF1 != pairs$ORF2, ]
#how many essen <--> essen intxn?
pairs$inxnEE = pairs$essen1 & pairs$essen2
table(pairs$inxnEE)
# FALSE TRUE
# 4521 1758
#How many essen genes?
tmp = essentialORFs %in% unique(c(pairs$ORF1, pairs$ORF2))
table(tmp)
#FALSE TRUE
#194 958 #So, intxn among essential genes appear to be suppressed.
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; #######!!!!!!
tmp = essentialORFsPPI %in% degreeTb$ORF[degreeTb$Freq>degreeCutoff]
GooddEssentialORFsPPI = essentialORFsPPI[tmp]
###########################
if(debug >= 5){GooddEssentialORFsPPI = GooddEssentialORFsPPI[1:100] }
###########################
# simulate aging
# -> exponential age to all pairs
# -> maximal age for each essential gene
# -> minimal age for all essential modules
set.seed(2013)
lambda_v = 1 / 3^seq(2,5)
#lambda_v = 1 / s
#p_v = seq(0.1, 1.0, by=0.1) ; #the chance that each gene interaction is active at t=0
p_v = c(0.7, 0.7, 0.7, 0.8, 0.8, 0.8, 0.9, 0.9, 0.9, 1.0, 1.0, 1.0) ; #the chance that each gene interaction is active at t=0
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
################################## network simulations
# for debug, lambda = 1/10
# lambda = 1/35; p=0.9
for(lambda in lambda_v) {
# p=0.9, #for debug
for( p in p_v) {
popSize = 100 #too small pop size and too small p can lead to very few living individuals
popAges = numeric(popSize)
time1 = date()
for( j in 1:popSize) {
runningORFs = GooddEssentialORFsPPI
#I could add stochasticity into pairs here.
pairs$active = runif(length(pairs[,1]))
# tmp = pairs$active > 1-p
# table(tmp) / length(tmp) ; #double-check, very good.
ModuleTb = data.frame(runningORFs)
pairs$age = rexp( length(pairs[,1]), rate=lambda ) #exponential ages for pairs
pairs$age = ifelse(pairs$active > (1-p), pairs$age, NA )
#loop every essential genes to identify the module age
for (i in 1:length(runningORFs)) {
myORF = runningORFs[i]
pos1 = grep(myORF, pairs$ORF1)
pos2 = grep(myORF, pairs$ORF2)
if( length( c(pos1,pos2))>=1 ) {
ModuleTb$age.m[i] = max( pairs$age[c(pos1,pos2)], na.rm=T ) #maximal intxn age -> module age
} else {
ModuleTb$age.m[i] = NA;
}
}
#head(ModuleTb);
summary(ModuleTb)
ModuleTb$age.m[ ModuleTb$age.m== -Inf] = 0; #dead births occur when links are not active
currentNetworkAge = min(ModuleTb$age.m)
popAges[j] = currentNetworkAge
}# end of j loop, population loop
time2 = date()
hist(popAges)
summary(popAges)
popAges = popAges[popAges>0]; #remove dead-births, which can occur when p is low
#time1; time2;
#s.tb = calculate.s ( popAges )
#plot( s.tb$s ~ s.tb$t )
#plot( s.tb$s ~ s.tb$t, type='l', log='x' )
lifespanGomp = flexsurvreg(formula = Surv(popAges) ~ 1, dist = 'gompertz') ### Use the flexsurvreg package to fit lifespan data to gompertz or weibull distribution
lifespanWeib = flexsurvreg(formula = Surv(popAges) ~ 1, dist = 'weibull')
c(lifespanWeib$AIC, lifespanGomp$AIC, lifespanWeib$AIC - lifespanGomp$AIC )
sOject = Surv(popAges)
#sim_names = c( "cutoff","p", "lambda", "meanLS", "medianLS", "R","G", "GompAIC", "WeibAIC")
sim = rbind(sim, c( degreeCutoff, p, lambda, mean(popAges), median(popAges),
lifespanGomp$res[1,1], lifespanGomp$res[2,1], lifespanGomp$AIC, lifespanWeib$AIC))
} # end of p-loop
write.csv(sim, file="scePPIaging_sim_20131217_v2_tmp.csv", row.names=F)
} #end of lambda loop
#conditions$avgLS[r] = mean(lifespansTemp)
#conditions$medianLS[r] = median(lifespansTemp)
#conditions$gompShape[r] = lifespanGomp$res[1,1]
#conditions$gompRate[r] = lifespanGomp$res[2,1]
#conditions$gompLogLik[r] = lifespanGomp$loglik
#conditions$gompAIC[r] = lifespanGomp$AIC
#conditions$weibShape[r] = lifespanWeib$res[1,1]
#conditions$weibScale[r] = lifespanWeib$res[2,1]
#conditions$weibLogLik[r] = lifespanWeib$loglik
#conditions$weibAIC[r] = lifespanWeib$AIC
write.csv(sim, file="scePPIaging_sim_20131217_v2_end.csv", row.names=F)
#quite('yes')
20131216.ppi.sim.aging_v2_helen.R #for cutoff=4 for degree
20131216.ppi.sim.aging_v2_helen_cutoff5.R
# todo,
# permutation effect on aging?
# lambda ~ 1/connectivity of nodes
rm(list=ls())
require('flexsurv')
#source("/Users/hongqin/lib/R/lifespan.r")
source("lifespan.r")
#setwd("~/projects/0.network.aging.prj/0.ppi.reliability.simulation")
list.files(path='data', )
debug = 0;
#yeast PPI
pairs = read.csv('data/pairs.csv', colClasses=c('character','character'))
#this yeast ppi dataset is consistent with Taiwan group's report.
#essential gene info
essenTb = read.csv('data/SummaryRegressionHetHom2013Oct29.csv', colClasses=rep('character', 9))
#######################
# How do the two data set overlap? DIP seems to contain some questionable orfs
uniq.orf.from.pairs = unique(c(pairs$ORF1, pairs$ORF2))
matches = uniq.orf.from.pairs %in% unique(essenTb$orf)
table(matches)
#FALSE TRUE
# 261 4217
# YDL026W dubious, YAR009C transposable element
dubiousORF = uniq.orf.from.pairs[! matches]
matches = uniq.orf.from.pairs %in% unique(essenTb$orf[essenTb$essenflag=='essential'])
table(matches)
#FALSE TRUE
#3506 972 #this is amazingly consistent with the Taiwan group's report.
matches = uniq.orf.from.pairs %in% unique(essenTb$orf[essenTb$essenflag=='nonessential'])
table(matches)
# FALSE TRUE
# 1287 3191 #this is amazingly consistent with Taiwan group's report.
#remove dubious orfs from PPI
pairs$Removeflag = ifelse( pairs$ORF1 %in%dubiousORF | pairs$ORF2 %in%dubiousORF, T,F )
table(pairs$Removeflag)
#FALSE TRUE
#12180 1487
pairs = pairs[! pairs$Removeflag, ]
table(pairs$Removeflag)
pairs = pairs[,1:2] ##This set of pairs is read for analysis
###############################
# label essential nodes, remove nonesse-nonessen pairs
essentialORFs = essenTb$orf[essenTb$essenflag=='essential']
pairs$essen1 = pairs$ORF1 %in% essentialORFs
pairs$essen2 = pairs$ORF2 %in% essentialORFs
head(pairs)
#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,1:4 ] #only 6279 intxn left
#remove self-intxn, just to make sure
pairs = pairs[ pairs$ORF1 != pairs$ORF2, ]
#how many essen <--> essen intxn?
pairs$inxnEE = pairs$essen1 & pairs$essen2
table(pairs$inxnEE)
# FALSE TRUE
# 4521 1758
#How many essen genes?
tmp = essentialORFs %in% unique(c(pairs$ORF1, pairs$ORF2))
table(tmp)
#FALSE TRUE
#194 958 #So, intxn among essential genes appear to be suppressed.
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; #######!!!!!!
tmp = essentialORFsPPI %in% degreeTb$ORF[degreeTb$Freq>degreeCutoff]
GooddEssentialORFsPPI = essentialORFsPPI[tmp]
###########################
if(debug >= 5){GooddEssentialORFsPPI = GooddEssentialORFsPPI[1:100] }
###########################
# simulate aging
# -> exponential age to all pairs
# -> maximal age for each essential gene
# -> minimal age for all essential modules
set.seed(2013)
lambda_v = 1 / 3^seq(2,5)
#lambda_v = 1 / s
#p_v = seq(0.1, 1.0, by=0.1) ; #the chance that each gene interaction is active at t=0
p_v = c(0.7, 0.7, 0.7, 0.8, 0.8, 0.8, 0.9, 0.9, 0.9, 1.0, 1.0, 1.0) ; #the chance that each gene interaction is active at t=0
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
################################## network simulations
# for debug, lambda = 1/10
# lambda = 1/35; p=0.9
for(lambda in lambda_v) {
# p=0.9, #for debug
for( p in p_v) {
popSize = 100 #too small pop size and too small p can lead to very few living individuals
popAges = numeric(popSize)
time1 = date()
for( j in 1:popSize) {
runningORFs = GooddEssentialORFsPPI
#I could add stochasticity into pairs here.
pairs$active = runif(length(pairs[,1]))
# tmp = pairs$active > 1-p
# table(tmp) / length(tmp) ; #double-check, very good.
ModuleTb = data.frame(runningORFs)
pairs$age = rexp( length(pairs[,1]), rate=lambda ) #exponential ages for pairs
pairs$age = ifelse(pairs$active > (1-p), pairs$age, NA )
#loop every essential genes to identify the module age
for (i in 1:length(runningORFs)) {
myORF = runningORFs[i]
pos1 = grep(myORF, pairs$ORF1)
pos2 = grep(myORF, pairs$ORF2)
if( length( c(pos1,pos2))>=1 ) {
ModuleTb$age.m[i] = max( pairs$age[c(pos1,pos2)], na.rm=T ) #maximal intxn age -> module age
} else {
ModuleTb$age.m[i] = NA;
}
}
#head(ModuleTb);
summary(ModuleTb)
ModuleTb$age.m[ ModuleTb$age.m== -Inf] = 0; #dead births occur when links are not active
currentNetworkAge = min(ModuleTb$age.m)
popAges[j] = currentNetworkAge
}# end of j loop, population loop
time2 = date()
hist(popAges)
summary(popAges)
popAges = popAges[popAges>0]; #remove dead-births, which can occur when p is low
#time1; time2;
#s.tb = calculate.s ( popAges )
#plot( s.tb$s ~ s.tb$t )
#plot( s.tb$s ~ s.tb$t, type='l', log='x' )
lifespanGomp = flexsurvreg(formula = Surv(popAges) ~ 1, dist = 'gompertz') ### Use the flexsurvreg package to fit lifespan data to gompertz or weibull distribution
lifespanWeib = flexsurvreg(formula = Surv(popAges) ~ 1, dist = 'weibull')
c(lifespanWeib$AIC, lifespanGomp$AIC, lifespanWeib$AIC - lifespanGomp$AIC )
sOject = Surv(popAges)
#sim_names = c( "cutoff","p", "lambda", "meanLS", "medianLS", "R","G", "GompAIC", "WeibAIC")
sim = rbind(sim, c( degreeCutoff, p, lambda, mean(popAges), median(popAges),
lifespanGomp$res[1,1], lifespanGomp$res[2,1], lifespanGomp$AIC, lifespanWeib$AIC))
} # end of p-loop
write.csv(sim, file="scePPIaging_sim_20131217_v2_tmp.csv", row.names=F)
} #end of lambda loop
#conditions$avgLS[r] = mean(lifespansTemp)
#conditions$medianLS[r] = median(lifespansTemp)
#conditions$gompShape[r] = lifespanGomp$res[1,1]
#conditions$gompRate[r] = lifespanGomp$res[2,1]
#conditions$gompLogLik[r] = lifespanGomp$loglik
#conditions$gompAIC[r] = lifespanGomp$AIC
#conditions$weibShape[r] = lifespanWeib$res[1,1]
#conditions$weibScale[r] = lifespanWeib$res[2,1]
#conditions$weibLogLik[r] = lifespanWeib$loglik
#conditions$weibAIC[r] = lifespanWeib$AIC
write.csv(sim, file="scePPIaging_sim_20131217_v2_end.csv", row.names=F)
#quite('yes')
*** Binomial mortality model, survival function for network aging modeling !!!!
INVERSE Function of binomial CDF
Monday, December 16, 2013
sharpshoot and its symbotic bacteria, good model for bacteria evolution study (condon bias), small Ne of symbtiotic bacteria.
Nancy moran: sharpshoot and its symbotic bacteria, good model for bacteria evolution study (condon bias), small Ne of symbtiotic bacteria.
Taylor expansion of Gompertz function
Derivatives of u(t) near t=0 can be found by chain rules:
http://en.wikipedia.org/wiki/Taylor_series
Results repeated in Mathematica:
Friday, December 13, 2013
Yeast cell cycle and morphology (S. cerevisiae). Leland Hartwell, 1974
This diagram does not describe the relative time-duration of each phase. The relative durations will be influenced by growth state. In log-phase, based on the DNA content fluorescent signal, it can be estimated that 1N and 2N are the major peaks. So, G2-M phase should be relatively quick.
Hartwell 1974
Wednesday, December 11, 2013
Functions for life span distribution. GG91
page 43, two parameter Gompertz model
page 44, Gompertz-Makeham model
Quadratic form of Gompertz-Makeham model, El Shaarawi et al 1974
u(t) = M + R exp( at + bt^2)
Logistic equations like the Perks equation, which can be derived from GG91's avalanche-like aging model.
Weibull model (1951): u(t) = B * t^c
page 46, GG91's generalized bionomial mortality model:
u(t) = A + (b + c t)^n
GG91 argues that when b >>c, bionomial model lead to R=b^n, G=nc/b
When b<<c, it is the Weibull model.
page 44, Gompertz-Makeham model
Quadratic form of Gompertz-Makeham model, El Shaarawi et al 1974
u(t) = M + R exp( at + bt^2)
Logistic equations like the Perks equation, which can be derived from GG91's avalanche-like aging model.
Weibull model (1951): u(t) = B * t^c
page 46, GG91's generalized bionomial mortality model:
u(t) = A + (b + c t)^n
GG91 argues that when b >>c, bionomial model lead to R=b^n, G=nc/b
When b<<c, it is the Weibull model.
Estimation of mortality rate at time $t$, GG91
page 41, Sacher 1956, 1966
\mu(t) = 1/(2 dt) log( \frac{s(t-dt}{t+dt} )
Other methods discussed are: Cutler and Ederer 1958
GG91 argued that Saccher56 formula is better than CutlerEdere58 method.
\mu(t) = 1/(2 dt) log( \frac{s(t-dt}{t+dt} )
Other methods discussed are: Cutler and Ederer 1958
GG91 argued that Saccher56 formula is better than CutlerEdere58 method.
Variability in lifespan, GG91 notes
page 31,
There are two extreme positions to explain lifetime variability: heterogeneity hypothesis and stochastic hypothesis.
page 37
GG91 argued that stochastic variation should be considered in addition to environment and genetic factors. GG91 cited Sacher 1977 who tried to quantify stochastic variation from lifespan distribution.
GG91 viewed stochastic variation as 'kinetic' variation.
There are two extreme positions to explain lifetime variability: heterogeneity hypothesis and stochastic hypothesis.
page 37
GG91 argued that stochastic variation should be considered in addition to environment and genetic factors. GG91 cited Sacher 1977 who tried to quantify stochastic variation from lifespan distribution.
GG91 viewed stochastic variation as 'kinetic' variation.
Heritability of lifespan, GG91, notes
GG91, p33, cite Jacquard, 1982. Heritability coefficient for human lifespan is 0.16. GG91 interpret this as: Even if the lifespan of both parents exceed the mean lifespan by 20 years, their offspring will gain on average only 0.16x20=3.2 extra years from their genetic 'inheritance'.
Page 34. GG91 then argued that "selection for an increase in life span is ineffective". GG91 then used the heterozygote Aa example to make their point.
page 35,
GG91 cited that heteozygotes often live longer than homozygotes. Bileva et al 1978, 1981, Nekrasov 1981, Shakhbazov 1980.
Page 34. GG91 then argued that "selection for an increase in life span is ineffective". GG91 then used the heterozygote Aa example to make their point.
page 35,
GG91 cited that heteozygotes often live longer than homozygotes. Bileva et al 1978, 1981, Nekrasov 1981, Shakhbazov 1980.
History for the biology of life span, GG91, reading notes
John Graunt, 1662, first life table for London residents
Leonhard Euler, 1760, general investigations into mortality and the multiplication of the human race
Benjanmin Gompertz, 1825, Gompertz model. (Was Gompertz the first person to use the term "the force of mortality"?)
William Makeham, 1860, Gompertz-Makeham model
Rosset 1979 reviewed the history of biology of lifespan in a Polish article.
Leonhard Euler, 1760, general investigations into mortality and the multiplication of the human race
Benjanmin Gompertz, 1825, Gompertz model. (Was Gompertz the first person to use the term "the force of mortality"?)
William Makeham, 1860, Gompertz-Makeham model
Rosset 1979 reviewed the history of biology of lifespan in a Polish article.
Comp biol projects for undergraduates
SGU CLS-LOH data analysis
Yeast NGS analysis
Network reliability simulation
LOH peak simulation
-------------
Toggle switch simulation
SIR simulation
Yeast NGS analysis
Network reliability simulation
LOH peak simulation
-------------
Toggle switch simulation
SIR simulation
Tuesday, December 10, 2013
Useful iPad apps for education
https://itunes.apple.com/us/app/explain-everything/id431493086
http://www.educatorstechnology.com/2013/12/20-ipad-apps-to-showcase-students.html?utm_source=dlvr.it&utm_medium=twitter
Subscribe to:
Posts (Atom)