data<-read.csv("DNAdamage_Data.csv") head(data) ### Generate inferred melanin values based upon photo grey values melPhoto.df<-read.csv("MelaninPhotoRegression_DNAdamageExpt.csv") model1<-lm(melanin.per.length~photo.grey,data=melPhoto.df) inferred.melanin<-(predict(model1,newdata=data)) data<-cbind(data,inferred.melanin) # subset data sierra<-droplevels(subset(data,data$region=="sierra")) oly<-droplevels(subset(data,data$region=="oly")) ############# # TEST STATISTICS FOR RESULTS NARRATIVE ############# # Compare melanin between regions (Sierra/Oly) # First TEST FOR NORMALITY: shapiro.test(data$inferred.melanin) # NOT normal # Mann-Whitney U test wilcox.test(inferred.melanin~region,data=data) # W = 372, p < 2.2e-16 # WITHIN each region, do clones differ in melanin? summary(lm(inferred.melanin~clone,data=sierra)) #F3,83 = 4.382, p = 0.0065 summary(lm(inferred.melanin~clone,data=oly)) #F3,74 = 2.985, p = 0.037 #### ELISA RESULTS ## CHECK FOR NORMALITY IN ELISA absorbance FIRST shapiro.test(data$ELISA) # They are normal: # p = 0.34, W = 0.99 ######################## ### SIERRA MODELS ######################## Smelaninmodel<-lm(ELISA~inferred.melanin,data=sierra) Sclonemodel<-lm(ELISA~clone,data=sierra) Smodel2<-lm(ELISA~inferred.melanin+clone,data=sierra) Smodel3<-lm(ELISA~clone*inferred.melanin,data=sierra) summary(Smelaninmodel) # melanin significant predictor of DNA damage # F1,85 = 14.57 p=0.00026, R2 = 0.15 anova(Sclonemodel,Smodel2) # melanin still significant after accounting for clone # p= 6.86e-5, F1,82 = 17.599 # confirm correspondence between F val and p val: 1-pf(17.599,1,82) anova(Smodel2,Smodel3) # melanin-clone interaction not significant for Sierra # F3,82=1.5028, p=0.22 # confirm correspondence between F val and p val: 1-pf(1.5028,3,82) # So best Sierra model is Smodel2 (it includes clone) # But for the graph best-fit line that doesn't work, so we use Smelaninmodel coef(Smelaninmodel) ######################## ### END SIERRA MODELS ######################## ######################## ### OLYMPIC MODELS ######################## ### Omelaninmodel<-lm(ELISA~inferred.melanin,data=oly) Oclonemodel<-lm(ELISA~clone,data=oly) Omodel2<-lm(ELISA~clone+inferred.melanin,data=oly) Omodel3<-lm(ELISA~clone*inferred.melanin,data=oly) summary(Omelaninmodel) #F1,76 = 1.114, p=0.29, R2=0.014 anova(Oclonemodel,Omodel2) # melanin stil NOT significant after accounting for clone # p= 0.74, F1,73 = 0.1138 # confirm correspondence between F val and p val: 1-pf(0.1138,1,74) ######################## ### END OLYMPIC MODELS ######################## ####### CVs for DNA damage CVdamage.oly<-mean(sd(oly$ELISA)/mean(oly$ELISA)) CVdamage.sierra<-mean(sd(sierra$ELISA)/mean(sierra$ELISA)) CVdamage.oly; CVdamage.sierra # 0.36; 0.43 # So based upon these coefficients of variaiton, Sierra is more variable than olympic #### To plot each Sierra clone fit separately: Sclone1<-droplevels(subset(sierra,sierra$clone=="C1")) Sclone2<-droplevels(subset(sierra,sierra$clone=="P2")) Sclone3<-droplevels(subset(sierra,sierra$clone=="U3")) Sclone4<-droplevels(subset(sierra,sierra$clone=="W3")) ##### # Clone-specific models nullclone1<-lm(ELISA~1,data=Sclone1) clone1<-lm(ELISA~inferred.melanin,data=Sclone1) anova(nullclone1,clone1) # p=0.11 nullclone2<-lm(ELISA~1,data=Sclone2) clone2<-lm(ELISA~inferred.melanin,data=Sclone2) anova(nullclone2,clone2) # p=0.09 nullclone3<-lm(ELISA~1,data=Sclone3) clone3<-lm(ELISA~inferred.melanin,data=Sclone3) anova(nullclone3,clone3) # p=0.005 nullclone4<-lm(ELISA~1,data=Sclone4) clone4<-lm(ELISA~inferred.melanin,data=Sclone4) anova(nullclone4,clone4) # p=0.016 ### ###################### ### FIGURE 2 ###################### col1=rgb(213,94,0,maxColorValue=255) #Vermillion col2=rgb(0,114,178,maxColorValue=255) #Blue col3=rgb(0,158,155,maxColorValue=255) #Bluish green col4=rgb(204,121,167,maxColorValue=255) #Reddish purple col5=rgb(86,180,233,maxColorValue=255) #Sky Blue quartz(width=7.12,height=6) # Use the full range of data to figure out xlim: summary(data$inferred.melanin) xmin=0.2; xmax=0.55 par(las=1,bty="l") par(mgp=c(3.5,1,0),cex.axis=1.5,cex.lab=2,las=1,mar=c(4.5,5.5,1,1)) plot(ELISA~inferred.melanin,data=oly, xlim=c(xmin,xmax),ylim=c(0,1), pch=18, xlab=expression(paste("Inferred melanin (",mu,"g/mm)")), ylab="DNA damage (ELISA abs.)", cex=1.5,col=col4,lwd=2, bty="l") # OVERALL SIERRA LINE curve( coef(Smelaninmodel)[1]+coef(Smelaninmodel)[2]*x, from=0.25,to=0.55, add=T, col="black",lwd=5 ) # # Oly points again # points(ELISA~inferred.melanin,data=oly, # pch=18,cex=1.2,col=col4,lwd=2) # Sierra lines (4 gts) points(ELISA~inferred.melanin,data=Sclone1, pch=1,cex=1.5,col=col1,lwd=2) points(ELISA~inferred.melanin,data=Sclone2, pch=2,cex=1.5,col=col2,lwd=2) points(ELISA~inferred.melanin,data=Sclone3, pch=3,cex=1.5,col=col3,lwd=2) points(ELISA~inferred.melanin,data=Sclone4, pch=0,cex=1.5,col=col5,lwd=2) # add fitted lines for sierra clones, separate line for each clone # Plot clone-specific LMs curve( coef(clone1)[1]+coef(clone1)[2]*x, from=0.25,to=0.55, add=T, col=col1,lty=2,lwd=3 ) curve( coef(clone2)[1]+coef(clone2)[2]*x, from=0.25,to=0.55, add=T, col=col2,lty=3,lwd=3 ) curve( coef(clone3)[1]+coef(clone3)[2]*x, from=0.25,to=0.55, add=T, col=col3,lty=4,lwd=3 ) curve( coef(clone4)[1]+coef(clone4)[2]*x, from=0.25,to=0.55, add=T, col=col5,lty=5,lwd=3 ) # OVERALL SIERRA LINE again, so it is plotted on top curve( coef(Smelaninmodel)[1]+coef(Smelaninmodel)[2]*x, from=0.25,to=0.55, add=T, col="black",lwd=5 ) legend("topright",legend=c("Sierra, gt=C1","Sierra, gt=P2","Sierra, gt=U3","Sierra, gt=W3","Olympic (all gts)","Sierra (fit for all gts)"), #title=expression(underline(Region)), bty="n",cex=1.0,pt.cex=1.5, pch=c(1,2,3,0,18,NA), lty=c(2,3,4,5,0,1), col=c(col1,col2,col3,col5,col4,"black"), pt.lwd=2,lwd=c(3,3,3,3,3,5), seg.len=6) # Slopes for each Sierra clone's line coef(clone1);coef(clone2);coef(clone3);coef(clone4);