Pyrosequencing Data

Color Manipulation Experiment

This is the full set of pyrosequencing data for 371 samples submitted to EpigenDX. Note that two samples (18044 and 18045) were confused in the raw data. Since I couldn’t tell which was which and only had data for one I have removed that sample. I also removed samples entirely when greater than 1/2 of the assay was filled with 0s. I believe that those sample have problems with the assay and it isn’t clear if the non-zero numbers are reliably. Similarly, some sample in the dataset do include 0s that may reflect missing data rather than true 0 values. I remove those in the code below, but it is likely worth checking both with and without 0s in some cases (e.g., some particular CpGs have very low average methylation values and 0s might be real at those sites).

Read and organize the data

First I’m just reading in and organizing the data. For now, I only want to work with the 2017 color manipulation experiment. This file also includes the 2018 color-stressor experiment and the 2018 brain collection experiment, so I’ve removed those here.

library(plyr)
ds <- read.delim("Input_Samples.txt")
ds <- subset(ds, ds$Year == 2017)
ds[ds == 0] <- NA

Generate correlation matrix heatmaps

The next thing I want to do is to make plots that illustrate the degree of correlation between sequential CPG sites for each of the four genes. In this plots I’m ignoring the fact that some individuals are represented by multiple samples, but I don’t think that is really a problem for making this point.

library(reshape2)
library(ggplot2)

## Write a little function to extract upper triangle from correlation matrix
  get_upper_tri <- function(cormat){
    cormat[lower.tri(cormat)] <- NA
    return(cormat)
  }

CRHmat <- as.matrix(ds[, 7:17])
CRHRmat <- as.matrix(ds[, 30:48])
GRmat <- as.matrix(ds[, 49:62])
FKmat <- as.matrix(ds[, 19:29])
mats<-list(CRHmat, CRHRmat, GRmat, FKmat)
labs<-c("CRH","CRHR1","GR","FKBP5")

for(i in 1:4){
    USE <- mats[[i]]

    corm <- round(cor(USE, use = "pairwise.complete.obs"), 2)
    cormUpper <- get_upper_tri(corm)
    
    ## Melt the correlation matrix
      melted_cormat <- melt(cormUpper, na.rm = TRUE)
      colnames(melted_cormat) <- c("CpG_Num1", "CpG_Num2", "value")
    
    ## Plot the heatmap
      gghm <- ggplot(data = melted_cormat, aes(CpG_Num2, CpG_Num1, fill = value)) +
        geom_tile(color = "white") +
        scale_fill_gradient2(low = "red", high = "blue", mid = "white",
            midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Pearson\nCorrelation") +
        theme_minimal() + 
        theme(axis.text.x = element_text(angle = 45, vjust = 1,
              size = 9, hjust = 1)) +
        coord_fixed()
      
      gghm2 <- gghm + 
        geom_text(aes(CpG_Num2, CpG_Num1, label = value), color = "white", size = 1.5) +
        theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank(),
          axis.ticks = element_blank(),
          legend.justification = c(1, 0),
          legend.position = c(0.6, 0.7),
          legend.direction = "horizontal")+
        guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
              title.position = "top", title.hjust = 0.5))
    
      assign(paste(labs[i],"gghm",sep="_"),gghm2)
    
}
CRH_gghm
CRHR1_gghm
GR_gghm
FKBP5_gghm

Plot overall methylation percentages for each gene

The methylation percentages seem to be very highly correlated among CpG sites within a single assay, both when looking just at consecutive sites (second diagonal row) and when looking at the whole assay wide level. The exceptions are GR -18 and GR -7 to -5, which seem to vary more or less independently from the other sites. Overall I think it suggests that it is probably reasonable to average over all of the sites in an assay (with possible exception of GR) and then to use that average as the main response for each individual. One problem with averaging, though, is that not all individuals have data for each CpG. Before averaging, we should check how methylation values differ across the different sites because very different values will make it trickier to do weighted averages when some data are missing.

  CRHmeth <- rep(NA, ncol(CRHmat))
  CRHsd <- CRHmeth
  for(i in 1:ncol(CRHmat)){
    CRHmeth[i] <- mean(na.omit(CRHmat[,i]))
    CRHsd[i] <- sd(na.omit(CRHmat[,i]))
    }

  CRHRmeth <- rep(NA, ncol(CRHRmat))
  CRHRsd <- CRHRmeth
  for(i in 1:ncol(CRHRmat)){
    CRHRmeth[i] <- mean(na.omit(CRHRmat[,i]))
    CRHRsd[i] <- sd(na.omit(CRHRmat[,i]))
  }
  
  GRmeth <- rep(NA, ncol(GRmat))
  GRsd <- GRmeth
  for(i in 1:ncol(GRmat)){
    GRmeth[i] <- mean(na.omit(GRmat[,i]))
    GRsd[i] <- sd(na.omit(GRmat[,i]))
  }
  
  FKmeth <- rep(NA, ncol(FKmat))
  FKsd <- FKmeth
  for(i in 1:ncol(FKmat)){
    FKmeth[i] <- mean(na.omit(FKmat[,i]))
    FKsd[i] <- sd(na.omit(FKmat[,i]))
  }
  
  plot(1,1,type="n",yaxs="i",xaxs="i",ylim=c(0,100),xlim=c(0,20),yaxt="n",xaxt="n",
       xlab="CpG Sites",ylab="Methylation",bty="n")
  axis(1,seq(-2,22,2),cex=0.9)
  axis(2,seq(-20,120,20),las=2)
  
  CRHo <- .1
  CRHRo <- -.1
  GRo <- .2
  FKo <- -.2
  
  r<-seq(1,20,1)
  
  for(i in 1:length(CRHmeth)){
    lines(rep(r[i],2)+CRHo,c(CRHmeth[i]+CRHsd[i],CRHmeth[i]-CRHsd[i]))
  }
  
  for(i in 1:length(GRmeth)){
    lines(rep(r[i],2)+GRo,c(GRmeth[i]+GRsd[i],GRmeth[i]-GRsd[i]))
  }
  
  for(i in 1:length(CRHRmeth)){
    lines(rep(r[i],2)+CRHRo,c(CRHRmeth[i]+CRHRsd[i],CRHRmeth[i]-CRHRsd[i]))
  }
  
  for(i in 1:length(FKmeth)){
    lines(rep(r[i],2)+FKo,c(FKmeth[i]+FKsd[i],FKmeth[i]-FKsd[i]))
  }
  
  lines(seq(1,length(CRHmeth),1)+CRHo,CRHmeth,lwd=2,col="orange")
  points(seq(1,length(CRHmeth),1)+CRHo,CRHmeth,pch=21,bg="orange")
  
  lines(seq(1,length(GRmeth),1)+GRo,GRmeth,lwd=2,col="lightblue")
  points(seq(1,length(GRmeth),1)+GRo,GRmeth,pch=21,bg="lightblue")
  
  lines(seq(1,length(CRHRmeth),1)+CRHRo,CRHRmeth,lwd=2,col="coral3")
  points(seq(1,length(CRHRmeth),1)+CRHRo,CRHRmeth,pch=21,bg="coral3")
  
  lines(seq(1,length(FKmeth),1)+FKo,FKmeth,lwd=2,col="chartreuse3")
  points(seq(1,length(FKmeth),1)+FKo,FKmeth,pch=21,bg="chartreuse3")
  
  legend("topright",c("Crh","Crhr1","GR","Fkbp5"),lwd=2,pch=21,bty="n",
         pt.bg=c("orange","coral3","lightblue","chartreuse3"),
         col=c("orange","coral3","lightblue","chartreuse3"))

The x-axis here doesn’t really mean anything since it is just showing position within the assay. I’ve overlaid all four genes to look at the variation, but remember that these are from totally different locatoins and the same x-axis value for the four genes here is meaningless. The point here is just to show that for each gene, some individual CpG sites have overall higher or lower average methylation values. That means that if we just do a raw average across all CpG sites, missing values at high methylation CpGs will lower the average more than missing values at low methylation CpG sites (and vice versa). A simple solution to this is to scale all the methylation values within each CpG site so that the average value is always 0 and one standard deviation is always 1 unit. Essentially what this would do is fill in missing values with the average (0) for that CpG. I’ve taken that approach below because it is easy…BUT, I think a better approach would probably be to fit this in a multi-level model with missing values imputed as parameters and then used to generate individual measures with uncertainty to use as predictors in the next level up. That is obviously a bit more complicated, but is probably something that I can think about doing. In the end it may not make too much difference since the number of zeros isn’t huge. I’m also not sure that one is necessarily more or less conservative since substituting average values could increase or decrease our chances of finding effects depending on the details of the data.

Scale and calculate assay average methylation

Here I am going back to the first data frame and scaling every individual column (after excluding 0s as above) so that the mean = 0 and SD = 1. Then I calculate average values for each individual gene, except that for GR I’m analyzing -18, -17 to -8, -7, and -6 to -5 as four separate groups because of the lack of correlation in the heatmaps above.

dss <- as.data.frame(scale(ds[,7:62]))
for(i in 1:nrow(dss)){
  dss$CRH[i] <- mean(na.omit(t(dss[i,1:12])))
  dss$CRHR1[i] <- mean(na.omit(t(dss[i,24:42])))
  dss$FK[i] <- mean(na.omit(t(dss[i,13:23])))
  dss$GRn17t8[i] <- mean(na.omit(t(dss[i,44:53])))
  dss$GRn6t5[i] <- mean(na.omit(t(dss[i,55:56])))
}

dss <- cbind(ds[,1:6],dss)

Plot values by groups

This is just a quick plot to visualize what the values look like at 1st and 3rd capture for control vs. dulled females and for each of the genes (with 4 comparisons for GR). This isn’t the formal analysis and I’m not running any stats here since we are planning to include initial color in formal models. I just want to see what the distribution of data look like and make sure nothing weird is going on before modeling.

library(scales)
plot(1,1,type="n",bty="n",yaxt="n",xaxt="n",xlab="",ylab="Methylation (sd)",xlim=c(0,21.5),ylim=c(-3,3))
axis(1,at=c(-5,1.5,4.5,9.5,12.5,17.5,20.5,30),labels=c("","Pre","Post","Pre","Post","Pre","Post",""))
axis(2,at=seq(-5,5,1))
abline(h=0,lty=2,col="gray40")
nogr <- list(subset(dss$CRH, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$CRH, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$CRH, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$CRH, dss$Capture == 3 & dss$Treatment == "Dulled"),
             subset(dss$CRHR1, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$CRHR1, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$CRHR1, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$CRHR1, dss$Capture == 3 & dss$Treatment == "Dulled"),
             subset(dss$FK, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$FK, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$FK, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$FK, dss$Capture == 3 & dss$Treatment == "Dulled"))
xs<-c(1,2,4,5,9,10,12,13,17,18,20,21)
gap<-0.17
cc<-rep(c("orange","lightblue"),8)
for(i in 1:length(xs)){
  points(rep(xs[i],length(nogr[[i]]))+runif(length(nogr[[i]]),min= -gap, max = gap-.05),nogr[[i]],pch=21,bg=alpha(cc[i],0.5))
  lines(rep(xs[i],2)+.33,c(mean(na.omit(nogr[[i]]))+sd(na.omit(nogr[[i]])),
                           mean(na.omit(nogr[[i]]))-sd(na.omit(nogr[[i]]))))
  points(xs[i]+.33,mean(na.omit(nogr[[i]])),pch=22,col=cc[i],bg="black",cex=1.2)
}
cs<-1.5
xs<- -2.5
text(3,xs,"Crh",cex=cs)
text(11,xs,"Crhr1",cex=cs)
text(19,xs,"Fkbp5",cex=cs)

legend("topright",c("Control","Dulled"),pch=21,pt.bg=c("orange","lightblue"),bty="n")


plot(1,1,type="n",bty="n",yaxt="n",xaxt="n",xlab="",ylab="Methylation (sd)",xlim=c(0,30),ylim=c(-3,3))
axis(1,at=c(-5,1.5,4.5,9.5,12.5,17.5,20.5,25.5,28.5,35),
     labels=c("","Pre","Post","Pre","Post","Pre","Post","Pre","Post",""))
axis(2,at=seq(-5,5,1))
abline(h=0,lty=2,col="gray40")
nogr <- list(subset(dss$GRn18, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$GRn18, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$GRn18, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$GRn18, dss$Capture == 3 & dss$Treatment == "Dulled"),
             subset(dss$GRn17t8, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$GRn17t8, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$GRn17t8, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$GRn17t8, dss$Capture == 3 & dss$Treatment == "Dulled"),
             subset(dss$GRn7, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$GRn7, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$GRn7, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$GRn7, dss$Capture == 3 & dss$Treatment == "Dulled"),
             subset(dss$GRn6t5, dss$Capture == 1 & dss$Treatment == "Control"),
             subset(dss$GRn6t5, dss$Capture == 1 & dss$Treatment == "Dulled"),
             subset(dss$GRn6t5, dss$Capture == 3 & dss$Treatment == "Control"),
             subset(dss$GRn6t5, dss$Capture == 3 & dss$Treatment == "Dulled"))
xs<-c(1,2,4,5,9,10,12,13,17,18,20,21,25,26,28,29)
gap<-0.17
cc<-rep(c("orange","lightblue"),8)
for(i in 1:length(xs)){
  points(rep(xs[i],length(nogr[[i]]))+runif(length(nogr[[i]]),min= -gap, max = gap-.05),nogr[[i]],pch=21,bg=alpha(cc[i],0.5))
  lines(rep(xs[i],2)+.33,c(mean(na.omit(nogr[[i]]))+sd(na.omit(nogr[[i]])),
                           mean(na.omit(nogr[[i]]))-sd(na.omit(nogr[[i]]))))
  points(xs[i]+.33,mean(na.omit(nogr[[i]])),pch=22,col=cc[i],bg="black",cex=1.2)
}
cs<-1.5
xs<- -2.5
text(3,xs,"GRn8",cex=cs)
text(11,xs,"GRn17t8",cex=cs)
text(19,xs,"GRn7",cex=cs)
text(27,xs,"GRn6t5",cex=cs)

legend("topright",c("Control","Dulled"),pch=21,pt.bg=c("orange","lightblue"),bty="n")

Nothing super obvious seems to be going on here in terms of changes, but this doesn’t yet account for initial coloration or initial methylation status. One weird thing that is visible here is that GRn7 (which was weird already in not being spatially correlated) has a very strange bimodal distribution. We saw some CpGs like this in the NGS data too. I’m not sure if it is some kind of artifact or a real case of where methylation tends to be bimodal for some reason. It is possible that more individual CpGs are like this too since I haven’t plotted these distributions for every individual site and it would get smoothed over by the averages. One other interesting thing is that some of the GR sites do seem to have an overall seasonal change in methylation, although it doesn’t seem like that slope varies by treatment.

Fit models

Now I’m going to fit the models that we proposed in the pre-registration for each of the four genes (fitting the 4 GR measures separately). The structure of these models is like this:

Methylation ~ Capture# x Treatment x Brightness + (1|Individual)

In the tables below each bold heading is one model and the coefficients, CI, and p-value are listed for each model with descriptions of the random effects below. I can make point estimates with confidence intervals for all differences or plot relationships (e.g. color vs. methylation at 3rd capture for dull vs. control), but it doesn’t look like there is much to see.

dss$Band <- as.factor(dss$Band)
dss$Capture <- as.factor(dss$Capture)
dss$BibB1 <- scale(dss$BibB1)

library(lme4)
library(lmerTest)
library(sjPlot)

CRH.m <- lmer(CRH ~ Capture*Treatment*BibB1 + (1|Band), data = dss)
CRHR1.m <- lmer(CRHR1 ~ Capture*Treatment*BibB1 + (1|Band), data = dss)
FK.m <- lmer(FK ~ Capture*Treatment*BibB1 + (1|Band), data = dss)

GRn18.m <- lmer(GRn18 ~ Capture*Treatment*BibB1 + (1|Band), data = dss)
GRn17t8.m <- lmer(GRn17t8 ~ Capture*Treatment*BibB1 + (1|Band), data = dss)
GRn7.m <- lmer(GRn7 ~ Capture*Treatment*BibB1 + (1|Band), data = dss)
GRn6t5.m <- lmer(GRn6t5 ~ Capture*Treatment*BibB1 + (1|Band), data = dss)

tab_model(CRH.m,CRHR1.m,FK.m,p.style="asterisk")
  CRH CRHR 1 FK
Predictors Estimates CI Estimates CI Estimates CI
(Intercept) 0.17 -0.12 – 0.45 -0.23 -0.52 – 0.06 0.15 -0.15 – 0.44
Capture 3 -0.08 -0.26 – 0.10 0.32 ** 0.14 – 0.51 -0.11 ** -0.18 – -0.03
Dulled -0.22 -0.62 – 0.19 0.26 -0.16 – 0.67 -0.13 -0.56 – 0.30
Bib B 1 0.19 -0.08 – 0.45 0.05 -0.22 – 0.32 -0.08 -0.35 – 0.20
Capture3:TreatmentDulled -0.03 -0.30 – 0.24 -0.19 -0.46 – 0.08 -0.10 -0.21 – 0.01
Capture3:BibB1 -0.07 -0.24 – 0.10 0.00 -0.17 – 0.17 -0.00 -0.07 – 0.07
TreatmentDulled:BibB1 -0.13 -0.55 – 0.28 0.02 -0.41 – 0.45 0.15 -0.28 – 0.59
Capture3:TreatmentDulled:BibB1 0.04 -0.24 – 0.31 -0.10 -0.38 – 0.17 -0.03 -0.14 – 0.08
Random Effects
σ2 0.12 0.12 0.02
Ï„00 0.61 Band 0.65 Band 0.80 Band
ICC 0.84 Band 0.85 Band 0.98 Band
Observations 118 118 118
Marginal R2 / Conditional R2 0.039 / 0.845 0.033 / 0.851 0.022 / 0.979
  • p<0.05   ** p<0.01   *** p<0.001 
tab_model(GRn18.m,GRn17t8.m,GRn7.m,GRn6t5.m,p.style="asterisk")
  G Rn 18 G Rn 17 t 8 G Rn 7 G Rn 6 t 5
Predictors Estimates CI Estimates CI Estimates CI Estimates CI
(Intercept) -0.06 -0.43 – 0.32 -0.21 -0.48 – 0.06 -0.17 -0.52 – 0.17 -0.66 *** -0.91 – -0.42
Capture 3 0.24 -0.07 – 0.55 0.14 -0.05 – 0.34 0.25 ** 0.10 – 0.41 0.94 *** 0.71 – 1.18
Dulled -0.19 -0.71 – 0.33 0.27 -0.12 – 0.65 0.11 -0.39 – 0.61 0.54 ** 0.19 – 0.89
Bib B 1 0.10 -0.25 – 0.46 -0.10 -0.35 – 0.16 0.07 -0.25 – 0.40 -0.08 -0.30 – 0.15
Capture3:TreatmentDulled 0.29 -0.17 – 0.74 0.09 -0.20 – 0.37 0.01 -0.21 – 0.24 -0.09 -0.44 – 0.26
Capture3:BibB1 0.11 -0.18 – 0.40 0.02 -0.16 – 0.20 -0.07 -0.22 – 0.08 0.03 -0.19 – 0.26
TreatmentDulled:BibB1 -0.14 -0.66 – 0.39 0.15 -0.24 – 0.55 -0.15 -0.66 – 0.36 0.16 -0.20 – 0.52
Capture3:TreatmentDulled:BibB1 -0.11 -0.55 – 0.34 0.01 -0.28 – 0.30 0.12 -0.12 – 0.35 -0.30 -0.66 – 0.06
Random Effects
σ2 0.28 0.13 0.08 0.21
Ï„00 0.74 Band 0.53 Band 1.01 Band 0.33 Band
ICC 0.72 Band 0.80 Band 0.93 Band 0.62 Band
Observations 100 118 116 118
Marginal R2 / Conditional R2 0.052 / 0.737 0.056 / 0.813 0.020 / 0.929 0.325 / 0.741
  • p<0.05   ** p<0.01   *** p<0.001

So…there doesn’t seem to be much going on. A few of the values change between 1st and 3rd capture but don’t differ between groups. There is a main effect of dulling for GRn6t5, but that seems be driven largely by the fact that the dulled group had higher pre and post treatment methylation values rather than any kind of different change between captures across the groups. Of course there are a lot of other things we can look at with these data, but for these main effects of coloring there doesn’t seem to be much happening. One other interesting thing here is that the ICC values for individuals seem to be very high. Perhaps that suggests that these values are really just not changing much on the time scale that we are looking at here. I’m plotting individual pairs of 1st to 3rd capture values (ignoring treatment for now) below to take a look at this.

Individual consistency in methylation

These are just some plots to ask how similar individuals are in their 1st and 3rd capture methylation scores for the different assays that we are analyzing above.

Note that I am plotting these on the standardized scale for convenience so that the mean is always 0 and SD is always 1 unit. In reality these different assays vary between each other in the average methylation (see plot above).

library(plyr)
dssj <- dss[,c("Band","Capture","Treatment","CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")]
dss1 <- subset(dssj,dssj$Capture == 1)
dss3 <- subset(dssj,dssj$Capture == 3)
dss3 <- dss3[,c("Band","CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")]
colnames(dss3) <- c("Band",paste(colnames(dss3)[2:8],".3",sep=""))

d.j <- join(dss1,dss3,"Band")

palette(c("lightblue","orange"))
p.comp<-function(x,y,cols,tit="GENE",X="Pre-Methylation",Y="Post-Methylation"){
  pal<-c()
  plot(1,1,type="n",bty="n",ylab=Y,xlab=X,xaxt="n",yaxt="n",xaxs="i",yaxs="i",
       xlim=c(-3,3),ylim=c(-3,3),main=tit)
  abline(lm(y~x))
  axis(1,seq(-5,5))
  axis(2,seq(-5,5),las=2)
  points(x,y,pch=21,bg=cols)
  text(2,-2.5,paste("r = ",round(cor(x,y,use="pairwise.complete.obs"),2)))
}
p.comp(d.j$CRH,d.j$CRH.3,d.j$Treatment,tit="Crh")
p.comp(d.j$CRHR1,d.j$CRHR1.3,d.j$Treatment,tit="Crhr1")
p.comp(d.j$FK,d.j$FK.3,d.j$Treatment,tit="Fkbp5")
p.comp(d.j$GRn18,d.j$GRn18.3,d.j$Treatment,tit="GRn18")
p.comp(d.j$GRn17t8,d.j$GRn17t8.3,d.j$Treatment,tit="GRn17t8")
p.comp(d.j$GRn7,d.j$GRn7.3,d.j$Treatment,tit="GRn7")
p.comp(d.j$GRn6t5,d.j$GRn6t5.3,d.j$Treatment,tit="GRn6t5")

So these plots match up pretty well with the very high ICC estimates from the mixed models above. It seems that percent methylation is very stable within individuals going from 1st to 3rd captures. Most of these are looking at averages across a bunch of CpGs so maybe it isn’t that surprising. I can certainly repeat any of these analyses at an individual CpG level and I expect there would be more variation there if only just from random noise, but remember that they tend to be all very highly correlated with each other within an assay. In one sense, this result is encouraging because it suggests a lot of confidence in the methylation scores themselves (they are clearly measuring something about individuals), but it may mean that we aren’t likely to see changes on this short time scale. It is also interesting to note that the CpG with a bimodal distribution (GRn7) is completely consistent here with maybe only two individuals switching from the low to the high modal group between captures.

Color Stressor Experiment

I’m going to run through the same sequence of steps described above for the 2018 color-stressor experiment here. I won’t describe as much since the logic is basically the same. The only differences here are that we have a 2x3 treatment (dull vs. control and stress vs. predator vs. control) and we have samples for 1st, 2nd, and 3rd capture.

Read and organize the data

Still excluding zeros. This time including only 2018 color-stress data.

ds <- read.delim("Input_Samples.txt")
ds <- subset(ds, ds$Year == 2018)
ds <- subset(ds, ds$Treatment != "Blood" & ds$Treatment != "Hippocampus" & ds$Treatment != "Hypothalamus")
ds[ds == 0] <- NA
for(i in 1:nrow(ds)){
  ds$Treat1[i] <- strsplit(as.character(ds$Treatment[i]),"_")[[1]][1]
  ds$Treat2[i] <- strsplit(as.character(ds$Treatment[i]),"_")[[1]][2]
}

Generate correlation matrix heatmaps

Same as in the first experiment above. The plots below look very similar overall to the correlation structure we saw in the 2017 experiment.

library(reshape2)
library(ggplot2)

## Write a little function to extract upper triangle from correlation matrix
  get_upper_tri <- function(cormat){
    cormat[lower.tri(cormat)] <- NA
    return(cormat)
  }

CRHmat <- as.matrix(ds[, 7:17])
CRHRmat <- as.matrix(ds[, 30:48])
GRmat <- as.matrix(ds[, 49:62])
FKmat <- as.matrix(ds[, 19:29])
mats<-list(CRHmat, CRHRmat, GRmat, FKmat)
labs<-c("CRH","CRHR1","GR","FKBP5")

for(i in 1:4){
    USE <- mats[[i]]

    corm <- round(cor(USE, use = "pairwise.complete.obs"), 2)
    cormUpper <- get_upper_tri(corm)
    
    ## Melt the correlation matrix
      melted_cormat <- melt(cormUpper, na.rm = TRUE)
      colnames(melted_cormat) <- c("CpG_Num1", "CpG_Num2", "value")
    
    ## Plot the heatmap
      gghm <- ggplot(data = melted_cormat, aes(CpG_Num2, CpG_Num1, fill = value)) +
        geom_tile(color = "white") +
        scale_fill_gradient2(low = "red", high = "blue", mid = "white",
            midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Pearson\nCorrelation") +
        theme_minimal() + 
        theme(axis.text.x = element_text(angle = 45, vjust = 1,
              size = 9, hjust = 1)) +
        coord_fixed()
      
      gghm2 <- gghm + 
        geom_text(aes(CpG_Num2, CpG_Num1, label = value), color = "white", size = 1.5) +
        theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank(),
          axis.ticks = element_blank(),
          legend.justification = c(1, 0),
          legend.position = c(0.6, 0.7),
          legend.direction = "horizontal")+
        guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
              title.position = "top", title.hjust = 0.5))
    
      assign(paste(labs[i],"gghm",sep="_"),gghm2)
    
}
CRH_gghm
CRHR1_gghm
GR_gghm
FKBP5_gghm

Plot overall methylation percentages for each gene

This is also plotted exactly the same as above. The methylation patterns here are shockingly similar to the 2017 data with an entirely separate dataset. I guess this is also encouraging for how well the method worked if nothing else…

  CRHmeth <- rep(NA, ncol(CRHmat))
  CRHsd <- CRHmeth
  for(i in 1:ncol(CRHmat)){
    CRHmeth[i] <- mean(na.omit(CRHmat[,i]))
    CRHsd[i] <- sd(na.omit(CRHmat[,i]))
    }

  CRHRmeth <- rep(NA, ncol(CRHRmat))
  CRHRsd <- CRHRmeth
  for(i in 1:ncol(CRHRmat)){
    CRHRmeth[i] <- mean(na.omit(CRHRmat[,i]))
    CRHRsd[i] <- sd(na.omit(CRHRmat[,i]))
  }
  
  GRmeth <- rep(NA, ncol(GRmat))
  GRsd <- GRmeth
  for(i in 1:ncol(GRmat)){
    GRmeth[i] <- mean(na.omit(GRmat[,i]))
    GRsd[i] <- sd(na.omit(GRmat[,i]))
  }
  
  FKmeth <- rep(NA, ncol(FKmat))
  FKsd <- FKmeth
  for(i in 1:ncol(FKmat)){
    FKmeth[i] <- mean(na.omit(FKmat[,i]))
    FKsd[i] <- sd(na.omit(FKmat[,i]))
  }
  
  plot(1,1,type="n",yaxs="i",xaxs="i",ylim=c(0,100),xlim=c(0,20),yaxt="n",xaxt="n",
       xlab="CpG Sites",ylab="Methylation",bty="n")
  axis(1,seq(-2,22,2),cex=0.9)
  axis(2,seq(-20,120,20),las=2)
  
  CRHo <- .1
  CRHRo <- -.1
  GRo <- .2
  FKo <- -.2
  
  r<-seq(1,20,1)
  
  for(i in 1:length(CRHmeth)){
    lines(rep(r[i],2)+CRHo,c(CRHmeth[i]+CRHsd[i],CRHmeth[i]-CRHsd[i]))
  }
  
  for(i in 1:length(GRmeth)){
    lines(rep(r[i],2)+GRo,c(GRmeth[i]+GRsd[i],GRmeth[i]-GRsd[i]))
  }
  
  for(i in 1:length(CRHRmeth)){
    lines(rep(r[i],2)+CRHRo,c(CRHRmeth[i]+CRHRsd[i],CRHRmeth[i]-CRHRsd[i]))
  }
  
  for(i in 1:length(FKmeth)){
    lines(rep(r[i],2)+FKo,c(FKmeth[i]+FKsd[i],FKmeth[i]-FKsd[i]))
  }
  
  lines(seq(1,length(CRHmeth),1)+CRHo,CRHmeth,lwd=2,col="orange")
  points(seq(1,length(CRHmeth),1)+CRHo,CRHmeth,pch=21,bg="orange")
  
  lines(seq(1,length(GRmeth),1)+GRo,GRmeth,lwd=2,col="lightblue")
  points(seq(1,length(GRmeth),1)+GRo,GRmeth,pch=21,bg="lightblue")
  
  lines(seq(1,length(CRHRmeth),1)+CRHRo,CRHRmeth,lwd=2,col="coral3")
  points(seq(1,length(CRHRmeth),1)+CRHRo,CRHRmeth,pch=21,bg="coral3")
  
  lines(seq(1,length(FKmeth),1)+FKo,FKmeth,lwd=2,col="chartreuse3")
  points(seq(1,length(FKmeth),1)+FKo,FKmeth,pch=21,bg="chartreuse3")
  
  legend("topright",c("Crh","Crhr1","GR","Fkbp5"),lwd=2,pch=21,bty="n",
         pt.bg=c("orange","coral3","lightblue","chartreuse3"),
         col=c("orange","coral3","lightblue","chartreuse3"))

Scale and calculate assay average methylation

Nothing different to see here. Just doing the same as above to prepare for plotting.

dss <- as.data.frame(scale(ds[,7:62]))
for(i in 1:nrow(dss)){
  dss$CRH[i] <- mean(na.omit(t(dss[i,1:12])))
  dss$CRHR1[i] <- mean(na.omit(t(dss[i,24:42])))
  dss$FK[i] <- mean(na.omit(t(dss[i,13:23])))
  dss$GRn17t8[i] <- mean(na.omit(t(dss[i,44:53])))
  dss$GRn6t5[i] <- mean(na.omit(t(dss[i,55:56])))
}

dss <- cbind(ds[,1:6],dss,ds[,63:64])

Plot values by groups

Here I have to diverge a bit since we have three captures and a 2x3 treatment structure. Because there are so many groups, I’m going to plot just the means with lines and error bars for each group.

  dsx <- dss
  gr <- c("CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")
  for(i in 1:7){
    dsx$GG <- dsx[,c(gr[i])]
    groups<-list(subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Control_Control"),
               subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Control_Stress"),
               subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Control_Predator"),
               subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Dull_Control"),
               subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Dull_Stress"),
               subset(dsx$GG, dsx$Capture == 1 & dsx$Treatment == "Dull_Predator"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Control_Control"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Control_Stress"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Control_Predator"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Dull_Control"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Dull_Stress"),
               subset(dsx$GG, dsx$Capture == 2 & dsx$Treatment == "Dull_Predator"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Control_Control"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Control_Stress"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Control_Predator"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Dull_Control"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Dull_Stress"),
               subset(dsx$GG, dsx$Capture == 3 & dsx$Treatment == "Dull_Predator"))
    mus <- rep(NA, length(groups))
    sds <- rep(NA, length(groups))
    for(k in 1:length(mus)){
      mus[k] <- mean(na.omit(groups[[k]]))
      sds[k] <- sd(na.omit(groups[[k]]))
    }
    
    plot(1,1,type="n",xaxt="n",yaxt="n",bty="n",yaxs="i",xaxs="i",ylab="Methylation (sd)",
        xlab="Capture number",main=gr[i],ylim=c(-2,2),xlim=c(.5,4.5))
    abline(h=0,lty=2,col="gray40")
    axis(1,at=c(-2,1,2,3,10),labels=c("","1","2","3",""))
    axis(2,seq(-10,10,1),las=2)
    
    offs<-c(-.05,-.03,-.01,.01,.03,.05)
    cols<-c("darkgoldenrod1","darkorange2","coral3","lightblue","chartreuse3","darkorchid4")
    for(j in 1:6){
      lines(c(1,2,3)+offs[j],c(mus[j],mus[j+6],mus[j+12]),lwd=2,col=cols[j])
        lines(rep(1,2)+offs[j],c(mus[j]+sds[j],mus[j]-sds[j]))
        lines(rep(2,2)+offs[j],c(mus[j+6]+sds[j+6],mus[j+6]-sds[j+6]))
        lines(rep(3,2)+offs[j],c(mus[j+12]+sds[j+12],mus[j+12]-sds[j+12]))
      points(c(1,2,3)+offs[j],c(mus[j],mus[j+6],mus[j+12]),pch=21,bg=cols[j])
    }
    legend("bottomright",lwd=2,col=cols,c("Control Control","Control Stress","Control Pred",
                           "Dull Control","Dull Stress","Dull Pred"),pch=21,pt.bg=cols,bty="n")
    
  }

From these plots it doesn’t seem like much is going on. We see the same seasonal shift in GRn6t5, which is interesting, but again it doesn’t differ by treatment group. I colored all the control groups in warm colors and dull groups in cool colors and it also doesn’t seem like there is much evidence for differences in those groups at either time point. Next I’ll fit the models like the ones fit for the color manipulation experiment above.

Fit models

I’m doing this just the same as above except that we have the two different treatments. The full model is like this:

Methylation ~ Capture x Treat1 x Treat2 x BibB1 + (1|Band)

dss$Band <- as.factor(dss$Band)
dss$Capture <- as.factor(dss$Capture)
dss$BibB1 <- scale(dss$BibB1)
dss$Treat1 <- as.factor(dss$Treat1)
dss$Treat2 <- as.factor(dss$Treat2)

library(lme4)
library(lmerTest)
library(sjPlot)

CRH.m <- lmer(CRH ~ Capture*Treat1*Treat2 + (1|Band), data = dss)
CRHR1.m <- lmer(CRHR1 ~ Capture*Treat1*Treat2 + (1|Band), data = dss)
FK.m <- lmer(FK ~ Capture*Treat1*Treat2 + (1|Band), data = dss)

GRn18.m <- lmer(GRn18 ~ Capture*Treat1*Treat2 + (1|Band), data = dss)
GRn17t8.m <- lmer(GRn17t8 ~ Capture*Treat1*Treat2 + (1|Band), data = dss)
GRn7.m <- lmer(GRn7 ~ Capture*Treat1*Treat2 + (1|Band), data = dss)
GRn6t5.m <- lmer(GRn6t5 ~ Capture*Treat1*Treat2 + (1|Band), data = dss)

tab_model(CRH.m,CRHR1.m,FK.m,p.style="asterisk")
  CRH CRHR 1 FK
Predictors Estimates CI Estimates CI Estimates CI
(Intercept) 0.16 -0.30 – 0.62 -0.01 -0.62 – 0.60 -0.11 -0.63 – 0.41
Capture 2 -0.20 -0.44 – 0.05 -0.12 -0.33 – 0.09 -0.01 -0.14 – 0.12
Capture 3 -0.02 -0.25 – 0.20 -0.03 -0.25 – 0.18 -0.06 -0.19 – 0.07
Dull -0.16 -0.86 – 0.54 0.19 -0.74 – 1.12 0.32 -0.48 – 1.11
Predator -0.24 -0.90 – 0.43 0.02 -0.86 – 0.90 0.05 -0.70 – 0.80
Stress -0.04 -0.70 – 0.63 -0.26 -1.12 – 0.60 -0.20 -0.93 – 0.54
Capture2:Treat1Dull 0.02 -0.39 – 0.42 0.04 -0.28 – 0.37 0.10 -0.09 – 0.29
Capture3:Treat1Dull 0.12 -0.27 – 0.50 -0.09 -0.41 – 0.23 0.08 -0.11 – 0.27
Capture2:Treat2Predator 0.10 -0.28 – 0.48 -0.10 -0.39 – 0.19 0.10 -0.08 – 0.28
Capture3:Treat2Predator 0.00 -0.36 – 0.37 -0.14 -0.45 – 0.17 0.13 -0.06 – 0.32
Capture2:Treat2Stress -0.02 -0.45 – 0.40 0.03 -0.27 – 0.34 -0.03 -0.21 – 0.16
Capture3:Treat2Stress -0.20 -0.61 – 0.21 0.06 -0.25 – 0.37 0.01 -0.18 – 0.20
Treat1Dull:Treat2Predator 0.74 -0.24 – 1.72 0.18 -1.12 – 1.47 0.11 -1.00 – 1.21
Treat1Dull:Treat2Stress -0.33 -1.31 – 0.65 0.39 -0.89 – 1.68 -0.15 -1.25 – 0.94
Capture2:Treat1Dull:Treat2Predator -0.09 -0.71 – 0.52 -0.01 -0.44 – 0.43 -0.17 -0.43 – 0.10
Capture3:Treat1Dull:Treat2Predator 0.00 -0.56 – 0.57 0.20 -0.26 – 0.65 -0.19 -0.47 – 0.09
Capture2:Treat1Dull:Treat2Stress 0.22 -0.40 – 0.85 0.07 -0.39 – 0.53 -0.07 -0.34 – 0.21
Capture3:Treat1Dull:Treat2Stress 0.17 -0.43 – 0.77 0.40 -0.06 – 0.86 -0.04 -0.31 – 0.23
Random Effects
σ2 0.06 0.06 0.02
Ï„00 0.59 Band 1.10 Band 0.82 Band
ICC 0.91 Band 0.95 Band 0.97 Band
Observations 127 172 176
Marginal R2 / Conditional R2 0.098 / 0.915 0.056 / 0.953 0.063 / 0.976
  • p<0.05   ** p<0.01   *** p<0.001 
tab_model(GRn18.m,GRn17t8.m,GRn7.m,GRn6t5.m,p.style="asterisk")
  G Rn 18 G Rn 17 t 8 G Rn 7 G Rn 6 t 5
Predictors Estimates CI Estimates CI Estimates CI Estimates CI
(Intercept) -0.37 -0.95 – 0.22 -0.09 -0.51 – 0.32 0.27 -0.28 – 0.83 -0.14 -0.64 – 0.36
Capture 2 0.33 * 0.05 – 0.60 0.46 * 0.10 – 0.83 0.07 -0.07 – 0.20 0.41 * 0.03 – 0.80
Capture 3 0.48 ** 0.20 – 0.75 0.57 ** 0.21 – 0.94 0.08 -0.06 – 0.22 0.64 ** 0.26 – 1.03
Dull 0.35 -0.55 – 1.24 -0.10 -0.73 – 0.53 0.26 -0.59 – 1.11 -0.32 -1.08 – 0.45
Predator -0.08 -0.92 – 0.77 -0.51 -1.11 – 0.09 -0.28 -1.09 – 0.52 -0.72 -1.44 – 0.00
Stress 0.17 -0.73 – 1.06 0.29 -0.30 – 0.87 -0.92 * -1.72 – -0.11 0.08 -0.62 – 0.79
Capture2:Treat1Dull -0.37 -0.79 – 0.05 0.05 -0.50 – 0.60 -0.11 -0.32 – 0.10 0.06 -0.52 – 0.64
Capture3:Treat1Dull -0.37 -0.78 – 0.05 -0.01 -0.56 – 0.54 -0.03 -0.24 – 0.18 -0.00 -0.59 – 0.58
Capture2:Treat2Predator 0.00 -0.38 – 0.39 -0.04 -0.55 – 0.47 -0.09 -0.28 – 0.10 0.23 -0.30 – 0.77
Capture3:Treat2Predator -0.04 -0.45 – 0.37 -0.07 -0.62 – 0.47 -0.01 -0.22 – 0.19 0.23 -0.34 – 0.80
Capture2:Treat2Stress 0.17 -0.28 – 0.62 -0.16 -0.69 – 0.37 -0.05 -0.26 – 0.17 0.01 -0.55 – 0.57
Capture3:Treat2Stress 0.08 -0.37 – 0.53 -0.24 -0.77 – 0.28 0.00 -0.21 – 0.21 0.26 -0.30 – 0.82
Treat1Dull:Treat2Predator -0.14 -1.46 – 1.18 0.63 -0.25 – 1.51 -0.52 -1.71 – 0.67 1.00 -0.06 – 2.06
Treat1Dull:Treat2Stress 0.13 -1.18 – 1.43 -0.41 -1.28 – 0.45 0.17 -1.02 – 1.35 0.22 -0.83 – 1.26
Capture2:Treat1Dull:Treat2Predator 0.30 -0.34 – 0.93 0.02 -0.74 – 0.78 0.20 -0.09 – 0.49 -0.26 -1.06 – 0.54
Capture3:Treat1Dull:Treat2Predator 0.59 -0.08 – 1.27 -0.33 -1.11 – 0.46 0.06 -0.24 – 0.36 0.07 -0.76 – 0.91
Capture2:Treat1Dull:Treat2Stress 0.05 -0.58 – 0.68 -0.28 -1.06 – 0.49 0.19 -0.11 – 0.50 -0.55 -1.37 – 0.27
Capture3:Treat1Dull:Treat2Stress -0.25 -0.88 – 0.38 -0.54 -1.32 – 0.23 0.28 -0.02 – 0.58 -0.90 * -1.72 – -0.08
Random Effects
σ2 0.10 0.18 0.03 0.20
Ï„00 0.95 Band 0.34 Band 0.95 Band 0.56 Band
ICC 0.90 Band 0.65 Band 0.97 Band 0.74 Band
Observations 151 176 172 176
Marginal R2 / Conditional R2 0.058 / 0.908 0.192 / 0.721 0.121 / 0.977 0.180 / 0.785
  • p<0.05   ** p<0.01   *** p<0.001

These model tables are harder to interpret without plotting point estimates with confidence intervals or estimating specific differences. I can go through and do that for each group, but overall it doesn’t look like there are any real differences between the treatment groups at capture 2 or 3. Just like above the ICC values for ID are very high, which makes it look like methylation is not changing much within individuals over this time frame.

Individual consistency in methylation

I’m making these basically the same as before except now that there are three captures I’m comparing 1 to 2 and 2 to 3 for each panel. I could also compare 1 to 3…but doesn’t seem totally necessary for now.

library(plyr)
dssj <- dss[,c("Band","Capture","Treatment","CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")]
dss1 <- subset(dssj,dssj$Capture == 1)
dss2 <- subset(dssj,dssj$Capture == 2)
dss2 <- dss2[,c("Band","CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")]
colnames(dss2) <- c("Band",paste(colnames(dss2)[2:8],".2",sep=""))
dss3 <- subset(dssj,dssj$Capture == 3)
dss3 <- dss3[,c("Band","CRH","CRHR1","FK","GRn18","GRn17t8","GRn7","GRn6t5")]
colnames(dss3) <- c("Band",paste(colnames(dss3)[2:8],".3",sep=""))

d.j <- join(dss1,dss2,"Band")
d.j <- join(d.j,dss3,"Band")

palette(c("lightblue","orange"))
p.comp<-function(x,y,cols,tit="GENE",X="Pre-Methylation",Y="Post-Methylation"){
  pal<-c()
  plot(1,1,type="n",bty="n",ylab=Y,xlab=X,xaxt="n",yaxt="n",xaxs="i",yaxs="i",
       xlim=c(-3,3),ylim=c(-3,3),main=tit)
  abline(lm(y~x))
  axis(1,seq(-5,5))
  axis(2,seq(-5,5),las=2)
  points(x,y,pch=21,bg=cols)
  text(2,-2.5,paste("r = ",round(cor(x,y,use="pairwise.complete.obs"),2)))
}
p.comp(d.j$CRH,d.j$CRH.2,d.j$Treatment,tit="Crh 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$CRH.2,d.j$CRH.3,d.j$Treatment,tit="Crh 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$CRHR1,d.j$CRHR1.2,d.j$Treatment,tit="Crhr1 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$CRHR1.2,d.j$CRHR1.3,d.j$Treatment,tit="Crhr1 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$FK,d.j$FK.2,d.j$Treatment,tit="Fkbp5 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$FK.2,d.j$FK.3,d.j$Treatment,tit="Fkbp5 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$GRn18,d.j$GRn18.2,d.j$Treatment,tit="GRn18 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$GRn18.2,d.j$GRn18.3,d.j$Treatment,tit="GRn18 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$GRn17t8,d.j$GRn17t8.2,d.j$Treatment,tit="GRn17t8 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$GRn17t8.2,d.j$GRn17t8.3,d.j$Treatment,tit="GRn17t8 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$GRn7,d.j$GRn7.2,d.j$Treatment,tit="GRn7 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$GRn7.2,d.j$GRn7.3,d.j$Treatment,tit="GRn7 2 to 3",X="Capture 2",Y="Capture 3")
p.comp(d.j$GRn6t5,d.j$GRn6t5.2,d.j$Treatment,tit="GRn6t5 1 to 2",X="Capture 1",Y="Capture 2")
p.comp(d.j$GRn6t5.2,d.j$GRn6t5.3,d.j$Treatment,tit="GRn6t5 2 to 3",X="Capture 2",Y="Capture 3")