zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Default R script now plots sensitivity/specificity curve 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4825 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
depristo 2010-12-13 16:55:11 +00:00
parent 7bc2049031
commit a6397ed8c3
1 changed files with 17 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
R/plot_Tranches.R
View File

@ -5,7 +5,8 @@ verbose = TRUE
tranchesFile = args[1] tranchesFile = args[1]
targetTITV = as.numeric(args[2]) targetTITV = as.numeric(args[2])
suppressLegend = ! is.na(args[3]) targetSensitivity = as.numeric(args[3])
suppressLegend = ! is.na(args[4])
# ----------------------------------------------------------------------------------------------- # -----------------------------------------------------------------------------------------------
# Useful general routines # Useful general routines
@ -38,7 +39,8 @@ leftShift <- function(x, leftValue = 0) {
# Tranches plot # Tranches plot
# ----------------------------------------------------------------------------------------------- # -----------------------------------------------------------------------------------------------
data2 = read.table(tranchesFile,sep=",",head=T) data2 = read.table(tranchesFile,sep=",",head=T)
data2 = data2[order(data2$FDRtranche, decreasing=T),] data2 = data2[order(data2$novelTiTv, decreasing=F),]
#data2 = data2[order(data2$FDRtranche, decreasing=T),]
cols = c("cornflowerblue", "cornflowerblue", "darkorange", "darkorange") cols = c("cornflowerblue", "cornflowerblue", "darkorange", "darkorange")
density=c(20, -1, -1, 20) density=c(20, -1, -1, 20)
outfile = paste(tranchesFile, ".pdf", sep="") outfile = paste(tranchesFile, ".pdf", sep="")
@ -64,10 +66,22 @@ barplot(d/1000,horiz=TRUE,col=cols,space=0.2,xlab="Number of Novel Variants (100
#abline(v= d[2,dim(d)[2]], lty=2) #abline(v= d[2,dim(d)[2]], lty=2)
#abline(v= d[1,3], lty=2) #abline(v= d[1,3], lty=2)
if ( ! suppressLegend ) if ( ! suppressLegend )
legend("topright", c('Cumulative TPs','Tranch-specific TPs', 'Tranch-specific FPs', 'Cumulative FPs' ), fill=cols, density=density, bg='white', cex=1.25) legend(5, 1.75, c('Cumulative TPs','Tranch-specific TPs', 'Tranch-specific FPs', 'Cumulative FPs' ), fill=cols, density=density, bg='white', cex=1.25)
mtext("Ti/Tv",2,line=2.25,at=length(data2$FDRtranche)*1.2,las=1, cex=1) mtext("Ti/Tv",2,line=2.25,at=length(data2$FDRtranche)*1.2,las=1, cex=1)
mtext("FDR",2,line=0,at=length(data2$FDRtranche)*1.2,las=1, cex=1) mtext("FDR",2,line=0,at=length(data2$FDRtranche)*1.2,las=1, cex=1)
axis(2,line=-1,at=0.7+(0:(length(data2$FDRtranche)-1))*1.2,tick=FALSE,labels=data2$FDRtranche, las=1, cex.axis=1.0) axis(2,line=-1,at=0.7+(0:(length(data2$FDRtranche)-1))*1.2,tick=FALSE,labels=data2$FDRtranche, las=1, cex.axis=1.0)
axis(2,line=1,at=0.7+(0:(length(data2$FDRtranche)-1))*1.2,tick=FALSE,labels=round(novelTiTv,3), las=1, cex.axis=1.0) axis(2,line=1,at=0.7+(0:(length(data2$FDRtranche)-1))*1.2,tick=FALSE,labels=round(novelTiTv,3), las=1, cex.axis=1.0)
# plot sensitivity vs. specificity
sensitivity = data2$truthSensitivity
if ( ! is.null(sensitivity) ) {
#specificity = titvFPEstV(targetTITV, novelTiTv)
specificity = novelTiTv
plot(sensitivity, specificity, type="b", col="cornflowerblue", xlab="Tranche truth sensitivity", ylab="Specificity (Novel Ti/Tv ratio)")
abline(h=targetTITV, lty=2)
abline(v=targetSensitivity, lty=2)
#text(max(sensitivity), targetTITV-0.05, labels="Expected novel Ti/Tv", pos=2)
}
dev.off() dev.off()