Skip to content
Permalink
Browse files

update references to supplementary figures in comments to final figur…

…e order
  • Loading branch information...
fhalbritter committed May 20, 2016
1 parent a1abf65 commit ff9ce5563759489de3df6f139b796c2bfc71c001
Showing with 33 additions and 30 deletions.
  1. +14 −12 analysis/methBench_additions.R
  2. +1 −1 analysis/methBench_fig5.R
  3. +18 −17 analysis/methBench_suppFigs.R
@@ -10,7 +10,7 @@ message("=== ADDITIONAL DATA AND ANALYSES ===")


##############################################################################
############ CLONAL BISULFITE SEQUENCING (Supplementary Figure 7) ############
############ CLONAL BISULFITE SEQUENCING (Supplementary Figure 8) ############
##############################################################################

dClonalBS <- read.table("data/results_validation_ClonalBS.txt", sep="\t", header=TRUE, comment.char="", quote="", check.names=FALSE)
@@ -25,17 +25,19 @@ dOther <- lapply(split(dOther,dOther$variable), function(x) { rownames(x) <- x$r
dClonalBS <- dClonalBS*100.0


### REPLICATES COMPARISON (Supplementary Figure 7.a) ###
### REPLICATES COMPARISON (Supplementary Figure 8.a) ###

svgPlot("clonal_reps", 3, 3)
svgPlot("clonal_reps", 6, 6)
par(mar=c(4,4,1,1)+0.1)
par(xpd=NA)
x <- c(dClonalBS["region_07",],dClonalBS["region_08",])
y <- c(dClonalBS["mandatory_7_b",],dClonalBS["mandatory_8_b",])
plot(as.numeric(x), as.numeric(y), xlab="ClonalBS 1", ylab="ClonalBS 2", pch=16, col=rgb(0.1,0.1,0.1,0.3), xlim=c(0,100), ylim=c(0,100), bty="l")
plot(as.numeric(x), as.numeric(y), xlab="ClonalBS 1", ylab="ClonalBS 2", pch=16, col=c(rep(rgb(0.5,0.1,0.1,0.3),ncol(dClonalBS)),rep(rgb(0.1,0.1,0.5,0.3),ncol(dClonalBS))), xlim=c(0,100), ylim=c(0,100), bty="l")
text(as.numeric(x), as.numeric(y), paste(c(rep("r7",ncol(dClonalBS)),rep("r8",ncol(dClonalBS))),colnames(dClonalBS)))
legend("topleft", paste0("r = ", round(cor(as.numeric(x), as.numeric(y), use="pairwise", method="pearson"),2)), bty="n")
dev.off()

### CORRELATION HEATMAP (Supplementary Figure 7.b) ###
### CORRELATION HEATMAP (Supplementary Figure 8.b) ###

tmpSelRegions <- c("region_07","region_08")
dCor <- dOther[intersect(names(dOther),datasetsByType$absolute)]
@@ -49,7 +51,7 @@ CairoSVG("results_analysis/plots/clonal_cors_reg78", width=13, height=13, points
pheatmap(data.matrix(dCor), ,col=colorRampPalette(rev(brewer.pal(5,"Oranges")),bias=0.1,space="Lab")(20*5),breaks=seq(0,1,0.01),border_color="white",cellwidth=22,cellheight=22, cluster_cols=TRUE, cluster_rows=TRUE, treeheight_col=0, number_color="black", display_numbers=TRUE)
dev.off()

### DEVIATION FROM CORRIDOR (Supplementary Figure 7.c) ###
### DEVIATION FROM CORRIDOR (Supplementary Figure 8.c) ###

# only proceed with selected regions and samples:
dClonalBS <- dClonalBS[selRegions, selSamples]
@@ -67,7 +69,7 @@ dClonalDev$crc <- substr(dClonalDev$variable,1,5)
message("ClonalBS: mean = ", round(mean(unlist(dClonalBS), na.rm=T),3), ", sd = ", round(sd(unlist(dClonalBS), na.rm=T),3), ", min = ", round(min(unlist(dClonalBS), na.rm=T),3), ", max =", round(max(unlist(dClonalBS), na.rm=T),3))
message("ClonalBS: quantiles (0.1,0.25,0.5,0.75,0.9) = ", paste(round(quantile(unlist(dClonalBS),c(0.1,0.25,0.5,0.75,0.9), na.rm=T),3), collapse=", "))
message("ClonalBS: deviation - mean abs. = ", round(mean(abs(dClonalDev$deviation), na.rm=T),3), ", median = ", round(median(abs(dClonalDev$deviation), na.rm=T),3), ", directional = ", round(mean(dClonalDev$deviation, na.rm=T),3))

# absolute and directional deviation plots:
tmp <- lsaData[ lsaData$sampleName%in%selSamples & lsaData$regionName%in%selRegions & lsaData$datasetName%in%datasetsByType$absolute,c("datasetName","deviationCorridor")]
tmp$datasetName <- datasetTable[tmp$datasetName,"prettyLabel"]
@@ -78,7 +80,7 @@ svgPlotGG(p, "clonal_box", 9, 19)
p <- ggplot(aggregate(deviationCorridor~datasetName,allDevData,mean), aes(x=datasetName, y=deviationCorridor, fill=deviationCorridor)) + geom_bar(colour="#333333", stat="identity") + defaultPlotTheme(flipX=FALSE) + xlab(NULL) + ylab("Directional deviation") + coord_flip() + scale_fill_gradient2(low=brewer.pal(3,"BrBG")[3],high=brewer.pal(3,"BrBG")[1],guide=F,space="Lab") + geom_hline(aes(yintercept=0),colour="#333333") + ylim(-5,5)
svgPlotGG(p, "clonal_box_dir", 6.5, 19)

### SCATTER PLOT PANEL (Supplementary Figure 7.d) ###
### SCATTER PLOT PANEL (Supplementary Figure 8.d) ###

dScatter <- dcast(rbind(data.frame(melt(cbind("region"=rownames(dClonalBS[selRegions,selSamples]),dClonalBS[selRegions,selSamples])),"L1"="ClonalBS 1"),melt(lapply(dOther, function(x) { cbind("region"=rownames(x[selRegions,selSamples]),x[selRegions,selSamples]) }))), L1~region+variable, value.var="value")
rownames(dScatter) <- dScatter[,1]
@@ -99,7 +101,7 @@ svgPlotGG(d, "clonal_scatters", 12, 7, units="in")


######################################################################################
############ LOW-INPUT TITRATION SERIES (Supplementary Figures 10 and 11) ############
############ LOW-INPUT TITRATION SERIES (Supplementary Figures 11 and 12) ############
######################################################################################

dataDir <- "data/"
@@ -180,22 +182,22 @@ dLowInput$corridorUpper <- apply(dLowInput,1,function(x) consensus$upper[x["regi
dLowInput$corridorDeviation <- ifelse(dLowInput$value>=dLowInput$corridorLower & dLowInput$value<=dLowInput$corridorUpper,0,ifelse(dLowInput$value<dLowInput$corridorLower, dLowInput$value-dLowInput$corridorLower, dLowInput$value-dLowInput$corridorUpper))


### DEVIATION FROM CONSENSUS (Supplementary Figures 10.b and 11.b) ###
### DEVIATION FROM CONSENSUS (Supplementary Figures 11.b and 12.b) ###

yRange <- ceiling(max(c(dLowInput$ownTargetDeviation,dLowInput$corridorDeviation),na.rm=TRUE))
for(curSeries in unique(dLowInput$series)) {
p <- ggplot(dLowInput[dLowInput$series==curSeries,], aes(x=as.numeric(quantity), y=corridorDeviation, fill=factor(prettyLabel))) + geom_point() + stat_summary(fun.y = mean, geom = 'ribbon', fun.ymax = max, fun.ymin = min, .alpha = 0.05, alpha = 0.5) + defaultPlotTheme(flipX=TRUE) + xlab("DNA quantity [ng]") + ylab("Deviation from corridor") + scale_fill_manual(values=plotColLookup$prettyLabel, guide=FALSE) + scale_x_continuous(breaks=1:length(levels(dLowInput$quantity)),labels=levels(dLowInput$quantity)) + facet_grid(prettyLabel~., drop=FALSE) + geom_hline(y_intercept=0) + theme(axis.line.x=element_blank()) + ylim(-yRange,yRange)
svgPlotGG(p, paste0("low_input_dev_consensus_",curSeries), 8,24)
}

### DEVIATION FROM OWN TARGET (Supplementary Figures 10.c and 11.c) ###
### DEVIATION FROM OWN TARGET (Supplementary Figures 11.c and 12.c) ###

for(curSeries in unique(dLowInput$series)) {
p <- ggplot(dLowInput[!is.na(dLowInput$ownTargetDeviation) & dLowInput$series==curSeries,], aes(x=as.numeric(quantity), y=ownTargetDeviation, fill=factor(prettyLabel))) + geom_point() + stat_summary(fun.y = mean, geom = 'ribbon', fun.ymax = max, fun.ymin = min, .alpha = 0.05, alpha = 0.5) + defaultPlotTheme(flipX=TRUE) + xlab("DNA quantity [ng]") + ylab("Deviation from own target value") + scale_fill_manual(values=plotColLookup$prettyLabel, guide=FALSE) + scale_x_continuous(breaks=1:length(levels(dLowInput$quantity)),labels=levels(dLowInput$quantity)) + facet_grid(prettyLabel~., drop=FALSE) + geom_hline(y_intercept=0) + theme(axis.line.x=element_blank()) + ylim(-yRange,yRange)
svgPlotGG(p, paste0("low_input_dev_owntarget_",curSeries), 8,24)
}

### STATUS (Supplementary Figures 10.a and 11.a) ###
### STATUS (Supplementary Figures 11.a and 12.a) ###

for(curSeries in unique(dLowInput$series)) {
p <- ggplot(dLowInputSuccessRate[dLowInputSuccessRate$series==curSeries,], aes(x=quantity, y=successRate)) + geom_bar(stat="identity") + defaultPlotTheme(flipX=TRUE) + xlab("DNA quantity [ng]") + ylab("Success rate (successful / attempted measurements)") + facet_grid(prettyLabel~., drop=FALSE) + geom_hline(y_intercept=0) + theme(axis.line.x=element_blank())
@@ -122,7 +122,7 @@ aucsNoise <- aggregate(AUC~dataset,data=aucData,mean,na.rm=T)
aucsNoise <- aucsNoise[order(aucsNoise$AUC),]


### SUPPLEMENTARY PANEL OF ROC CURVES (Supplementary Figure 13) ###
### SUPPLEMENTARY PANEL OF ROC CURVES (Supplementary Figure 14) ###

svgPlot("roc_small_all", 2*3,2.2*7, pointsize=18)
par(mfrow=c(7,3))
@@ -9,10 +9,10 @@ message("=== SUPPLEMENTARY FIGURE AND DATA ===")



##### HEATMAP-LIKE OVERVIEW OF AVAILABLE DATASETS (Fig. S1) #####
##### HEATMAP-LIKE OVERVIEW OF AVAILABLE DATASETS (Supplementary Figure 1) #####


message("Data overview table (Figure S1)...")
message("Data overview table (Supplementary Figure 1)...")

# calculate an overview of available assays:
tmp <- lsaData
@@ -66,10 +66,10 @@ svgPlotGG(d, "hm_regions_per_assay_sample", 26, 12)



### MDS PLOTS (Fig. Sx) ###
### MDS PLOTS (Supplementary Figure 13) ###

# generate a panel of MDS plots per dataset as a quick overview of the data:
message("Multi-dimensional scaling plots (Figure Sx)...")
message("Multi-dimensional scaling plots (Supplementary Figure 13)...")

svgPlot("mds_plots", 12, 10, pointsize=13)
par(mar=c(3,3,4,1))
@@ -113,9 +113,9 @@ dev.off()



### FULL PANEL OF BEESWARM PLOTS (Fig. S4) ###
### FULL PANEL OF BEESWARM PLOTS (Supplementary Figure 3) ###

message("Beeswarm plots panel for all assays (Figure S4.a)...")
message("Beeswarm plots panel for all assays (Supplementary Figure 3.a)...")

selRegions <- coreRegionNames
selSamples <- setdiff(sampleNames,sampleNames[grep("Titration",sampleNames)])
@@ -125,7 +125,7 @@ for(curRegion in selRegions) {

# determine which datasets define each consensus corridor:

message("Contribution to consensus heatmap (Figure S4.b)...")
message("Contribution to consensus heatmap (Supplementary Figure 3.b)...")

#selSamples <- sampleNames
inConsensus <- aggregate(assayGroup~datasetName+regionName+sampleName,consensusFreq[consensusFreq$sampleName%in%selSamples & consensusFreq$regionName%in%selRegions,],function(x){ 1 })
@@ -171,15 +171,16 @@ tmp$color <- plotColLookup$assayGroup[as.character(tmp[,"Assay type"])]
tmp$color[tmp$variable=="Not contributing to corridor"] <- makeTransparent(tmp$color[tmp$variable=="Not contributing to corridor"],0.25)

# and visualize the summary as a barplot:
message("Contribution to consensus bar plots (Figure S4.c)...")
message("Contribution to consensus bar plots (Supplementary Figure 3.c)...")
d <- ggplot(tmp, aes(x=`Assay type`,y=value,fill=color)) + geom_bar(width=0.5,colour="#333333",size=0.25,stat="identity") + geom_hline(yintercept=0) + ylim(0,1.2*max(tmp$height,na.rm=TRUE))+ xlab(NULL) + geom_text(aes(y=height,label=percent),angle=90,size=2,hjust=0) + ylab("Number of measurements") + scale_fill_identity() + defaultPlotTheme(flipX=TRUE,fontSize=8)
svgPlotGG(d, "contribution_to_consensus_summary", 4, 6)


### FRESH-FROZEN VS FFPE ###

### FRESH-FROZEN VS FFPE (Supplementary Figure 6) ###

message("Fresh-frozen vs FFPE (Figures Sxxx)...")

message("Fresh-frozen vs FFPE (Supplementary Figure 6)...")

selSamples <- sampleNamesByType$FrozenFFPE
selRegions <- regionNames
@@ -253,10 +254,10 @@ svgPlotGG(d,"performance_frozenffpe_bias_quant_v2", 18.3, 9, units="cm")



### FULL SCATTER PLOT PANELS OF ABSOLUTE AND RELATIVE ASSAYS (Supplementary Fig. 2 and 8) ###
### FULL SCATTER PLOT PANELS OF ABSOLUTE AND RELATIVE ASSAYS (Supplementary Figures 2 and 9) ###

selSamples <- sampleNames
message("Scatter plot panels (Supplementary Figures 2 and 8)...")
message("Scatter plot panels (Supplementary Figures 2 and 9)...")

scatterComparisons <- list("absVsAbs"=c("absolute","absolute"),"relVsAll"=c("relative","absoluteAndRelative"))
for(scatterCmp in scatterComparisons) {
@@ -322,9 +323,9 @@ for(scatterCmp in scatterComparisons) {



##### LINEAR MODELS FOR CAUSES OF DEVIATION (Supplementary Figure 4) #####
##### LINEAR MODELS FOR CAUSES OF DEVIATION (Supplementary Figure 5) #####

message("Causes of deviation (Supplementary Figure 4)...")
message("Causes of deviation (Supplementary Figure 5)...")

message("\t* build models")

@@ -448,9 +449,9 @@ for(curVar in selColsNumeric) {



##### LINEAR MODELS FOR INFLUENCE OF DNA AMOUNT ON DEVIATION (Supplementary Figure 9) #####
##### LINEAR MODELS FOR INFLUENCE OF DNA AMOUNT ON DEVIATION (Supplementary Figure 10) #####

message("Influence of DNA amount on deviations (Supplementary Figure 9)...")
message("Influence of DNA amount on deviations (Supplementary Figure 10)...")

message("\t* build models")

@@ -666,7 +667,7 @@ lmResults$dataset <- factor(datasetTable[as.character(lmResults$dataset),"pretty

message("\t* make plots")

ggData <- lmResults[!is.na(lmResults$p),] #corResults[!is.na(corResults$p),] #
ggData <- lmResults[!is.na(lmResults$p),]

# the plot contains a few very extreme p-values, which need to be cropped so that the other p-values can be
# better visualized. We therefore trim mLog10P values above or below the 95% quantile (these values will be

0 comments on commit ff9ce55

Please sign in to comment.
You can’t perform that action at this time.