inst/doc/RT-PCR-Script-ddCt.R

In ddCt: The ddCt Algorithm for the Analysis of Quantitative Real-Time PCR (qRT-PCR)

### R code from vignette source 'RT-PCR-Script-ddCt.Rnw'

###################################################
### code chunk number 1: setup
###################################################
options(width=50)


###################################################
### code chunk number 2: ddCtInstall (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install(ddCt)


###################################################
### code chunk number 3: specifyFiles
###################################################
load.path = "."
save.path = "."


###################################################
### code chunk number 4: The Data
###################################################

Ct.data.file.name       <-  c("")

sample.annotation.file.name  <- NULL
column.for.grouping          <- NULL
useOnlySamplesFromAnno       <- FALSE


###################################################
### code chunk number 5: transform
###################################################
#new.gene.names    <- c("Gen1"="Gene A",
#		        "Gen4"="Gene B",
#			"Gen5"="Gene C",
#			"Gen2"="HK1",
#			"Gen3"="HK2")
new.gene.names    <- NULL
new.sample.names  <- NULL


###################################################
### code chunk number 6: housekeeping
###################################################
name.referenz.sample <- c("Sample2")
name.referenz.gene   <- c("Gene3")


###################################################
### code chunk number 7: threshold
###################################################
Threshold.for.Ct <- 40


###################################################
### code chunk number 8: mean/median
###################################################
TYPEOFCALCULATION <- "mean"


###################################################
### code chunk number 9: eff (eval = FALSE)
###################################################
## EFFICIENCIES <- c("Gene A"=1.9,"Gene B"=1.8,"HK1"=2,"Gene C"=2,"HK2"=2)
## EFFICIENCIES.ERROR <- c("Gene A"=0.01,"Gene B"=0.1,"HK1"=0.05,"Gene C"=0.01,"HK2"=0.2)


###################################################
### code chunk number 10: effizienzen
###################################################
EFFICIENCIES <- NULL


###################################################
### code chunk number 11: eff.error
###################################################
EFFICIENCIES.ERROR <- NULL


###################################################
### code chunk number 12: PlotKind
###################################################
TheKindOfPlot <- c("level","Ct")


###################################################
### code chunk number 13: remaining
###################################################
#REMAIN
names.of.the.genes.REMAIN.in.Output   <- NULL
names.of.the.samples.REMAIN.in.Output <- NULL
#NOT
names.of.the.genes.NOT.in.Output      <- NULL
names.of.the.samples.NOT.in.Output    <- NULL


###################################################
### code chunk number 14: grouping
###################################################
GROUPINGBYSAMPLES <- TRUE


###################################################
### code chunk number 15: plot0
###################################################
own.plot.for.each.object     <- TRUE


###################################################
### code chunk number 16: plot1
###################################################
GroupingForPlot <- NULL
GroupingColor   <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3","#FF7F00")


###################################################
### code chunk number 17: plot2
###################################################
CUTOFF <- NULL


###################################################
### code chunk number 18: plot3
###################################################
BREWERCOL <- c("Set3","Set1","Accent","Dark2","Spectral","PuOr","BrBG")


###################################################
### code chunk number 19: plot4
###################################################
LEGENDE <- TRUE


###################################################
### code chunk number 20: ttest1
###################################################
perform.ttest <- FALSE


###################################################
### code chunk number 21: ttest
###################################################
column.for.ttest <- NULL


###################################################
### code chunk number 22: ttest3
###################################################
column.for.pairs     <- NULL      ## ansonsten NULL


###################################################
### code chunk number 23: ttest4
###################################################
names.of.the.samples.REMAIN.in.ttest <- names.of.the.samples.REMAIN.in.Output
names.of.the.samples.NOT.in.ttest    <- names.of.the.samples.NOT.in.Output



###################################################
### code chunk number 24: RT-PCR-Script-ddCt.Rnw:382-384
###################################################
names.of.the.samples.REMAIN.in.cor <- names.of.the.samples.REMAIN.in.Output
names.of.the.samples.NOT.in.cor    <- names.of.the.samples.NOT.in.Output 


###################################################
### code chunk number 25: lib
###################################################
library(Biobase)
library(lattice)
library(RColorBrewer)
library(ddCt)


###################################################
### code chunk number 26: runastest
###################################################
## in case key parameters were not specified, the script runs in test mode
testMode <- FALSE
if (length(Ct.data.file.name) == 1 & Ct.data.file.name[1]=="")
  testMode <- TRUE
if(testMode) {
  Ct.data.file.name <- "Experiment1.txt"
  load.path <- system.file("extdata", package="ddCt")
  name.referenz.sample <- c("Sample2")
  name.referenz.gene   <- c("Gene3","Gene2")
}


###################################################
### code chunk number 27: readData
###################################################
datadir <- function(x) file.path(load.path, x)
savedir <- function(x) file.path(save.path,x)
file.names <- datadir(Ct.data.file.name)
sho <- paste(gsub(".txt","",Ct.data.file.name),collapse="_")

warningFile <- savedir(paste("warning.output",sho,".txt",sep=""))

CtData <- SDMFrame(file.names)

if (!is.null(Threshold.for.Ct)){
 A <- Ct(CtData)>= Threshold.for.Ct
 coreData(CtData)$Ct[A] <- NA
}

if (! is.null(new.gene.names))   CtData[,2] <- new.gene.names[CtData[,2]]
if (! is.null(new.sample.names)) CtData[,1] <- new.sample.names[CtData[,1]]



if(! is.null(sample.annotation.file.name)){
   info <- datadir(sample.annotation.file.name)
   sampleInformation <- read.AnnotatedDataFrame(info,header=TRUE,sep="\t", row.names=NULL)
 }else{
   sampleInformation <- new("AnnotatedDataFrame",
                            data=data.frame(Sample=sampleNames(CtData)),
			    varMetadata=data.frame(labelDescription=c("unique identifier"),row.names=c("Sample")))
 }

 
 

if(useOnlySamplesFromAnno && !is.null(sample.annotation.file.name)){
  A <- CtData[,"Sample"] %in% pData(sampleInformation)[,"Sample"]
  warning( paste("Es werden folgende Samples entfernt:",paste(unique(CtData[!A,"Sample"]),collapse=", ")))
  CtData <- CtData[A,]
 }


if (is.null(EFFICIENCIES)){
  result <-      ddCtExpression(CtData,
                      calibrationSample=name.referenz.sample,
                      housekeepingGenes=name.referenz.gene,
                      type=TYPEOFCALCULATION,
                      sampleInformation=sampleInformation,
                      warningStream=warningFile)
  } else{
  result <- ddCtwithE(CtData,
                      calibrationSample=name.referenz.sample,
                      housekeepingGenes=name.referenz.gene,  
                      efficiencies=EFFICIENCIES,
                      efficiencies.error=EFFICIENCIES.ERROR,
                      type=TYPEOFCALCULATION,
                      sampleInformation=sampleInformation,
                      warningStream=warningFile)
  }
   
save(result,file=savedir(paste("Result",sho,".RData",sep="")))


###################################################
### code chunk number 28: write (HTML)Tables
###################################################
htmlName   <- paste("HTML",sho,sep="_")
tablesName <- paste("Tables",sho,sep="_")
getDir <- function(dir, ...) {
  if(!file.exists(dir)) {
    dir.create(dir,...)
  }
  return(dir)
}


html.path <- getDir(file.path(save.path,htmlName))
table.path <- getDir(file.path(save.path, tablesName))

if(!is.null(sample.annotation.file.name) & !is.null(column.for.grouping))
 {result<- result[,order(pData(result)[,column.for.grouping])]}

elistWrite(result,file=savedir(paste("allValues",sho,".csv",sep="")))
EE1 <- assayData(result)$exprs
FF1 <- assayData(result)$level.err
if(!is.null(GroupingForPlot)) 
  GFP1 <- pData(result)[,GroupingForPlot]

Ct  <- round(assayData(result)$Ct,2)
lv  <- round(EE1,2)

write.table(cbind(t(EE1),t(FF1)),file=file.path(table.path,"LevelPlusError.txt"),sep="\t",col.names=NA)
write.table(lv,file=file.path(table.path,"level_matrix.txt"),sep = "\t", col.names = NA)
write.table(Ct,file=file.path(table.path,"Ct_matrix.txt"),sep="\t", col.names = NA)

write.htmltable(cbind(rownames(lv),lv),title="Level",file=file.path(html.path,"level"))
write.htmltable(cbind(rownames(Ct),Ct),title="Ct",file=file.path(html.path,"Ct"))


if(is.null(EFFICIENCIES)){
 dCtValues  <- round(assayData(result)$dCt,2)
 ddCtValues <- round(assayData(result)$ddCt,2)
 write.table(dCtValues,file=file.path(table.path,"dCt_matrix.txt"),  sep="\t", col.names=NA)
 write.table(ddCtValues,file=file.path(table.path,"ddCt_matrix.txt"),sep="\t",col.names=NA)

 write.htmltable(cbind(rownames(dCtValues),dCtValues)  ,title="dCt",file=file.path(html.path,"dCt"))
 write.htmltable(cbind(rownames(ddCtValues),ddCtValues),title="ddCt",file=file.path(html.path,"ddCt"))
}



###################################################
### code chunk number 29: plot Level
###################################################
for(KindOfPlot in TheKindOfPlot){
if(KindOfPlot=="level"){
  EE1 <- assayData(result)$exprs
  FF1 <- assayData(result)$level.err
  theTitle <- "Level"
}

if(KindOfPlot=="ddCt"){
  EE1 <- assayData(result)$ddCt
  FF1 <- assayData(result)$ddCt.err
  theTitle <- "ddCt"
}

if(KindOfPlot=="dCt"){
  EE1 <- assayData(result)$dCt
  FF1 <- assayData(result)$dCt.err
  theTitle <- "dCt"
}

if(KindOfPlot=="Ct"){
  EE1 <- assayData(result)$Ct
  FF1 <- assayData(result)$Ct.err
  theTitle <- "Ct"
}



################
## order the set
################

if (!is.null(new.gene.names))
 {EE2  <- EE1[match(new.gene.names,rownames(EE1)),,drop=FALSE]
  FF2  <- FF1[match(new.gene.names,rownames(EE1)),,drop=FALSE]
 } else{
  EE2 <- EE1
  FF2 <- FF1
 }

if (!is.null(new.sample.names))
 {EE  <- EE2[,match(new.sample.names,colnames(EE2)),drop=FALSE]
  FF  <- FF2[,match(new.sample.names,colnames(EE2)),drop=FALSE]
  if(!is.null(GroupingForPlot)) GFP<- GFP1[match(new.sample.names,colnames(EE2))]
 } else{
  EE <- EE2
  FF <- FF2
  if(!is.null(GroupingForPlot)) GFP<- GFP1
 }

###################
## Reducing the set
###################


if (!is.null(names.of.the.genes.REMAIN.in.Output)){
   Gred <- (rownames(EE) %in% names.of.the.genes.REMAIN.in.Output)
}else{
   Gred <- !(rownames(EE) %in% names.of.the.genes.NOT.in.Output)
}


if (!is.null(names.of.the.samples.REMAIN.in.Output)){
   Sred <- (colnames(EE) %in% names.of.the.samples.REMAIN.in.Output)
}else{
   Sred <- !(colnames(EE) %in% names.of.the.samples.NOT.in.Output)
}

EEN <- EE[Gred,Sred,drop=FALSE]
FFN <- FF[Gred,Sred,drop=FALSE]
if(!is.null(GroupingForPlot)) GFPN <- as.factor(as.character(GFP[Sred]))

############
# the color
############

COLORS <- c()
for (i in 1:length(BREWERCOL))
  COLORS <- c(COLORS,brewer.pal(brewer.pal.info[BREWERCOL[i],]$maxcolors,BREWERCOL[i]))
if(GROUPINGBYSAMPLES){
  THECO  <- COLORS[1:sum(Sred)]
} else {
  THECO  <- COLORS[1:sum(Gred)]
}

##########
# plotting
##########
getFileName <- function(prefix="Ct", name) {
  sprintf("%s_Result_%s.pdf", prefix, sho)
}
if(own.plot.for.each.object){
 pdf(w=15,h=15,file=savedir(paste(theTitle,"Result",sho,".pdf",sep="")))
if(GROUPINGBYSAMPLES){
  for(k in 1:dim(EEN)[1]){
    EENN <- EEN[k,,drop=FALSE]
    FFNN <- FFN[k,,drop=FALSE]
    barploterrbar(EENN,EENN-FFNN,EENN+FFNN,barcol=THECO,legend=LEGENDE,columnForDiffBars=GROUPINGBYSAMPLES,theCut=CUTOFF,ylab=theTitle)
  }}else{
    for(k in 1:dim(EEN)[2]){
      EENN <- EEN[,k,drop=FALSE]
      FFNN <- FFN[,k,drop=FALSE]
      barploterrbar(EENN,EENN-FFNN,EENN+FFNN,barcol=THECO,legend=LEGENDE,columnForDiffBars=GROUPINGBYSAMPLES,theCut=CUTOFF,ylab=theTitle)
    }}
 if(GROUPINGBYSAMPLES & !is.null(GroupingForPlot)){
   for(k in 1:dim(EEN)[1]){
     EENN <- EEN[k,,drop=FALSE]
     FFNN <- FFN[k,,drop=FALSE]
     barploterrbar(EENN,EENN-FFNN,EENN+FFNN,barcol=GroupingColor,legend=LEGENDE,columnForDiffBars=GROUPINGBYSAMPLES,theCut=CUTOFF,ylab=theTitle,las=2,names.arg=colnames(EENN),main=rownames(EENN),groups=GFPN)
   }}
 dev.off()
}else{
  barploterrbar(EEN,EEN-FFN,EEN+FFN,barcol=THECO,legend=LEGENDE,las=2,columnForDiffBars=GROUPINGBYSAMPLES,theCut=CUTOFF,ylab=theTitle)
  dev.copy(pdf,w=15,h=15,file=savedir(paste(theTitle,"Result",sho,".pdf",sep="")))
  dev.off()}
}
result.bySample <- errBarchart(result)
print(result.bySample)
pdf(getFileName("errbarplot_bySample_", sho), width=15, height=15)
print(result.bySample)
dev.off()
result.byDetector <- errBarchart(result, by="Detector")
print(result.byDetector)
pdf(getFileName("errbarplot_byDetector_", sho), width=15, height=15)
print(result.byDetector)
dev.off()


###################################################
### code chunk number 30: Korrelation
###################################################
A <- name.referenz.gene 
if (length(A) > 1) { 
 if (!is.null(names.of.the.samples.REMAIN.in.cor)){
   corRed <- (rownames(pData(result)) %in% names.of.the.samples.REMAIN.in.cor)
 }else{
   corRed <- !(rownames(pData(result)) %in% names.of.the.samples.NOT.in.cor)
 }
 result2 <- assayData(result[,corRed])$Ct

 K <- combn(1:length(A),2)
 U <- ncol(K)
 par(mfrow=c(ceiling(sqrt(U)),ceiling(sqrt(U))))
 for (i in 1:U){
  Gen1 <- A[K[1,i]] # der Name des ersten Housekeepinggenes
  Gen2 <- A[K[2,i]] # der Name des zweiten Housekeepinggenes
  BART <- cor(result2[Gen1,],result2[Gen2,],use="pairwise.complete.obs")
  if(!is.na(BART)) {
    plot(result2[Gen1,],result2[Gen2,],xlab=Gen1,ylab=Gen2,pch="*",col="red", main=paste("Correlation:",round(BART,3)))
  } else {
    gen1.allna <- all(is.na(result2[Gen1,]))
    warn <- sprintf("Correlation does not exist, since %s in all Samples are 'Undetermined'\n",ifelse(gen1.allna, Gen1, Gen2))

    plot.new()
    text(0.5,0.5, warn)
  }
  }} else{
   if (!is.null(names.of.the.samples.REMAIN.in.cor)){
   corRed <- (rownames(pData(result)) %in% names.of.the.samples.REMAIN.in.cor)
 }else{
   corRed <- !(rownames(pData(result)) %in% names.of.the.samples.NOT.in.cor)
 }
 result2 <- assayData(result[,corRed])$Ct
  plot(result2[A,],pch="*",col="red", main="Expression HK Gene")
 }


###################################################
### code chunk number 31: TTest
###################################################

myFoldChange <- function(x){
  x <- as.numeric(x)
  if(length(x) == 2 ){
    resu <- 2^(x[1]-x[2])
  }else if( length(x) == 1  ){
    resu <- 2^x
  }else{
    stop("unexpected situation in myFoldChange")
  }
  return(resu)
}


if(perform.ttest){
ttestName   <- paste("tTests",sho,sep="_")
ttest.path <- getDir(file.path(save.path, ttestName))

if (!is.null(names.of.the.samples.REMAIN.in.ttest)){
  ttestRed <- (rownames(pData(result)) %in% names.of.the.samples.REMAIN.in.ttest)
}else{
  ttestRed <- !(rownames(pData(result)) %in% names.of.the.samples.NOT.in.ttest)
}


result3 <- result[,ttestRed]


daten  <- assayData(result3)$dCt # der t-Test wird immer mit den normalisierten Werten gemacht
if( ! column.for.ttest %in% colnames(pData(result3)) )
  stop(paste(" did not find :", column.for.ttest,": in pData ",sep=""))
faktor <- as.character(pData(result3)[,column.for.ttest])
mmm <- nlevels(as.factor(faktor))
if( mmm == 1 ){
  stop( " found only a single group for t-test ")
}else if( mmm == 2 ){
  aa <- matrix(c(1,2), ncol=1) # aa = die paarweisen Vergleiche
}else{
  aa <- combn(1:mmm,2)
}

 
for (k in 1:ncol(aa)){
  
  Groups  <- levels(as.factor(faktor))[aa[,k]]
  subs    <- faktor %in%  Groups
  datenS  <- daten[,subs]
  faktorS <- as.factor(faktor[subs])  # hier wird gewaerleistet dass der neue Faktor genau zwei Elemente hat
  
  if( ! is.null(column.for.pairs) ){
    if( ! column.for.pairs %in% colnames(pData(result)) )
      stop(paste(" did not find :", column.for.pairs,": in pData ",sep=""))
    paarungS <- as.character(pData(result3)[,column.for.pairs])
    paarungS <- paarungS[subs]
    wenigerRes <- 1
    optTest <- "paired "
  }else{
    wenigerRes <- 0
    optTest <- ""
  }



  res  <- matrix(NA,nrow=nrow(datenS),ncol=8-wenigerRes)
  res2 <- matrix(NA,nrow=nrow(datenS),ncol=2-wenigerRes)

  for (i in 1:nrow(datenS)){ 
    a  <- datenS[i,]
    b  <- is.na(a)
    cc <- a[!b]
    d  <- faktorS[!b]
    if( ! is.null(column.for.pairs) ){ 
      paar <- as.character(paarungS[!b])
      ## restrict to valid pair data only
      validPaarItems <- paar[duplicated(paar)]
      valid <- which( as.character(paar) %in% validPaarItems )
      paar <- paar[valid]
      cc <- cc[valid]
      d <- d[valid]
      if(all(table(d) >1)) {
        sel1 <- which(as.character(d) == Groups[1])
        sel2 <- which(as.character(d) == Groups[2])
        stopifnot( all( paar[sel1][order(paar[sel1])] ==paar[sel2][order(paar[sel2])] )   )
        group1 <- cc[sel1][order(paar[sel1])]
        group2 <- cc[sel2][order(paar[sel2])]
        ff <- t.test(x=group1, y=group2, paired=TRUE)
        ff2<- wilcox.test(x=group1, y=group2, paired=TRUE)
      }else{
        ff <- NULL
        ff2 <- NULL
      }
    }else{
      if(all(table(d)>1)) {
        ff <- t.test(cc~d)
        ff2<- wilcox.test(cc~d)
      }else{
        ff  <- NULL
        ff2 <- NULL
      }
    }
    if( !is.null(ff) ){
      res[i,]  <-c(signif(ff$statistic),
	           signif(ff$p.value),
		   ff2$statistic,
		   signif(ff2$p.value),
		   myFoldChange(ff$estimate),
                   ff$parameter["df"],
                   ff$estimate)          ## 2 Stueck oder 1 Stueck (paired)
      
      res2[i,] <-c(names(ff$estimate))  ## 2 Stueck oder 1 Stueck (paired)
     }
  }
  AllGe      <- rownames(assayData(result)$exprs)
  theHKGenes <- rep("",length(AllGe))
  theHKGenes[AllGe %in% name.referenz.gene] <- "X"
  gg <- cbind(AllGe,theHKGenes,res)
  
   myColnames <- c("Name",
		   "Housekeeping Gene",
                   paste("statistic(", optTest,"t.test)", sep=""),
                   paste("pvalue(",    optTest,"t.test)", sep=""),
                   paste("statistic(", optTest,"Wilcox)", sep=""),
                   paste("pvalue(",    optTest,"Wilcox)", sep=""),
                   "foldChange",
                   "degreeOfFreedom"
                   )
  
  if( ! is.null(column.for.pairs) ){
    Mr1 <- res2[,1]
    extraName <- unique(Mr1[!is.na(Mr1)])
    if( length(extraName) < 1 ){
      extraName <- "fehlenUnklar"
    }
    myColnames <- c(myColnames, extraName) 
  }else{
    Mr1 <- res2[,1]
    Mr2 <- res2[,2]
    stopifnot(length(unique(Mr1[!is.na(Mr1)]))==1 & length(unique(Mr2[!is.na(Mr2)]))==1   )
    myColnames <- c(myColnames,
                    unique(Mr1[!is.na(Mr1)]),
                    unique(Mr2[!is.na(Mr2)])
                    ) 
  }
  
  colnames(gg) <- myColnames
  pVSpalte <- paste("pvalue(",    optTest,"t.test)", sep="")
  gg <-gg[order(as.numeric(gg[,pVSpalte])),]
  
  SAVED <- paste("ttest",Groups[1],Groups[2],sep="")
  write.table(gg,file=file.path(ttest.path, paste(SAVED,"txt",sep=".")),sep="\t",row.names=FALSE)
  write.htmltable(gg,file=file.path(ttest.path,SAVED))
}
}


###################################################
### code chunk number 32: BoxplotsHousekkepingGenes
###################################################
if(perform.ttest){
 if (!is.null(names.of.the.samples.REMAIN.in.ttest)){
   ttestRed <- (rownames(pData(result)) %in% names.of.the.samples.REMAIN.in.ttest)
 }else{
   ttestRed <- !(rownames(pData(result)) %in% names.of.the.samples.NOT.in.ttest)
 }
result3 <- result[,ttestRed]
daten   <- assayData(result3)$Ct  # es geht ja hier um die Expressionen der
				 # Housekeeping Gene vor Normalisierung
faktor	<- as.character(pData(result3)[,column.for.ttest])
mmm     <- levels(as.factor((faktor)))
N <- length(mmm)

daten2 <- daten[name.referenz.gene,,drop=FALSE] 
BoxPl <- list()
for (i in 1:N){
 BoxPl[[i]] <- t(daten2[,mmm[i]==faktor])
}

res <- list()
for(i in 1:length(name.referenz.gene)){
 A <- lapply(BoxPl, function(x) x[,i])
 names(A) <- rep(name.referenz.gene[i], N)
 res      <- c(res,A)
}

theColor <- 2 + (1:N)
pdf(file=savedir(paste("HKGenesPerGroup_",sho,".pdf",sep="")),w=15,h=15)
boxplot(res,col=theColor,main="Ct expression of housekeeping genes per group")
dev.off()
}


###################################################
### code chunk number 33: hh
###################################################
if(file.exists(warningFile)){
  bart   <- unlist(read.delim(warningFile,as.is=TRUE,header=FALSE))
  fehler <- sapply(bart, function(y) gsub("simpleWarning in withCallingHandlers\\(\\{: ","",y))
}


###################################################
### code chunk number 34: fehler
###################################################
if(file.exists(warningFile))
  fehler


###################################################
### code chunk number 35: session
###################################################
toLatex(sessionInfo())