From 139ae79d9e6f21543f624bcac7db264eedc319bb Mon Sep 17 00:00:00 2001 From: corin Date: Mon, 2 May 2011 13:21:23 +0000 Subject: [PATCH] Update to work with latest eval format git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5719 348d0f76-0448-11de-a6fe-93d51630548a --- R/DataProcessingReport/GetTearsheetStats.R | 179 ++++++++++----------- 1 file changed, 86 insertions(+), 93 deletions(-) diff --git a/R/DataProcessingReport/GetTearsheetStats.R b/R/DataProcessingReport/GetTearsheetStats.R index e7ac7a436..119cf64bb 100644 --- a/R/DataProcessingReport/GetTearsheetStats.R +++ b/R/DataProcessingReport/GetTearsheetStats.R @@ -20,14 +20,13 @@ cmdargs = gsa.getargs( bamlist = list(value=NA, doc="list of BAM files"), evalroot = list(value=NA, doc="VariantEval root"), tearout = list(value=NA, doc="Output path for tearsheet PDF")#, - #plotout = list(value=NA, doc="Output path for PDF") + plotout = list(value=NA, doc="Output path for PDF") ), doc="Creates a tearsheet" ); bamlist = scan(cmdargs$bamlist, "character"); -#print(paste("grep SQUID ", cmdargs$yaml, ' |grep "C..." -o', sep="")) squids <- system(paste("grep SQUID ", cmdargs$yaml, ' |grep "C..." -o', sep=""), intern=TRUE) indexed = c(); nonindexed = c(); @@ -88,24 +87,22 @@ d2proj = d2[which(d2$"Project" %in% unique(dproj$Project) & d2$"Sample" %in% dpr tearsheet<-function(){ - tearsheetdrop <- "data/tearsheetdrop.jpg" #put the path to the tearsheet backdrop here + tearsheetdrop <- "~Documents/Sting/R/gsalib/data/tearsheetdrop.jpg" #put the path to the tearsheet backdrop here pdf(file= cmdargs$tearout, width=22, height=17, pagecentre=TRUE, pointsize=24) - + #define layout postable<-matrix(c(1, 1, 1, 1, 1, 1, rep(c(2, 2, 2, 4, 4, 4), 5), rep(c(3, 3, 3, 4, 4, 4), 3), rep(c(3,3,3,5,5,5), 5), 6,6,6,7,7,7), nrow=15, ncol=6, byrow=TRUE) layout(postable, heights=c(1, rep(.18, 13), 2), respect=FALSE) #prep for title bar - full<-strsplit(cmdargs$evalroot, "/") - name<-strsplit(full[[1]][length(full[[1]])], ".",fixed=TRUE)[[1]][1] drop<-read.jpeg(system.file(tearsheetdrop, package="gsalib")) #plot title bar par(mar=c(0,0,0,0)) plot(drop) - text(155, 50, name, family="serif", adj=c(0,0), cex=3, col=gray(.25)) + text(155, 50, "NonAutism_Walsh", family="serif", adj=c(0,0), cex=3, col=gray(.25)) # Project summary @@ -130,12 +127,12 @@ tearsheet<-function(){ callable_target = paste(na.omit(unique(dproj$"Target Territory")), collapse=", "); table1<-rbind(paste(used_samples," used samples/", unused_samples + used_samples," total samples", sep=""), sequencing_protocol, bait_design, callable_target) - print(nrow(table1)) rownames(table1)<-c("Samples","Sequencing Protocol", "Bait Design","Callable Target") par(mar=c(0,0,1,0)) textplot(table1, col.rownames="darkblue", show.colnames=FALSE, cex=1.25, valign="top") - title(main=sprintf("Project Summary (%s)\n", projects), family="sans", cex.main=1.25, line=-1) - # Bases summary + title(main=sprintf("Project Summary (%s)\n", projects), family="sans", cex.main=1.25, line=-1) + + # Bases summary reads_per_lane_mean = format(mean(dproj$"PF Reads (HS)", na.rm=TRUE), 8, 3,1, scientific=TRUE); reads_per_lane_sd = format(sd(dproj$"PF Reads (HS)", na.rm=TRUE), 8, 3,1, scientific=TRUE); @@ -188,7 +185,6 @@ tearsheet<-function(){ table2<-cbind(lanes, samps) colnames(table2)<-c("Per lane", "Per sample") - print(nrow(table2)) rownames(table2)<-c("Reads", "Used bases", "Average target coverage", "% loci covered to 10x", "% loci covered to 20x","% loci covered to 30x") par(mar=c(0,0,1,0)) @@ -246,31 +242,30 @@ tearsheet<-function(){ table3<-rbind(paste(instrument), used_lanes, sprintf("%s rejected by sequencing, %s by analysis\n", unused_lanes_by_sequencing, unused_lanes_by_analysis), sprintf("%0.1f +/- %0.1f lanes (median=%0.1f)\n", lanes_per_sample_mean, lanes_per_sample_sd, lanes_per_sample_median), sprintf("%s paired, %s widowed, %s single\n", lanes_paired, lanes_widowed, lanes_single), sprintf("%0.1f +/- %0.1f bases (median=%0.1f)\n", read_length_mean, read_length_sd, read_length_median), sprintf("\tSequencing dates: %s to %s\n", start_date, end_date)) - print(nrow(table3)) - print(table3) + rownames(table3)<-c("Sequencer", "Used lanes", "Unused lanes","Used lanes/sample", "Lane parities", "Read lengths", "Sequencing dates") par(mar=c(0,0,1,0)) textplot(table3, rmar=1, col.rownames="dark blue", show.colnames=FALSE, cex=1.25, valign="top") title(main="Sequencing Summary", family="sans", cex.main=1.25, line=0) +eval = gsa.read.gatkreport("NonAutism_Walsh.cleaned.snps_and_indels.filtered.annotated.eval") + + # Variant summary - - eval.counts = read.csv(paste(cmdargs$evalroot, ".Count_Variants.csv", sep=""), header=TRUE, comment.char="#"); - eval.counts.called = subset(eval.counts, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); - eval.counts.called.all = subset(eval.counts.called, novelty_name == "all")$nVariantLoci; - eval.counts.called.known = subset(eval.counts.called, novelty_name == "known")$nVariantLoci; - eval.counts.called.novel = subset(eval.counts.called, novelty_name == "novel")$nVariantLoci; +##TODO: Fix this csv reader + eval.counts = eval$CountVariants + eval.counts.all = subset(eval.counts, Novelty == "all")$nVariantLoci; + eval.counts.known = subset(eval.counts, Novelty == "known")$nVariantLoci; + eval.counts.novel = subset(eval.counts, Novelty == "novel")$nVariantLoci; - eval.titv = read.csv(paste(cmdargs$evalroot, ".Ti_slash_Tv_Variant_Evaluator.csv", sep=""), header=TRUE, comment.char="#"); - eval.titv.called = subset(eval.titv, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); - eval.titv.called.all = subset(eval.titv.called, novelty_name == "all")$ti.tv_ratio; - eval.titv.called.known = subset(eval.titv.called, novelty_name == "known")$ti.tv_ratio; - eval.titv.called.novel = subset(eval.titv.called, novelty_name == "novel")$ti.tv_ratio; + eval.titv = eval$TiTvVariantEvaluator + eval.titv.all = subset(eval.titv, Novelty == "all")$tiTvRatio; + eval.titv.known = subset(eval.titv, Novelty == "known")$tiTvRatio; + eval.titv.novel = subset(eval.titv, Novelty == "novel")$tiTvRatio; + + table4 = matrix(c(eval.counts.all, eval.counts.known, eval.counts.novel, eval.titv.all, eval.titv.known, eval.titv.novel, "3.0 - 3.2", "3.2 - 3.4", "2.7 - 3.0"), nrow=3); - table4 = matrix(c(eval.counts.called.all, eval.counts.called.known, eval.counts.called.novel, eval.titv.called.all, eval.titv.called.known, eval.titv.called.novel, "3.0 - 3.2", "3.2 - 3.4", "2.7 - 3.0"), nrow=3); - print(nrow(table4)) - print(paste("columns should be three, actually are:", ncol(table4))) rownames(table4) = c("All", "Known", "Novel"); colnames(table4) = c("Found", "Ti/Tv ratio", "Expected Ti/Tv ratio"); @@ -279,46 +274,45 @@ tearsheet<-function(){ par(mar=c(0,0,0,0)) textplot(table4, rmar=1, col.rownames="dark blue", cex=1.25, valign="top") title(main="Variant Summary", family="sans", cex.main=1.25, line=-2) - -#plots +# +# #plots +# #fix this reader +# eval.bysample = read.csv(paste(cmdargs$evalroot, ".SimpleMetricsBySample.csv", sep=""), header=TRUE, comment.char="#"); +# eval.bysample.called = subset(eval.bysample, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); +# eval.bysample.all = subset(eval.bysample.called, novelty_name == "all"); +# eval.bysample.known = subset(eval.bysample.called, novelty_name == "known"); +# eval.bysample.novel = subset(eval.bysample.called, novelty_name == "novel"); - eval.bysample = read.csv(paste(cmdargs$evalroot, ".SimpleMetricsBySample.csv", sep=""), header=TRUE, comment.char="#"); - eval.bysample.called = subset(eval.bysample, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); - eval.bysample.called.all = subset(eval.bysample.called, novelty_name == "all"); - eval.bysample.called.known = subset(eval.bysample.called, novelty_name == "known"); - eval.bysample.called.novel = subset(eval.bysample.called, novelty_name == "novel"); - - eval.ac = read.csv(paste(cmdargs$evalroot, ".MetricsByAc.csv", sep=""), header=TRUE, comment.char="#"); - eval.ac.called = subset(eval.ac, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); - eval.ac.called.all = subset(eval.ac.called, novelty_name == "all"); - eval.ac.called.known = subset(eval.ac.called, novelty_name == "known"); - eval.ac.called.novel = subset(eval.ac.called, novelty_name == "novel"); - - eval.func = read.csv(paste(cmdargs$evalroot, ".Functional_Class_Counts_by_Sample.csv", sep=""), header=TRUE, comment.char="#"); - eval.func.called = subset(eval.func, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); - eval.func.called.all = subset(eval.func.called, novelty_name == "all"); - eval.func.called.known = subset(eval.func.called, novelty_name == "known"); - eval.func.called.novel = subset(eval.func.called, novelty_name == "novel"); + eval.ac = eval$SimpleMetricsByAC.metrics + eval.ac.all = subset(eval.ac, Novelty == "all"); + eval.ac.known = subset(eval.ac, Novelty == "known"); + eval.ac.novel = subset(eval.ac, Novelty == "novel"); +# +# eval.func = read.csv(paste(cmdargs$evalroot, ".Functional_Class_Counts_by_Sample.csv", sep=""), header=TRUE, comment.char="#"); +# eval.func.called = subset(eval.func, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); +# eval.func.all = subset(eval.func.called, novelty_name == "all"); +# eval.func.known = subset(eval.func.called, novelty_name == "known"); +# eval.func.novel = subset(eval.func.called, novelty_name == "novel"); - #boxplot(eval.bysample.called.all$CountVariants, eval.bysample.called.known$CountVariants, eval.bysample.called.novel$CountVariants, names=c("All", "Known", "Novel"), ylab="Variants per sample", main="", cex=1.3, cex.lab=1.3, cex.axis=1.3); + #boxplot(eval.bysample.all$CountVariants, eval.bysample.known$CountVariants, eval.bysample.novel$CountVariants, names=c("All", "Known", "Novel"), ylab="Variants per sample", main="", cex=1.3, cex.lab=1.3, cex.axis=1.3); + +# par(mar=c(5, 4, 4, 2) + 0.1) +# ind = order(eval.bysample.all$CountVariants); +# plot(c(1:length(eval.bysample.all$CountVariants)), eval.bysample.all$CountVariants[ind], col="black", cex=1.1, cex.lab=1.1, cex.axis=1.1, main="Variants per Sample", xlab="Sample", ylab="Number of variants", bty="n", ylim=c(0, max(eval.bysample.all$CountVariants))); +# points(c(1:length(eval.bysample.known$CountVariants)), eval.bysample.known$CountVariants[ind], col="blue", cex=1.3); +# points(c(1:length(eval.bysample.novel$CountVariants)), eval.bysample.novel$CountVariants[ind], col="red", cex=1.3); +# legend("right", max(eval.bysample.all$CountVariants)/2, c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); par(mar=c(5, 4, 4, 2) + 0.1) - ind = order(eval.bysample.called.all$CountVariants); - plot(c(1:length(eval.bysample.called.all$CountVariants)), eval.bysample.called.all$CountVariants[ind], col="black", cex=1.1, cex.lab=1.1, cex.axis=1.1, main="Variants per Sample", xlab="Sample", ylab="Number of variants", bty="n", ylim=c(0, max(eval.bysample.called.all$CountVariants))); - points(c(1:length(eval.bysample.called.known$CountVariants)), eval.bysample.called.known$CountVariants[ind], col="blue", cex=1.3); - points(c(1:length(eval.bysample.called.novel$CountVariants)), eval.bysample.called.novel$CountVariants[ind], col="red", cex=1.3); - legend("right", max(eval.bysample.called.all$CountVariants)/2, c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); - - par(mar=c(5, 4, 4, 2) + 0.1) - plot(eval.ac.called.all$AC, eval.ac.called.all$n, col="black", type="l", lwd=2, cex=1.1, cex.lab=1.1, cex.axis=1.1, xlab="Allele count", ylab="Number of variants", main="Variants by Allele Count", log="xy", bty="n"); - points(eval.ac.called.known$AC, eval.ac.called.known$n, col="blue", type="l", lwd=2); - points(eval.ac.called.novel$AC, eval.ac.called.novel$n, col="red", type="l", lwd=2); + plot(eval.ac.all$AC, eval.ac.all$n, col="black", type="l", lwd=2, cex=1.1, cex.lab=1.1, cex.axis=1.1, xlab="Allele count", ylab="Number of variants", main="Variants by Allele Count", log="xy", bty="n"); + points(eval.ac.known$AC, eval.ac.known$n, col="blue", type="l", lwd=2); + points(eval.ac.novel$AC, eval.ac.novel$n, col="red", type="l", lwd=2); legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), lwd=2); - #plot(eval.func.called.all$Synonymous[ind] / (eval.func.called.all$Missense + eval.func.called.all$Nonsense)[ind], ylim=c(0, 2), cex=1.3, cex.lab=1.3, cex.axis=1.3, bty="n", xlab="Sample", ylab="Ratio of synonymous to non-synonymous variants", col="black"); - #points(eval.func.called.known$Synonymous[ind] / (eval.func.called.known$Missense + eval.func.called.known$Nonsense)[ind], cex=1.3, col="blue"); - #points(eval.func.called.novel$Synonymous[ind] / (eval.func.called.novel$Missense + eval.func.called.novel$Nonsense)[ind], cex=1.3, col="red"); + #plot(eval.func.all$Synonymous[ind] / (eval.func.all$Missense + eval.func.all$Nonsense)[ind], ylim=c(0, 2), cex=1.3, cex.lab=1.3, cex.axis=1.3, bty="n", xlab="Sample", ylab="Ratio of synonymous to non-synonymous variants", col="black"); + #points(eval.func.known$Synonymous[ind] / (eval.func.known$Missense + eval.func.known$Nonsense)[ind], cex=1.3, col="blue"); + #points(eval.func.novel$Synonymous[ind] / (eval.func.novel$Missense + eval.func.novel$Nonsense)[ind], cex=1.3, col="red"); #legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); @@ -330,44 +324,43 @@ tearsheet() # Plots plots<-function(){ -eval.bysample = read.csv(paste(cmdargs$evalroot, ".SimpleMetricsBySample.csv", sep=""), header=TRUE, comment.char="#"); -eval.bysample.called = subset(eval.bysample, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); -eval.bysample.called.all = subset(eval.bysample.called, novelty_name == "all"); -eval.bysample.called.known = subset(eval.bysample.called, novelty_name == "known"); -eval.bysample.called.novel = subset(eval.bysample.called, novelty_name == "novel"); +# eval.bysample = read.csv(paste(cmdargs$evalroot, ".SimpleMetricsBySample.csv", sep=""), header=TRUE, comment.char="#"); +# eval.bysample.called = subset(eval.bysample, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); +# eval.bysample.all = subset(eval.bysample.called, novelty_name == "all"); +# eval.bysample.known = subset(eval.bysample.called, novelty_name == "known"); +# eval.bysample.novel = subset(eval.bysample.called, novelty_name == "novel"); -eval.ac = read.csv(paste(cmdargs$evalroot, ".MetricsByAc.csv", sep=""), header=TRUE, comment.char="#"); -eval.ac.called = subset(eval.ac, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); -eval.ac.called.all = subset(eval.ac.called, novelty_name == "all"); -eval.ac.called.known = subset(eval.ac.called, novelty_name == "known"); -eval.ac.called.novel = subset(eval.ac.called, novelty_name == "novel"); -eval.func = read.csv(paste(cmdargs$evalroot, ".Functional_Class_Counts_by_Sample.csv", sep=""), header=TRUE, comment.char="#"); -eval.func.called = subset(eval.func, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); -eval.func.called.all = subset(eval.func.called, novelty_name == "all"); -eval.func.called.known = subset(eval.func.called, novelty_name == "known"); -eval.func.called.novel = subset(eval.func.called, novelty_name == "novel"); + eval.ac = eval$SimpleMetricsByAC.metrics + eval.ac.all = subset(eval.ac.called, Novelty == "all"); + eval.ac.known = subset(eval.ac.called, Novelty == "known"); + eval.ac.novel = subset(eval.ac.called, Novelty == "novel"); +# +# eval.func = read.csv(paste(cmdargs$evalroot, ".Functional_Class_Counts_by_Sample.csv", sep=""), header=TRUE, comment.char="#"); +# eval.func.called = subset(eval.func, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); +# eval.func.all = subset(eval.func.called, novelty_name == "all"); +# eval.func.known = subset(eval.func.called, novelty_name == "known"); +# eval.func.novel = subset(eval.func.called, novelty_name == "novel"); -pdf(cmdargs$plotout); + pdf(file= cmdargs$plotout, width=22, height=17, pagecentre=TRUE, pointsize=24) +# +# boxplot(eval.bysample.all$CountVariants, eval.bysample.known$CountVariants, eval.bysample.novel$CountVariants, names=c("All", "Known", "Novel"), ylab="Variants per sample", main="", cex=1.3, cex.lab=1.3, cex.axis=1.3); +# +# ind = order(eval.bysample.all$CountVariants); +# plot(c(1:length(eval.bysample.all$CountVariants)), eval.bysample.all$CountVariants[ind], col="black", cex=1.3, cex.lab=1.3, cex.axis=1.3, xlab="Sample", ylab="Number of variants", bty="n", ylim=c(0, max(eval.bysample.all$CountVariants))); +# points(c(1:length(eval.bysample.known$CountVariants)), eval.bysample.known$CountVariants[ind], col="blue", cex=1.3); +# points(c(1:length(eval.bysample.novel$CountVariants)), eval.bysample.novel$CountVariants[ind], col="red", cex=1.3); +# legend(0, max(eval.bysample.all$CountVariants)/2, c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); -boxplot(eval.bysample.called.all$CountVariants, eval.bysample.called.known$CountVariants, eval.bysample.called.novel$CountVariants, names=c("All", "Known", "Novel"), ylab="Variants per sample", main="", cex=1.3, cex.lab=1.3, cex.axis=1.3); - -ind = order(eval.bysample.called.all$CountVariants); -plot(c(1:length(eval.bysample.called.all$CountVariants)), eval.bysample.called.all$CountVariants[ind], col="black", cex=1.3, cex.lab=1.3, cex.axis=1.3, xlab="Sample", ylab="Number of variants", bty="n", ylim=c(0, max(eval.bysample.called.all$CountVariants))); -points(c(1:length(eval.bysample.called.known$CountVariants)), eval.bysample.called.known$CountVariants[ind], col="blue", cex=1.3); -points(c(1:length(eval.bysample.called.novel$CountVariants)), eval.bysample.called.novel$CountVariants[ind], col="red", cex=1.3); -legend(0, max(eval.bysample.called.all$CountVariants)/2, c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); - -plot(eval.ac.called.all$AC, eval.ac.called.all$n, col="black", type="l", lwd=2, cex=1.3, cex.lab=1.3, cex.axis=1.3, xlab="Allele count", ylab="Number of variants", main="", log="xy", bty="n"); -points(eval.ac.called.known$AC, eval.ac.called.known$n, col="blue", type="l", lwd=2); -points(eval.ac.called.novel$AC, eval.ac.called.novel$n, col="red", type="l", lwd=2); +plot(eval.ac.all$AC, eval.ac.all$n, col="black", type="l", lwd=2, cex=1.3, cex.lab=1.3, cex.axis=1.3, xlab="Allele count", ylab="Number of variants", main="", log="xy", bty="n"); +points(eval.ac.known$AC, eval.ac.known$n, col="blue", type="l", lwd=2); +points(eval.ac.novel$AC, eval.ac.novel$n, col="red", type="l", lwd=2); legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), lwd=2); - -plot(eval.func.called.all$Synonymous[ind] / (eval.func.called.all$Missense + eval.func.called.all$Nonsense)[ind], ylim=c(0, 2), cex=1.3, cex.lab=1.3, cex.axis=1.3, bty="n", xlab="Sample", ylab="Ratio of synonymous to non-synonymous variants", col="black"); -points(eval.func.called.known$Synonymous[ind] / (eval.func.called.known$Missense + eval.func.called.known$Nonsense)[ind], cex=1.3, col="blue"); -points(eval.func.called.novel$Synonymous[ind] / (eval.func.called.novel$Missense + eval.func.called.novel$Nonsense)[ind], cex=1.3, col="red"); -legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); +# +# plot(eval.func.all$Synonymous[ind] / (eval.func.all$Missense + eval.func.all$Nonsense)[ind], ylim=c(0, 2), cex=1.3, cex.lab=1.3, cex.axis=1.3, bty="n", xlab="Sample", ylab="Ratio of synonymous to non-synonymous variants", col="black"); +# points(eval.func.known$Synonymous[ind] / (eval.func.known$Missense + eval.func.known$Nonsense)[ind], cex=1.3, col="blue"); +# points(eval.func.novel$Synonymous[ind] / (eval.func.novel$Missense + eval.func.novel$Nonsense)[ind], cex=1.3, col="red"); +# legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21); dev.off(); } -