library(arules) #for discretize function library(olsrr) #for ols_test_breusch_pagan() ### READING DATA ### dat <- read.csv(file="Luecke_DrosophilaLegSize_Discs_LarvaData.csv", header=T) # third instar disc data formatting and partitioning # parameters combining timing methods, always in order kb866F, kb866M, prolF, prolM. Used in make.discs.df() full_larva_days <- c(6.7, 7.0, 7.2, 8.1) # time to pupa formation percent_L3 <- c(0.5, 0.5, 0.5, 0.5) # assume L3 at 50% of larva development discs0 <- make.discs.df(dat, full_larva_days, percent_L3) discs1 <- discs0[complete.cases(discs0$T1DAPI),] discs <- discs1[complete.cases(discs1$T2DAPI),] discs.ph3 <- complete.cols(discs) discs.dcp1 <- discs[complete.cases(discs),] discs.time0 <- discs[complete.cases(discs$relativeTP),] discs.time <- discs.time0[order(discs.time0$relativeTP),] discs.time.ph3 <- complete.cols(discs.time) discs.wpp <- discs[discs$timepoint=='wpp',] # embryo disc precursor data formatting embryo <- read.csv(file = "Luecke_DrosophilaLegSize_Discs_EmbryoData.csv") embryo$ratio <- embryo$T1/embryo$T2 embryo$diff <- embryo$T1 - embryo$T2 embryo$ID <- paste(embryo$species, embryo$sex) ### FIGURES AND STATISTICAL TESTS ### # Fig 3A: cell count difference in embryonic precursors boxplot(embryo$diff ~ embryo$ID, ylab="T1 Excess Precursor Cells") abline(h=0, col='red', lty=3) # associated statistical tests embryo.diff <- aov(diff ~ sex*species, data=embryo) summary(embryo.diff) coefficients(embryo.diff) # Fig 3A': cell count difference at wpp plot(factor(discs.wpp$ID), discs.wpp$DAPI.diff, ylab="T1 Total Cell Excess at Pupariation") abline(h=0, col='red', lty=3) # associated statistical tests by(discs.wpp, discs.wpp$ID, function(x) c(nrow(x), mean(x$DAPI.ratio), mean(x$DAPI.diff))) by(discs.wpp, discs.wpp$ID, function(x) t.test(x$DAPI.diff)) by(discs.wpp, discs.wpp$species, function(x) t.test(x$DAPI.diff ~ x$sex)) anova.wppdiff <- aov(DAPI.diff ~ species*sex, data = discs.wpp) summary(anova.wppdiff) coefficients(anova.wppdiff) #downsample prolongata samples to see if significant dimorphism at smallest sample size (9 for kb866 f) smallest_n <- min(by(discs.wpp, discs.wpp$ID, function(x) c(nrow(x)))) ps.wd <- vector() dfs.wd <- vector() ts.wd <- vector() for(i in 1:1000){ ps.wd[i] <- t.test(sample(discs.wpp[discs.wpp$ID=="prol m",]$DAPI.diff, smallest_n), sample(discs.wpp[discs.wpp$ID=="prol f",]$DAPI.diff, smallest_n))$p.value ts.wd[i] <- t.test(sample(discs.wpp[discs.wpp$ID=="prol m",]$DAPI.diff, smallest_n), sample(discs.wpp[discs.wpp$ID=="prol f",]$DAPI.diff, smallest_n))$statistic dfs.wd[i] <- t.test(sample(discs.wpp[discs.wpp$ID=="prol m",]$DAPI.diff, smallest_n), sample(discs.wpp[discs.wpp$ID=="prol f",]$DAPI.diff, smallest_n))$parameter } mean(ps.wd) mean(ts.wd) mean(dfs.wd) # multiple test correction p.adjust(c(by(discs.wpp, discs.wpp$species, function(x) t.test(x$DAPI.diff ~ x$sex)$p.value)[1], mean(ps.wd)), method = "holm") # test for embiggened T1 in males and females of both samples, again downsampled p.holder <- vector() wpp.diff.ps <- by(discs.wpp, discs.wpp$ID, function(x) { for(i in 1:1000){p.holder[i]<-t.test(sample(x$DAPI.diff, smallest_n))$p.value} return(mean(p.holder)) }) t.holder <- vector() wpp.diff.ts <- by(discs.wpp, discs.wpp$ID, function(x) { for(i in 1:1000){t.holder[i]<-t.test(sample(x$DAPI.diff, smallest_n))$statistic} return(mean(t.holder)) }) df.holder <- vector() wpp.diff.dfs <- by(discs.wpp, discs.wpp$ID, function(x) { for(i in 1:1000){df.holder[i]<-t.test(sample(x$DAPI.diff, smallest_n))$parameter} return(mean(df.holder)) }) wpp.diff.ts wpp.diff.dfs p.adjust(wpp.diff.ps, method="holm") # Fig 3B: cell difference over time, uses functions from end of script plot.a.byreltime.splitbyID(discs.time, "DAPI.diff", "Difference (T1 - T2)") # associated statistical tests discs.time.kbF <- discs.time[discs.time$ID=="kb866 f",] discs.time.kbM <- discs.time[discs.time$ID=="kb866 m",] discs.time.prF <- discs.time[discs.time$ID=="prol f",] discs.time.prM <- discs.time[discs.time$ID=="prol m",] reg.diffreltime.kbF <- lm(DAPI.diff ~ relativeTP, discs.time.kbF) reg.diffreltime.kbM <- lm(DAPI.diff ~ relativeTP, discs.time.kbM) reg.diffreltime.prF <- lm(DAPI.diff ~ relativeTP, discs.time.prF) reg.diffreltime.prM <- lm(DAPI.diff ~ relativeTP, discs.time.prM) summary(reg.diffreltime.kbF) summary(reg.diffreltime.kbM) summary(reg.diffreltime.prF) summary(reg.diffreltime.prM) ols_test_breusch_pagan(reg.diffreltime.kbF) ols_test_breusch_pagan(reg.diffreltime.kbM) ols_test_breusch_pagan(reg.diffreltime.prF) ols_test_breusch_pagan(reg.diffreltime.prM) # Fig 3C, difference in mitotic cells and apoptotic cells par(mfrow=c(1,2), mar=c(4.1,4.1,2.5,1.1)) plot(factor(discs.ph3$ID), discs.ph3$PH3.diff, ylab="T1 excess of mitotic cells") abline(h=0, col='red', lty=3) plot(factor(discs.dcp1$ID), discs.dcp1$DCP1.diff, ylab="T1 excess of apoptotic cells") abline(h=0, col='red', lty=3) # associatd statistical tests anova.ph3diff <- aov(PH3.diff ~ species*sex, data=discs.ph3) summary(anova.ph3diff) coefficients(anova.ph3diff) # sex dimorphism by(discs.ph3, discs.ph3$species, function(x) t.test(x$PH3.diff~x$sex)) p.adjust(by(discs.ph3, discs.ph3$species, function(x) t.test(x$PH3.diff~x$sex)$p.value), method='holm') # T1 bias by(discs.ph3, discs.ph3$ID, function(x) t.test(x$PH3.diff)) by(discs.ph3, discs.ph3$ID, function(x) nrow(x)) # bootstrap downsampling the male prol to female sample size (86) and see the average p value, still significantly different from 0: ps.m <- vector() ts.m <- vector() dfs.m <- vector() for(i in 1:10000){ ps.m[i]<-t.test(sample(discs.ph3[discs.ph3$ID=="prol m",]$PH3.diff, 86))$p.value ts.m[i]<-t.test(sample(discs.ph3[discs.ph3$ID=="prol m",]$PH3.diff, 86))$statistic dfs.m[i]<-t.test(sample(discs.ph3[discs.ph3$ID=="prol m",]$PH3.diff, 86))$parameter } mean(ps.m) mean(ts.m) mean(dfs.m) # compared to females which are not: ps.f <- vector() for(i in 1:10000){ps.f[i]<-t.test(sample(discs.ph3[discs.ph3$ID=="prol f",]$PH3.diff, 86))$p.value} mean(ps.f) p.adjust(c(mean(ps.f), mean(ps.m)), method="holm") #but downsampling to kb866 male level (34) loses significance ps.m2 <- vector() for(i in 1:10000){ps.m2[i]<-t.test(sample(discs.ph3[discs.ph3$ID=="prol m",]$PH3.diff, 34))$p.value} mean(ps.m2) # statistics for apoptosis: anova.dcp1diff <- aov(DCP1.diff ~ species*sex, data=discs.dcp1) summary(anova.dcp1diff) coefficients(anova.dcp1diff) by(discs.dcp1, discs.dcp1$species, function(x) t.test(x$DCP1.diff~x$sex)) by(discs.dcp1, discs.dcp1$ID, function(x) t.test(x$DCP1.diff)) # Fig 3D timing of mitotic difference and size divergence par(mfrow=c(2,1), oma=c(1,0,0,0)) brks = 8 par(mar=c(1,5,3,2)) boxplot(PH3.diff ~ discretize(relativeTP, method = "interval", breaks = brks), data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], xaxt="n", ylab="T1 Mitotic Cell Excess") abline(h=0, col='red', lty=3) par(mar=c(4,5,0,2)) boxplot(DAPI.diff ~ discretize(relativeTP, method = "interval", breaks = brks), data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], ylab="T1 Total Cell Excess", xlab="Percent Larval Development") abline(h=0, col='red', lty=3) # associated statistics for each discritized bin by(data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], INDICES=discretize(discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',]$relativeTP, method = "interval", breaks = brks), function(x) t.test(x$PH3.diff)) by(data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], INDICES=discretize(discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',]$relativeTP, method = "interval", breaks = brks), function(x) t.test(x$DAPI.diff)) ph3diff.pvalues <- by(data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], INDICES=discretize(discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',]$relativeTP, method = "interval", breaks = brks), function(x) t.test(x$PH3.diff)[3]) celldiff.pvalues <- by(data=discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',], INDICES=discretize(discs.time.ph3[discs.time.ph3$species=='prol',][discs.time.ph3$sex=='m',]$relativeTP, method = "interval", breaks = brks), function(x) t.test(x$DAPI.diff)[3]) p.adjust(ph3diff.pvalues, method="holm") p.adjust(celldiff.pvalues, method="holm") # Fig 4 compenstatory dynamics par(mfrow=c(2,2), mar=c(4.1,4.1,2.5,1.1)) compensation.MIdiff.lms <- plot.compensation.splitbyID.MIdiff(discs.ph3) # associated statistics for(i in 1:4) { print(names(compensation.MIdiff.lms)[i]) print(summary(compensation.MIdiff.lms[[i]])) } # Supplemental Figure 1: divergence of cell number over time par(mfrow=c(2,2), mar=c(4.1,4.1,2.5,1.1)) plot.bothab.byreltime.splitbyID(discs.time, "T1DAPI", "T2DAPI", "Disc Cells") # statistics for regression slopes of above plot discs.time.kbF <- discs.time[discs.time$ID=="kb866 f",] discs.time.kbM <- discs.time[discs.time$ID=="kb866 m",] discs.time.prF <- discs.time[discs.time$ID=="prol f",] discs.time.prM <- discs.time[discs.time$ID=="prol m",] reg.T1reltime.kbF <- lm(T1DAPI ~ relativeTP, discs.time.kbF) reg.T1reltime.kbM <- lm(T1DAPI ~ relativeTP, discs.time.kbM) reg.T1reltime.prF <- lm(T1DAPI ~ relativeTP, discs.time.prF) reg.T1reltime.prM <- lm(T1DAPI ~ relativeTP, discs.time.prM) reg.T2reltime.kbF <- lm(T2DAPI ~ relativeTP, discs.time.kbF) reg.T2reltime.kbM <- lm(T2DAPI ~ relativeTP, discs.time.kbM) reg.T2reltime.prF <- lm(T2DAPI ~ relativeTP, discs.time.prF) reg.T2reltime.prM <- lm(T2DAPI ~ relativeTP, discs.time.prM) summary(reg.T1reltime.kbF) summary(reg.T1reltime.kbM) summary(reg.T1reltime.prF) summary(reg.T1reltime.prM) summary(reg.T2reltime.kbF) summary(reg.T2reltime.kbM) summary(reg.T2reltime.prF) summary(reg.T2reltime.prM) coefficients(reg.T1reltime.kbF) confint(reg.T1reltime.kbF, 'relativeTP', level=0.95) confint(reg.T1reltime.kbF, 'relativeTP', level=0.90) coefficients(reg.T2reltime.kbF) confint(reg.T2reltime.kbF, 'relativeTP', level=0.95) confint(reg.T2reltime.kbF, 'relativeTP', level=0.90) coefficients(reg.T1reltime.kbM) confint(reg.T1reltime.kbM, 'relativeTP', level=0.95) confint(reg.T1reltime.kbM, 'relativeTP', level=0.90) coefficients(reg.T2reltime.kbM) confint(reg.T2reltime.kbM, 'relativeTP', level=0.95) confint(reg.T2reltime.kbM, 'relativeTP', level=0.90) coefficients(reg.T1reltime.prF) confint(reg.T1reltime.prF, 'relativeTP', level=0.95) confint(reg.T1reltime.prF, 'relativeTP', level=0.90) coefficients(reg.T2reltime.prF) confint(reg.T2reltime.prF, 'relativeTP', level=0.95) confint(reg.T2reltime.prF, 'relativeTP', level=0.90) coefficients(reg.T1reltime.prM) confint(reg.T1reltime.prM, 'relativeTP', level=0.95) confint(reg.T1reltime.prM, 'relativeTP', level=0.90) coefficients(reg.T2reltime.prM) confint(reg.T2reltime.prM, 'relativeTP', level=0.95) confint(reg.T2reltime.prM, 'relativeTP', level=0.90) # Supplemental Figure 2: Difference in MI by disc size compensation.MIdiff.T1cellcount.lms <- plot.compensation.splitbyID.MIdiff.T1cellcount(discs.ph3) for(i in 1:4) { print(names(compensation.MIdiff.T1cellcount.lms)[i]) print(summary(compensation.MIdiff.T1cellcount.lms[[i]])) } ### FUNCTIONS ### make.discs.df <- function(dat, days, pL3) { dat$ID <- paste(dat$species, dat$sex, sep=" ") dat$complete <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {1} else {0}) dat$relativeTP <- sapply(1:nrow(dat), function(x) if(dat$method[x]=='AEL') {as.numeric(as.character(dat$timepoint[x])) / days[which(levels(factor(dat$ID))==dat$ID[x])]} else if(dat$method[x]=='AL3') {(as.numeric(as.character(dat$timepoint[x]))/24 + days[which(levels(factor(dat$ID))==dat$ID[x])]*pL3[which(levels(factor(dat$ID))==dat$ID[x])]) / days[which(levels(factor(dat$ID))==dat$ID[x])]} else if(dat$method[x]=='dev') { if(dat$timepoint[x]=='wpp') {1.0} else if(dat$timepoint[x]=='wL3') {0.95} } else {NA} ) dat$T1aCELLSIZE <- dat$T1aDAPI/dat$T1aSIZE dat$T1aNET <- dat$T1aPH3 - dat$T1aDCP1 dat$T1aMI <- dat$T1aPH3/dat$T1aDAPI dat$T1aAI <- dat$T1aDCP1/dat$T1aDAPI dat$T1aRATE <- dat$T1aNET/dat$T1aDAPI dat$T1bCELLSIZE <- dat$T1bDAPI/dat$T1bSIZE dat$T1bNET <- dat$T1bPH3 - dat$T1bDCP1 dat$T1bMI <- dat$T1bPH3/dat$T1bDAPI dat$T1bAI <- dat$T1bDCP1/dat$T1bDAPI dat$T1bRATE <- dat$T1bNET/dat$T1bDAPI dat$T2aCELLSIZE <- dat$T2aDAPI/dat$T2aSIZE dat$T2aNET <- dat$T2aPH3 - dat$T2aDCP1 dat$T2aMI <- dat$T2aPH3/dat$T2aDAPI dat$T2aAI <- dat$T2aDCP1/dat$T2aDAPI dat$T2aRATE <- dat$T2aNET/dat$T2aDAPI dat$T2bCELLSIZE <- dat$T2bDAPI/dat$T2bSIZE dat$T2bNET <- dat$T2bPH3 - dat$T2bDCP1 dat$T2bMI <- dat$T2bPH3/dat$T2bDAPI dat$T2bAI <- dat$T2bDCP1/dat$T2bDAPI dat$T2bRATE <- dat$T2bNET/dat$T2bDAPI dat$T1DAPI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aDAPI[x] + dat$T1bDAPI[x])/2} else {dat$T1aDAPI[x]}) dat$T1SIZE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aSIZE[x] + dat$T1bSIZE[x])/2} else {dat$T1aSIZE[x]}) dat$T1PH3 <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aPH3[x] + dat$T1bPH3[x])/2} else {dat$T1aPH3[x]}) dat$T1DCP1 <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aDCP1[x] + dat$T1bDCP1[x])/2} else {dat$T1aDCP1[x]}) dat$T1CELLSIZE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aCELLSIZE[x] + dat$T1bCELLSIZE[x])/2} else {dat$T1aCELLSIZE[x]}) dat$T1NET <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aNET[x] + dat$T1bNET[x])/2} else {dat$T1aNET[x]}) dat$T1MI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aMI[x] + dat$T1bMI[x])/2} else {dat$T1aMI[x]}) dat$T1AI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aAI[x] + dat$T1bAI[x])/2} else {dat$T1aAI[x]}) dat$T1RATE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T1bDAPI[x])) {(dat$T1aRATE[x] + dat$T1bRATE[x])/2} else {dat$T1aRATE[x]}) dat$T2DAPI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aDAPI[x] + dat$T2bDAPI[x])/2} else {dat$T2aDAPI[x]}) dat$T2SIZE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aSIZE[x] + dat$T2bSIZE[x])/2} else {dat$T2aSIZE[x]}) dat$T2PH3 <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aPH3[x] + dat$T2bPH3[x])/2} else {dat$T2aPH3[x]}) dat$T2DCP1 <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aDCP1[x] + dat$T2bDCP1[x])/2} else {dat$T2aDCP1[x]}) dat$T2CELLSIZE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aCELLSIZE[x] + dat$T2bCELLSIZE[x])/2} else {dat$T2aCELLSIZE[x]}) dat$T2NET <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aNET[x] + dat$T2bNET[x])/2} else {dat$T2aNET[x]}) dat$T2MI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aMI[x] + dat$T2bMI[x])/2} else {dat$T2aMI[x]}) dat$T2AI <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aAI[x] + dat$T2bAI[x])/2} else {dat$T2aAI[x]}) dat$T2RATE <- sapply(1:nrow(dat), function(x) if(!is.na(dat$T2bDAPI[x])) {(dat$T2aRATE[x] + dat$T2bRATE[x])/2} else {dat$T2aRATE[x]}) dat$DAPI.diff <- dat$T1DAPI - dat$T2DAPI dat$DAPI.ratio <- dat$T1DAPI/dat$T2DAPI dat$SIZE.ratio <- dat$T1SIZE/dat$T2SIZE dat$PH3.diff <- dat$T1PH3 - dat$T2PH3 dat$PH3.ratio <- dat$T1PH3/dat$T2PH3 dat$DCP1.diff <- dat$T1DCP1 - dat$T2DCP1 dat$DCP1.ratio <- dat$T1DCP1/dat$T2DCP1 dat$CELLSIZE.ratio <- dat$T1CELLSIZE/dat$T2CELLSIZE dat$MI.ratio <- dat$T1MI/dat$T2MI dat$AI.ratio <- dat$T1AI/dat$T2AI dat$NET.diff <- dat$T1NET - dat$T2NET dat$NET.ratio <- dat$T1NET/dat$T2NET dat$RATE.ratio <- dat$T1RATE/dat$T2RATE discs <- dat[,c('species', 'sex', 'batch', 'set', 'method', 'complete', 'timepoint', 'relativeTP', 'ID', 'T1DAPI', 'T1SIZE', 'T1CELLSIZE', 'T1PH3', 'T1DCP1', 'T1NET', 'T1MI', 'T1AI', 'T1RATE', 'T2DAPI', 'T2SIZE', 'T2CELLSIZE', 'T2PH3', 'T2DCP1', 'T2NET', 'T2MI', 'T2AI', 'T2RATE', 'DAPI.diff', 'DAPI.ratio', 'SIZE.ratio', 'CELLSIZE.ratio', 'PH3.diff', 'PH3.ratio', 'DCP1.diff', 'DCP1.ratio', 'MI.ratio', 'AI.ratio', 'NET.diff', 'NET.ratio', 'RATE.ratio')] return(discs) } complete.cols <- function(dat) { c <- c() for(i in 1:ncol(dat)) { if(sum(complete.cases(dat[,i])) == length(complete.cases(dat[,i]))){c <- c(c,i)} } return(dat[,c]) } plot.a.byreltime <- function(dat,astring,ylab,title,ylim) { a <- which(names(dat) == astring) lma <- lm(dat[,a]~dat$relativeTP) if(summary(lma)$coefficients[2,4] < 0.05) {linetype = 1} else {linetype = 3} plot(dat$relativeTP, dat[,a], main=title, ylab=ylab, xlab="percent larval development", ylim=ylim) abline(lma, lty=linetype) abline(h=0, col='red') } plot.a.byreltime.splitbyID <- function(dat, astring, ylab) { a <- which(names(dat) == astring) by(dat, dat$ID, function(x) plot.a.byreltime(x,astring,ylab,x$ID[1],c(min(dat[,a]),max(dat[,a]))) ) return(NULL) } plot.compensation.splitbyID.MIdiff <- function(dat) { by(dat, dat$ID, function(x) { identify <- x$ID[1] lmx <- lm(x$T1MI-x$T2MI ~ x$DAPI.ratio) if(summary(lmx)$coefficients[2,4] < 0.05) {linetype = 1} else {linetype = 3} plot(x$DAPI.ratio, x$T1MI-x$T2MI, main=identify, xlab="Cell Number: T1/T2", ylab="Mitotic Index: T1 - T2") abline(v=1) abline(h=0) abline(lmx, col='red', lty=linetype) } ) lms <- by(dat, dat$ID, function(x) lm(x$T1MI-x$T2MI ~ x$DAPI.ratio)) return(lms) } plot.compensation.splitbyID.MIdiff.T1cellcount <- function(dat) { by(dat, dat$ID, function(x) { identify <- x$ID[1] lmx <- lm(x$T1MI-x$T2MI ~ x$T1DAPI) if(summary(lmx)$coefficients[2,4] < 0.05) {linetype = 1} else {linetype = 3} plot(x$T1DAPI, x$T1MI-x$T2MI, main=identify, xlab="T1 Cell Number", ylab="Mitotic Index: T1 - T2") abline(h=0) abline(lmx, col='red', lty=linetype) } ) lms <- by(dat, dat$ID, function(x) lm(x$T1MI-x$T2MI ~ x$T1DAPI)) return(lms) }