Residual Disease Paper
========================================================
Exploring Default PCR Problems
-------------------------------------------------------


by Shelley Herbrich

*May 13, 2013*


```{r options, echo=FALSE}

opts_chunk$set(tidy=TRUE, message=FALSE, warning=FALSE)

```


## 1 Executive Summary

### 1.1 Introduction

Currently, PCR vendor software internally chooses a default threshold for each target gene on a plate (in this case, 18S, FABP4, and ADH1B). On revisiting the data obtained from Washington, we noticed that for a few of the plates the threshold chosen for FABP4 is quite low. The problem is that in going from one plate to another, these  cutoffs are largely stable for 18S and ADH1B, but not for FABP4. 

We suspect that the cutoffs for FABP4 may be being dragged out of place by some outlier wells, as might be the case  with one or two problematic samples. 

### 1.2 Methods

For this analysis, we focus primarily on the first plate run on Plate 2. Some of the samples were flagged as problematic (and were later rerun), but the threshold for FABP4 was applied to all samples, including those that were retained. 

We suspect a better threshold for FABP4 would be near the same level as those for ADH1B and 18S. We interpolate the revised Ct values based on this proposed threshold.  

### 1.3 Results

It appears that the software chooses the gene-specific threshold based on the first encountered curve for that particular target. This results in significantly underestimated thresholds when there are poor quality wells on a given plate. The implications of this change aren't subtle. The Ct values we would infer for FABP4 with both thresholds are different, and the mean difference is about 5.5 units (or about 2^5.5 = 46-fold). 

Therefore, we need to develop an alternative quantification method for our PCR data that is independent of the software-defined threshold values. 

## 3 Assessing the Impact of Incorrect Thresholds

```{r libraries}
library(RColorBrewer)
```
```{r loadPCRFiles}
load(file.path("RDataObjects","rawPCRWash_Unfiltered.RData"))
load(file.path("RDataObjects","PCRResults.RData"))
```

```{r extractAnnot}
#extract annotation info
rawPCRPlate2 <- rawPCRWash[which(rawPCRWash$Plate=="Plate.2"),]
wells <- rawPCRPlate2$Well
genes <- rawPCRPlate2$Target.Name
thresholds <- unique(rawPCRPlate2$Ct.Threshold)
names(thresholds) <- unique(genes)
```

```{r defineColors}
colors <- array(NA,length(genes))
colors[which(genes=="ADH1B")] <- rep(brewer.pal(9, "OrRd")[3:8],4)[1:23]
colors[which(genes=="FABP4")] <- rep(brewer.pal(9, "Purples")[3:8],4)[1:23]
colors[which(genes=="18S")] <- rep(brewer.pal(9, "Greens")[3:8],4)[1:23]
```

```{r deltaRnPlot}
par(mar=c(5.1,5.1,4.1,3.1))
plot(NA,xlim=c(1,40),ylim=c(-5,1.1), xlab="Cycle Number", ylab="",yaxt="n",main="Plate 2")
axis(2, at=seq(-5,2,by=2), labels=10^seq(-5,2,by=2),las=2, cex.axis=0.7)
mtext("Delta Rn", side=2, at=-2, line=4, las=0, cex=1)

for(i1 in 1:length(wells)){
  wn <- wells[i1]
  tempRn <- log10(rawPCRPlate2[which(rawPCRPlate2$Well==wn),grep("Cycle",colnames(rawPCRPlate2))])
  missing <- is.na(tempRn)
  points(c(1:40),tempRn,col=colors[i1],type="l")
}

#add thresholds
abline(h=log10(thresholds))
mtext(names(thresholds), side=4, at=log10(thresholds)+c(0.1,0,-0.1), line=0.5, las=2, cex=.75)

#add proposed threshold
otherFABP4 <- c(1.674,1.674,1.631,1.589)
estimatedThreshold <- mean(otherFABP4)
abline(h=log10(estimatedThreshold),col="red",lty=2)
mtext("FABP4", side=2, at=log10(estimatedThreshold), line=0.5, las=2, cex=.75,col="red")
```


```{r estimateFABP4Ct}
goi <- "FABP4"  
fabp4Wells <- wells[which(genes==goi)]
defaultCt <- as.numeric(rawPCRPlate2$Ct[which(genes==goi)])
mean(defaultCt,na.rm=TRUE)

estimatedCt <- array(NA,length(fabp4Wells))
for(i1 in 1:length(fabp4Wells)){
  wellID <- fabp4Wells[i1]
  subRn <- as.numeric(rawPCRPlate2[which(rawPCRPlate2$Well==wellID),grep("Cycle",colnames(rawPCRPlate2))])
  ub <- which(subRn>estimatedThreshold)[1]
  m <- (subRn[ub]-subRn[ub-1])
  b <- subRn[ub-1]
  estimatedCt[i1] <- (ub-1)+(estimatedThreshold-b)/m
}
mean(estimatedCt,na.rm=TRUE)
```

```{r plotProposedAndCurrent}
plot(defaultCt, ylim=c(min(c(defaultCt,estimatedCt),na.rm=TRUE),40), ylab="Number of Cycles",main="FABP4 PCR Results")
abline(h=40,lty=2)
points(c(1:23)[!is.na(defaultCt)],defaultCt[!is.na(defaultCt)],type="l")

points(estimatedCt,pch=8,col="red")
imputeCt <- estimatedCt
imputeCt[is.na(imputeCt)] <- 40
points(imputeCt,type="l",col="red")

legend(15,17, c("default Ct", "estimated Ct"), cex=0.8, col=c("black","red"), pch=c(1,8), bty="n")
```

```{r filterPoorQuality} 
#flagging bad rows
goodRow <- rep(TRUE, nrow(rawPCRWash))
## Plate Plate.1 '- 2nd set 021613' all fine
## Plate Plate.2 '021213', Rerun samples
goodRow[(rawPCRWash$Sample.Name == "W37") & (rawPCRWash$Plate == "Plate.2")] <- FALSE
goodRow[(rawPCRWash$Sample.Name == "W32") & (rawPCRWash$Plate == "Plate.2")] <- FALSE
goodRow[(rawPCRWash$Sample.Name == "W20") & (rawPCRWash$Plate == "Plate.2")] <- FALSE
## Plate Plate.3 '021613'
goodRow[(rawPCRWash$Well == "A6") & (rawPCRWash$Plate == "Plate.3")] <- FALSE
goodRow[(rawPCRWash$Well == "B6") & (rawPCRWash$Plate == "Plate.3")] <- FALSE
goodRow[(rawPCRWash$Well == "F3") & (rawPCRWash$Plate == "Plate.3")] <- FALSE
goodRow[(rawPCRWash$Well == "G2") & (rawPCRWash$Plate == "Plate.3")] <- FALSE
## Plate Plate.4 '021813'; all 80-4 (rerun)
goodRow[(rawPCRWash$Well == "C2") & (rawPCRWash$Plate == "Plate.4")] <- FALSE
goodRow[(rawPCRWash$Well == "C3") & (rawPCRWash$Plate == "Plate.4")] <- FALSE
goodRow[(rawPCRWash$Well == "C5") & (rawPCRWash$Plate == "Plate.4")] <- FALSE
goodRow[(rawPCRWash$Well == "C6") & (rawPCRWash$Plate == "Plate.4")] <- FALSE
## Plate Plate.5 '3rd batch 201313'
goodRow[(rawPCRWash$Well == "F6") & (rawPCRWash$Plate == "Plate.5")] <- FALSE
## Plate W6 '4th set 021413' all fine
## Plate W7 'second machine' row B (sample 81-2) may be weak overall
rawPCRWash <- rawPCRWash[goodRow,]
```

```{r sampleQPCR,echo=FALSE,eval=FALSE}
sampleQPCR <- function(x){
  tempNo <- array(NA,length(x$Well))
  for(i1 in 1:length(x$Well)){
    subDRn <- cbind(1:40, t(x[i1,grep("Cycle",colnames(x))]))
    m <- pcrfit(subDRn, cyc=1, fluo=2, model=l3)
    to <- tryCatch(takeoff(m)$top, error=function(e) NA)
    border <- c(0,0)
    if(to < x[i1,"Baseline.End"] & !is.na(to)) border <-c(x[i1,"Baseline.End"]-to,0)
    if(x$Target.Name[i1]=="18S"){
      sw <- tryCatch(sliwin(m,border=border,wsize=4:6,plot=FALSE), error=function(e) NULL)
    }else{
      sw <- tryCatch(sliwin(m,border=border,plot=FALSE), error=function(e) NULL) 
    }
    if(!is.null(sw)) tempNo[i1] <- unlist(sw["init"])
  }
  targetMeans <- sapply(split(tempNo,x$Target.Name),mean,na.rm=TRUE)
  qpcr <- log2(targetMeans["FABP4"])-log2(targetMeans["18S"])
  qpcr
}
```

```{r sampleCt}
sampleCt <- function(x){
  geneMeans <- sapply(split(x,factor(x$Target.Name)),function(y) mean(as.numeric(y$Ct),na.rm=TRUE))
  geneThresh <- unique(x$Ct.Threshold)
  names(geneThresh) <- unique(x$Target.Name)
  rv <- geneMeans["18S"]-geneMeans["FABP4"]
  if(geneThresh["FABP4"] < 1){
    tmp <- x[which(x$Target.Name=="FABP4"),]
    estimatedCt <- array(NA,nrow(tmp))
    for(i1 in 1:length(estimatedCt)){
      subRn <- as.numeric(tmp[i1,grep("Cycle",colnames(tmp))])
      ub <- which(subRn>estimatedThreshold)[1]
      m <- (subRn[ub]-subRn[ub-1])
      b <- subRn[ub-1]
      estimatedCt[i1] <- (ub-1)+(estimatedThreshold-b)/m
    }
    rv <- c(rv,geneMeans["18S"]-mean(estimatedCt,na.rm=TRUE))
  }else{
    rv <- c(rv,NA)
  }
  rv
}
```

```{r getCt}
ct <- NULL
for(i1 in 1:length(unique(rawPCRWash$Plate))){
  tmp <- rawPCRWash[which(rawPCRWash$Plate==unique(rawPCRWash$Plate)[i1]),]
  tmpCt <- t(sapply(split(tmp,factor(tmp$Sample.Name)),sampleCt))
  ct <- rbind(ct,tmpCt)
}
```
```{r arrowPlot}
plot(ct[PCRResults$Sample.Name[which(PCRResults$Source=="Washington")],1], PCRResults$FABP4[which(PCRResults$Source=="Washington")],main="FABP4 PCR Results",xlab="Delta Ct",ylab="Alternative Quantification",pch=21,bg="grey",cex=0.75)
points(ct[PCRResults$Sample.Name[which(PCRResults$Source=="Washington")],2], PCRResults$FABP4[which(PCRResults$Source=="Washington")], col="red",pch=8)
arrows(ct[PCRResults$Sample.Name[which(PCRResults$Source=="Washington")],1],PCRResults$FABP4[which(PCRResults$Source=="Washington")],x1=ct[PCRResults$Sample.Name[which(PCRResults$Source=="Washington")],2],col="pink",length=.15)
legend("topleft", c("default Ct", "estimated Ct"), cex=0.8, col=c("grey","red"), pch=c(1,8), bty="n")
```

## Appendix

```{r getLocation}
getwd()
```

```{r sessionInfo}
sessionInfo()
```