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

Fixed a bug where the script looked for the wrong column name. Also, all results are now returned in a single plot. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4216 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2010-09-06 14:19:57 +00:00
parent 4d4ef5b42c
commit 19e22cfa87
2 changed files with 16 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
R/VariantReport/VariantReport.R
View File

@ -13,51 +13,27 @@ if (interactive()) {
eval = read.eval(evalRoot);
plot.begin(plotRoot, "variantReport");
# Venn diagram
plot.begin(plotRoot, "venn");
plot.callsetConcordance(eval);
plot.end(plotRoot);
# Venn by AC
plot.begin(plotRoot, "venn_by_ac.all", width=12, height=8);
plot.callsetConcordanceByAC(eval, novelty_name="all");
plot.end(plotRoot);
plot.begin(plotRoot, "venn_by_ac.known", width=12, height=8);
plot.callsetConcordanceByAC(eval, novelty_name="known");
plot.end(plotRoot);
plot.begin(plotRoot, "venn_by_ac.novel", width=12, height=8);
plot.callsetConcordanceByAC(eval, novelty_name="novel");
plot.end(plotRoot);
# Allele count spectrum
plot.begin(plotRoot, "acs.all", width=12, height=8);
plot.alleleCountSpectrum(eval, novelty_name="all");
plot.end(plotRoot);
plot.begin(plotRoot, "acs.known", width=12, height=8);
plot.alleleCountSpectrum(eval, novelty_name="known");
plot.end(plotRoot);
plot.begin(plotRoot, "acs.novel", width=12, height=8);
plot.alleleCountSpectrum(eval, novelty_name="novel");
plot.end(plotRoot);
# Ti/Tv spectrum
plot.begin(plotRoot, "titv.all", width=12, height=8);
plot.titvSpectrum(eval, novelty_name="all");
plot.end(plotRoot);
plot.begin(plotRoot, "titv.known", width=12, height=8);
plot.titvSpectrum(eval, novelty_name="known");
plot.end(plotRoot);
plot.begin(plotRoot, "titv.novel", width=12, height=8);
plot.titvSpectrum(eval, novelty_name="novel");
plot.end(plotRoot);
# Per-sample
#plot.begin(plotRoot, "variants_per_sample", width=12, height=8);
#plot.variantsPerSample(eval);
#plot.end(plotRoot);
plot.end(plotRoot);

31
R/gsacommons.R
View File

@ -134,16 +134,16 @@ eval.getMetrics <- function(eval, jexl_expression) {
callset.counts.titv = eval$TiTv[which(eval$TiTv$evaluation_name == "eval" & eval$TiTv$comparison_name == "dbsnp" & eval$TiTv$jexl_expression == jexl_expression),];
callset.calledCounts = callset.counts[which(callset.counts$filter_name == "called" & callset.counts$novelty_name == "all"),]$nVariantLoci;
callset.calledCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "all"),]$ti.tv.ratio;
callset.calledCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "all"),]$ti.tv_ratio;
callset.knownCounts = callset.counts[which(callset.counts$filter_name == "called" & callset.counts$novelty_name == "known"),]$nVariantLoci;
callset.knownCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "known"),]$ti.tv.ratio;
callset.knownCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "known"),]$ti.tv_ratio;
callset.novelCounts = callset.counts[which(callset.counts$filter_name == "called" & callset.counts$novelty_name == "novel"),]$nVariantLoci;
callset.novelCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "novel"),]$ti.tv.ratio;
callset.novelCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "called" & callset.counts.titv$novelty_name == "novel"),]$ti.tv_ratio;
callset.filteredCounts = callset.counts[which(callset.counts$filter_name == "filtered" & callset.counts$novelty_name == "all"),]$nVariantLoci;
callset.filteredCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "filtered" & callset.counts.titv$novelty_name == "all"),]$ti.tv.ratio;
callset.filteredCounts.titv = callset.counts.titv[which(callset.counts.titv$filter_name == "filtered" & callset.counts.titv$novelty_name == "all"),]$ti.tv_ratio;
metrics = list(
all = callset.calledCounts,
@ -181,8 +181,8 @@ plot.callsetConcordance <- function(eval, col=c("#FF6342", "#63C6DE", "#ADDE63")
filtered = intersection$filtered + aonly$filtered + bonly$filtered
);
par.def = par(no.readonly = TRUE);
par(mar=c(0.1, 0.1, 0.1, 0.1));
#par.def = par(no.readonly = TRUE);
#par(mar=c(0.1, 0.1, 0.1, 0.1));
plot.venn(aonly$all + intersection$all, bonly$all + intersection$all, 0, intersection$all, 0, 0, pos=c(0.32, 0.32, 0.68, 0.70), col=col);
@ -190,14 +190,10 @@ plot.callsetConcordance <- function(eval, col=c("#FF6342", "#63C6DE", "#ADDE63")
text(0.5, 0.75, cex=1.2, adj=c(0.5, 0.33), .plot.callsetConcordance.getLabelText("Intersection", intersection, union));
text(1, 0.5, cex=1.2, pos=2, .plot.callsetConcordance.getLabelText(eval$CallsetNames[2], bonly, union));
par(par.def);
#par(par.def);
}
plot.callsetConcordanceByAC <- function(eval, normalize=TRUE, novelty_name="all", col=c("#FF6342", "#5C92A4", "#55BBFF")) {
#aonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[1]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#intersection = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & eval$MetricsByAc$jexl_expression == "Intersection" & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#bonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[2]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
aonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[1], novelty_name);
intersection = eval.getMetricsByAc(eval, "Intersection", novelty_name);
bonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[2], novelty_name);
@ -214,6 +210,9 @@ plot.callsetConcordanceByAC <- function(eval, normalize=TRUE, novelty_name="all"
bonly$n = 0;
}
#par.def = par(no.readonly = TRUE);
#par(mar=c(5, 5, 3, 5));
if (normalize == TRUE) {
norm = aonly$n + intersection$n + bonly$n;
matnorm = rbind(aonly$n/norm, intersection$n/norm, bonly$n/norm);
@ -226,13 +225,11 @@ plot.callsetConcordanceByAC <- function(eval, normalize=TRUE, novelty_name="all"
}
legend("topright", c(eval$CallsetOnlyNames[1], "Intersection", eval$CallsetOnlyNames[2]), fill=col);
#par(par.def);
}
plot.alleleCountSpectrum <- function(eval, novelty_name="all", col=c("#FF6342", "#5C92A4", "#55BBFF")) {
#aonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[1]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#intersection = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & eval$MetricsByAc$jexl_expression == "Intersection" & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#bonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[2]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
aonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[1], novelty_name);
intersection = eval.getMetricsByAc(eval, "Intersection", novelty_name);
bonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[2], novelty_name);
@ -262,10 +259,6 @@ eval.getMetricsByAc <- function(eval, jexl_expression, novelty_name="all") {
}
plot.titvSpectrum <- function(eval, novelty_name="all", col=c("#FF6342", "#5C92A4", "#55BBFF")) {
#aonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[1]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#intersection = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & eval$MetricsByAc$jexl_expression == "Intersection" & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
#bonly = eval$MetricsByAc[which(eval$MetricsByAc$evaluation_name == "eval" & eval$MetricsByAc$comparison_name == "dbsnp" & as.character(eval$MetricsByAc$jexl_expression) == as.character(eval$CallsetOnlyNames[2]) & eval$MetricsByAc$filter_name == "called" & eval$MetricsByAc$novelty_name == novelty_name),];
aonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[1], novelty_name);
intersection = eval.getMetricsByAc(eval, "Intersection", novelty_name);
bonly = eval.getMetricsByAc(eval, eval$CallsetOnlyNames[2], novelty_name);