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

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
This commit is contained in:
corin 2011-05-02 13:21:23 +00:00
parent 2d81262f87
commit 139ae79d9e
1 changed files with 86 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

179
R/DataProcessingReport/GetTearsheetStats.R
View File

@ -20,14 +20,13 @@ cmdargs = gsa.getargs(
bamlist = list(value=NA, doc="list of BAM files"), bamlist = list(value=NA, doc="list of BAM files"),
evalroot = list(value=NA, doc="VariantEval root"), evalroot = list(value=NA, doc="VariantEval root"),
tearout = list(value=NA, doc="Output path for tearsheet PDF")#, 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" doc="Creates a tearsheet"
); );
bamlist = scan(cmdargs$bamlist, "character"); 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) squids <- system(paste("grep SQUID ", cmdargs$yaml, ' |grep "C..." -o', sep=""), intern=TRUE)
indexed = c(); indexed = c();
nonindexed = c(); nonindexed = c();
@ -88,24 +87,22 @@ d2proj = d2[which(d2$"Project" %in% unique(dproj$Project) & d2$"Sample" %in% dpr
tearsheet<-function(){ 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) pdf(file= cmdargs$tearout, width=22, height=17, pagecentre=TRUE, pointsize=24)
#define layout #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) 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) layout(postable, heights=c(1, rep(.18, 13), 2), respect=FALSE)
#prep for title bar #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")) drop<-read.jpeg(system.file(tearsheetdrop, package="gsalib"))
#plot title bar #plot title bar
par(mar=c(0,0,0,0)) par(mar=c(0,0,0,0))
plot(drop) 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 # Project summary
@ -130,12 +127,12 @@ tearsheet<-function(){
callable_target = paste(na.omit(unique(dproj$"Target Territory")), collapse=", "); 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) 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") rownames(table1)<-c("Samples","Sequencing Protocol", "Bait Design","Callable Target")
par(mar=c(0,0,1,0)) par(mar=c(0,0,1,0))
textplot(table1, col.rownames="darkblue", show.colnames=FALSE, cex=1.25, valign="top") 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) title(main=sprintf("Project Summary (%s)\n", projects), family="sans", cex.main=1.25, line=-1)
# Bases summary
# Bases summary
reads_per_lane_mean = format(mean(dproj$"PF Reads (HS)", na.rm=TRUE), 8, 3,1, scientific=TRUE); 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); 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) table2<-cbind(lanes, samps)
colnames(table2)<-c("Per lane", "Per sample") 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") 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)) 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)) 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") rownames(table3)<-c("Sequencer", "Used lanes", "Unused lanes","Used lanes/sample", "Lane parities", "Read lengths", "Sequencing dates")
par(mar=c(0,0,1,0)) par(mar=c(0,0,1,0))
textplot(table3, rmar=1, col.rownames="dark blue", show.colnames=FALSE, cex=1.25, valign="top") 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) 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 # Variant summary
##TODO: Fix this csv reader
eval.counts = read.csv(paste(cmdargs$evalroot, ".Count_Variants.csv", sep=""), header=TRUE, comment.char="#"); eval.counts = eval$CountVariants
eval.counts.called = subset(eval.counts, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); eval.counts.all = subset(eval.counts, Novelty == "all")$nVariantLoci;
eval.counts.called.all = subset(eval.counts.called, novelty_name == "all")$nVariantLoci; eval.counts.known = subset(eval.counts, Novelty == "known")$nVariantLoci;
eval.counts.called.known = subset(eval.counts.called, novelty_name == "known")$nVariantLoci; eval.counts.novel = subset(eval.counts, Novelty == "novel")$nVariantLoci;
eval.counts.called.novel = subset(eval.counts.called, novelty_name == "novel")$nVariantLoci;
eval.titv = read.csv(paste(cmdargs$evalroot, ".Ti_slash_Tv_Variant_Evaluator.csv", sep=""), header=TRUE, comment.char="#"); eval.titv = eval$TiTvVariantEvaluator
eval.titv.called = subset(eval.titv, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); eval.titv.all = subset(eval.titv, Novelty == "all")$tiTvRatio;
eval.titv.called.all = subset(eval.titv.called, novelty_name == "all")$ti.tv_ratio; eval.titv.known = subset(eval.titv, Novelty == "known")$tiTvRatio;
eval.titv.called.known = subset(eval.titv.called, novelty_name == "known")$ti.tv_ratio; eval.titv.novel = subset(eval.titv, Novelty == "novel")$tiTvRatio;
eval.titv.called.novel = subset(eval.titv.called, novelty_name == "novel")$ti.tv_ratio;
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"); rownames(table4) = c("All", "Known", "Novel");
colnames(table4) = c("Found", "Ti/Tv ratio", "Expected Ti/Tv ratio"); colnames(table4) = c("Found", "Ti/Tv ratio", "Expected Ti/Tv ratio");
@ -279,46 +274,45 @@ tearsheet<-function(){
par(mar=c(0,0,0,0)) par(mar=c(0,0,0,0))
textplot(table4, rmar=1, col.rownames="dark blue", cex=1.25, valign="top") 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) 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.ac = eval$SimpleMetricsByAC.metrics
eval.bysample.called = subset(eval.bysample, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); eval.ac.all = subset(eval.ac, Novelty == "all");
eval.bysample.called.all = subset(eval.bysample.called, novelty_name == "all"); eval.ac.known = subset(eval.ac, Novelty == "known");
eval.bysample.called.known = subset(eval.bysample.called, novelty_name == "known"); eval.ac.novel = subset(eval.ac, Novelty == "novel");
eval.bysample.called.novel = subset(eval.bysample.called, novelty_name == "novel"); #
# eval.func = read.csv(paste(cmdargs$evalroot, ".Functional_Class_Counts_by_Sample.csv", sep=""), header=TRUE, comment.char="#");
eval.ac = read.csv(paste(cmdargs$evalroot, ".MetricsByAc.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.ac.called = subset(eval.ac, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); # eval.func.all = subset(eval.func.called, novelty_name == "all");
eval.ac.called.all = subset(eval.ac.called, novelty_name == "all"); # eval.func.known = subset(eval.func.called, novelty_name == "known");
eval.ac.called.known = subset(eval.ac.called, novelty_name == "known"); # eval.func.novel = subset(eval.func.called, novelty_name == "novel");
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");
#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) par(mar=c(5, 4, 4, 2) + 0.1)
ind = order(eval.bysample.called.all$CountVariants); 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");
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(eval.ac.known$AC, eval.ac.known$n, col="blue", type="l", lwd=2);
points(c(1:length(eval.bysample.called.known$CountVariants)), eval.bysample.called.known$CountVariants[ind], col="blue", cex=1.3); points(eval.ac.novel$AC, eval.ac.novel$n, col="red", type="l", lwd=2);
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);
legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), 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"); #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.called.known$Synonymous[ind] / (eval.func.called.known$Missense + eval.func.called.known$Nonsense)[ind], cex=1.3, col="blue"); #points(eval.func.known$Synonymous[ind] / (eval.func.known$Missense + eval.func.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"); #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); #legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21);
@ -330,44 +324,43 @@ tearsheet()
# Plots # Plots
plots<-function(){ plots<-function(){
eval.bysample = read.csv(paste(cmdargs$evalroot, ".SimpleMetricsBySample.csv", sep=""), header=TRUE, comment.char="#"); # 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 = 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.all = subset(eval.bysample.called, novelty_name == "all");
eval.bysample.called.known = subset(eval.bysample.called, novelty_name == "known"); # eval.bysample.known = subset(eval.bysample.called, novelty_name == "known");
eval.bysample.called.novel = subset(eval.bysample.called, novelty_name == "novel"); # 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.ac = eval$SimpleMetricsByAC.metrics
eval.func.called = subset(eval.func, evaluation_name == "eval" & comparison_name == "dbsnp" & jexl_expression == "none" & filter_name == "called"); eval.ac.all = subset(eval.ac.called, Novelty == "all");
eval.func.called.all = subset(eval.func.called, novelty_name == "all"); eval.ac.known = subset(eval.ac.called, Novelty == "known");
eval.func.called.known = subset(eval.func.called, novelty_name == "known"); eval.ac.novel = subset(eval.ac.called, Novelty == "novel");
eval.func.called.novel = subset(eval.func.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.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); 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);
ind = order(eval.bysample.called.all$CountVariants); points(eval.ac.novel$AC, eval.ac.novel$n, col="red", type="l", lwd=2);
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);
legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), 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"); # 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.called.known$Synonymous[ind] / (eval.func.called.known$Missense + eval.func.called.known$Nonsense)[ind], cex=1.3, col="blue"); # points(eval.func.known$Synonymous[ind] / (eval.func.known$Missense + eval.func.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"); # 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); # legend("topright", c("All", "Known", "Novel"), col=c("black", "blue", "red"), pt.cex=1.3, pch=21);
dev.off(); dev.off();
} }