From 19ee6116fc2f04d1e30e7025e340ff6bed0f646c Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Tue, 9 Oct 2012 17:20:04 +0200 Subject: [PATCH 1/7] Minor bugfix of the name of the plot function and a hardcoded ylim within the plot fcn --- bam/bedtoolshist2pdf.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/bam/bedtoolshist2pdf.R b/bam/bedtoolshist2pdf.R index 9ec9e85..3a6375a 100644 --- a/bam/bedtoolshist2pdf.R +++ b/bam/bedtoolshist2pdf.R @@ -42,8 +42,8 @@ main <- function(histfile2id, tab.outfile, hist.pdf.file){ hist.tabs = read.files(histfile2id) dp2frac = get.dp2frac(hist.tabs) dp2bases = get.dp2bases(hist.tabs) - plot.dp2frac(dp2frac, ltypes = histfile2id[names(dp2frac), 'ltypes'], colors = histfile2id[names(dp2frac), 'colors']) - plot.dp2bases(dp2bases) + plot.dp2y(dp2frac, ltypes = histfile2id[names(dp2frac), 'ltypes'], colors = histfile2id[names(dp2frac), 'colors']) + plot.dp2y(dp2bases) #Plot @@ -66,11 +66,11 @@ main <- function(histfile2id, tab.outfile, hist.pdf.file){ par(mfrow = c(1, 2)) unstim.ind = grep('unstim', names(dp2frac)) dp2frac.unstim = dp2frac[unstim.ind] - plot.dp2frac(dp2frac.unstim, ltypes = histfile2id[names(dp2frac.unstim), 'ltypes'], colors = histfile2id[names(dp2frac.unstim), 'colors']) + plot.dp2y(dp2frac.unstim, ltypes = histfile2id[names(dp2frac.unstim), 'ltypes'], colors = histfile2id[names(dp2frac.unstim), 'colors']) lps.ind = grep('LPS', names(dp2frac)) dp2frac.lps = dp2frac[lps.ind] - plot.dp2frac(dp2frac.lps, ltypes = histfile2id[names(dp2frac.lps), 'ltypes'], colors = histfile2id[names(dp2frac.lps), 'colors']) + plot.dp2y(dp2frac.lps, ltypes = histfile2id[names(dp2frac.lps), 'ltypes'], colors = histfile2id[names(dp2frac.lps), 'colors']) dev.off() } @@ -145,7 +145,7 @@ plot.dp2y <- function(dp2frac, max.plot.dp = 100, ylab = 'fraction covered', lty n.hist = length(dp2frac) ylim = range(unlist(dp2frac)) - ylim = c(0, 5e8) + #ylim = c(0, 5e8) jhist = 1 plot(dp2frac[[jhist]][min.plot.dp:max.plot.dp], xlab = 'depth', ylab = ylab, col=colors[jhist], type = 'l', ylim = ylim, lty = ltypes[jhist]) for(jhist in 2:n.hist){ From 66247bd603a725e5591865bca5c59087b7bd9ab7 Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Wed, 10 Oct 2012 14:21:13 +0200 Subject: [PATCH 2/7] changed the hardcoded paths in snv.ase.R --- ase/snv.ase.R | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/ase/snv.ase.R b/ase/snv.ase.R index 5132ffb..a9f6326 100644 --- a/ase/snv.ase.R +++ b/ase/snv.ase.R @@ -4,12 +4,13 @@ sys = 'local' if(sys == 'kalkyl'){ } -if(sys == 'local'){ - varcalls.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/varcalls' - basecount.datadir='~/Dropbox/postdoc/projects/ase/data/sim/basecount' - ase.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/ase' - resdir = '~/Dropbox/postdoc/projects/ase/res/sim' - annovar.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/annovar' +if(sys == 'local'){ + basedir = '~/Dropbox/postdoc/projects/ase_msc_benchmarking/data/gsnap' + varcalls.datadir = file.path(basedir, 'varcalls') + basecount.datadir = file.path(basedir, 'basecount') + annovar.datadir = file.path(basedir, 'annovar') + ase.datadir = file.path(basedir, 'ase') + resdir = '~/Dropbox/postdoc/projects/ase_msc_benchmarking/res/gsnap/ase' } #mpileup calls @@ -66,7 +67,7 @@ main <- function(){ bc.filt = lapply(bc.filt, function(bc){bc[which(bc[, 'site.dp'] >= min.dp), ]}) #total number of vars - length(unique(unlist(lapply(bc.filt, '[[', 'site')))) #112756 + length(unique(unlist(lapply(bc.filt, '[[', 'site')))) #97628 ### #DP dist From 028f6392b3d0d1fd66ef78380c73ea13b7dafa14 Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Thu, 11 Oct 2012 14:09:01 +0200 Subject: [PATCH 3/7] Added dbSNP filter step of called variants. --- ase/snv.ase.R | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/ase/snv.ase.R b/ase/snv.ase.R index a9f6326..9063121 100644 --- a/ase/snv.ase.R +++ b/ase/snv.ase.R @@ -5,12 +5,13 @@ sys = 'local' if(sys == 'kalkyl'){ } if(sys == 'local'){ - basedir = '~/Dropbox/postdoc/projects/ase_msc_benchmarking/data/gsnap' + basedir = '~/Documents/postdoc/ase/asebenchmark/novercontrol/data/gsnap' varcalls.datadir = file.path(basedir, 'varcalls') basecount.datadir = file.path(basedir, 'basecount') annovar.datadir = file.path(basedir, 'annovar') ase.datadir = file.path(basedir, 'ase') resdir = '~/Dropbox/postdoc/projects/ase_msc_benchmarking/res/gsnap/ase' + dbsnpdir = '~/Documents/postdoc/ase/asebenchmark/novercontrol/data/annot/dbsnp' } #mpileup calls @@ -22,11 +23,14 @@ vars.regex = '.*hetvars$' #basecounts bc.regex = '.*basecount$' -#annovar +#annot annovar.file = file.path(annovar.datadir, 'annovar_all_variants.RData') +dbsnp.file = file.path(dbsnpdir, 'dbsnp.v135.common.pos.txt') + #Res bc.file = file.path(basecount.datadir, 'bc.RData') +bc.dbsnp.file = file.path(basecount.datadir, 'bc.dbsnp.RData') ase.res.file = file.path(ase.datadir, 'ase.res.RData') ase.sample.cond.diff.file = file.path(ase.datadir, 'induced.ase.res.RData') induced.ase.sig.file = file.path(ase.datadir, 'induced.ase.sig.RData') @@ -54,7 +58,7 @@ main <- function(){ #Dump save(basecount.list, file = bc.file) - + ### #Filter @@ -69,6 +73,19 @@ main <- function(){ #total number of vars length(unique(unlist(lapply(bc.filt, '[[', 'site')))) #97628 + #Filter on dbSNP presence + dbsnp.vars = read.table(dbsnp.file, stringsAsFactors = FALSE) + colnames(dbsnp.vars) = 'site' + bc.filt = lapply(bc.filt, function(jsample.vars, dbsnp.vars){merge(jsample.vars, dbsnp.vars, by = 'site')}, dbsnp.vars = dbsnp.vars) + + #total number of vars + length(unique(unlist(lapply(bc.filt, '[[', 'site')))) #77969 + + #Dump + save(bc.filt, file = bc.dbsnp.file) + + + ### #DP dist ### From 14a5f37c68b29eba2fda3c6e8233479bc79f6450 Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Fri, 12 Oct 2012 10:17:50 +0200 Subject: [PATCH 4/7] Added FDR estimation for the case of not applying the filter of alternative allele direction --- ase/fdr.R | 74 +++++++++++++++++++++++++-------------------------- ase/snv.ase.R | 19 +++++++------ 2 files changed, 45 insertions(+), 48 deletions(-) diff --git a/ase/fdr.R b/ase/fdr.R index 440eff7..e6ee29c 100644 --- a/ase/fdr.R +++ b/ase/fdr.R @@ -1,71 +1,70 @@ sys = 'local' + if(sys == 'kalk'){ } if(sys == 'local'){ simdir = '~/Dropbox/postdoc/projects/ase/data/sim/synt' snpdir = '~/Documents/postdoc/ase/sim/snps' - exprdir = file.path(simdir, 'expr') - asedir = file.path(simdir, 'ase') - number2sample.map.file = file.path(simdir, 'customref2enstcounts.filemapping.tab') - snp2count.file = file.path(simdir, 'snps', 'snp2count.RData') - snp2count.uniq.file = file.path(simdir, 'snps', 'snp2count.uniq.RData') - ase.res.file = file.path(asedir, 'ase.res.RData') + ase.method.resdir = '~/Dropbox/postdoc/projects/ase/res/sim' + ase.method.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/ase' annovardir = '~/Dropbox/postdoc/projects/ase/data/sim/synt/annovar' - #annovar.tab.file = file.path(annovardir, 'synt.snps.tab.variant_function') - #annovar.file = file.path(annovardir, 'vars.refgene.annot.RData') - annovar.file = file.path(annovardir, 'synt.snps.geneannot.RData') - vars.annot.list.file = file.path(asedir, 'ase.res.annot.RData') - genes.pass.file = file.path(asedir, 'multisnp.genes.pass.RData') - nonuniq.genes.file = file.path(asedir, 'nonuniq.genes.RData') - induced.ase.res.file = file.path(asedir, 'induced.ase.res.RData') - induced.ase.sig.file = file.path(asedir, 'induced.ase.sig.RData') } +exprdir = file.path(simdir, 'expr') +asedir = file.path(simdir, 'ase') + +ase.res.file = file.path(asedir, 'ase.res.RData') +annovar.file = file.path(annovardir, 'synt.snps.geneannot.RData') +vars.annot.list.file = file.path(asedir, 'ase.res.annot.RData') +genes.pass.file = file.path(asedir, 'multisnp.genes.pass.RData') +nonuniq.genes.file = file.path(asedir, 'nonuniq.genes.RData') +induced.ase.res.file = file.path(asedir, 'induced.ase.res.RData') +induced.ase.sig.file = file.path(asedir, 'induced.ase.sig.RData') + main <- function(){ + ###### #1. Condition-INDEPENDENT ASE ###### - ### #1.1 Variants ### - #Without test for sig diff number of treat vs untreat + + #Load "TRUTH" ase.noninduced.vars.file = file.path(asedir, 'ase.noninduced.vars.RData') - load(ase.noninduced.vars.file) #alt.sig.vars + load(ase.noninduced.vars.file) #sig.vars, alt.sig.vars + true.sig.vars = unique(unlist(lapply(sig.vars, '[[', 'chrpos'))) true.alt.sig.vars = alt.sig.vars - resdir = '~/Dropbox/postdoc/projects/ase/res/sim' - ase.cond.diff.file = file.path(resdir, 'ase.cond.diff.tab') - load(ase.cond.diff.file) #alt.sig.vars - - #get fdr + + #Load significant variants after application of ASE method + sig.ase.res.file = file.path(ase.method.resdir, 'sig.ase.res.RData') + load(sig.ase.res.file) #sig.vars, alt.sig.vars + sig.vars = unique(unlist(lapply(sig.vars, '[[', 'chrpos'))) + + #get fdr, alt allele direction filtered: NO + fdr = get.fdr(sig.vars, true.sig.vars) + print(fdr) # + + #get fdr, alt allele direction filtered: YES fdr = get.fdr(alt.sig.vars, true.alt.sig.vars) print(fdr) #29%, 6049, 21092 - - #With test for sig diff treat vs untreat - ase.noninduced.vars.file = file.path(asedir, 'ase.noninduced.vars.RData') - load(ase.noninduced.vars.file) - true.sig.alt.ase = sig.alt.ase - resdir = '~/Dropbox/postdoc/projects/ase/res/sim' - ase.cond.diff.file = file.path(resdir, 'ase.cond.diff.tab') - load(ase.cond.diff.file) #sig.alt.ase - #get fdr - fdr = get.fdr(rownames(sig.alt.ase), rownames(true.sig.alt.ase)) - print(fdr) #96%, 25, 26 ### #1.2 Genes ### - #load data + + #Load "TRUTH" load(genes.pass.file) true.genes.pass = genes.pass true.gene2nsamples = gene2nsamples - run.ase.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/ase' - genes.pass.file = file.path(run.ase.datadir, 'multisnp.genes.pass.RData') + + #Load significant variants after application of ASE method + genes.pass.file = file.path(ase.method.datadir, 'multisnp.genes.pass.RData') load(genes.pass.file) #genes.pass #get fdr, n.samples = 2 @@ -87,8 +86,7 @@ main <- function(){ true.genes.pass = genes.pass true.sig.vars.nsamples1 = sig.vars.nsamples1 true.genes.pass.nsamples1 = genes.pass.nsamples1 - ase.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/ase' - induced.ase.sig.file = file.path(ase.datadir, 'induced.ase.sig.RData') + induced.ase.sig.file = file.path(ase.method.datadir, 'induced.ase.sig.RData') load(induced.ase.sig.file) #sig.vars, genes.pass ## diff --git a/ase/snv.ase.R b/ase/snv.ase.R index 9063121..9b58abd 100644 --- a/ase/snv.ase.R +++ b/ase/snv.ase.R @@ -32,8 +32,7 @@ dbsnp.file = file.path(dbsnpdir, 'dbsnp.v135.common.pos.txt') bc.file = file.path(basecount.datadir, 'bc.RData') bc.dbsnp.file = file.path(basecount.datadir, 'bc.dbsnp.RData') ase.res.file = file.path(ase.datadir, 'ase.res.RData') -ase.sample.cond.diff.file = file.path(ase.datadir, 'induced.ase.res.RData') -induced.ase.sig.file = file.path(ase.datadir, 'induced.ase.sig.RData') +sig.ase.res.file = file.path(resdir, 'sig.ase.res.RData') main <- function(){ @@ -133,7 +132,14 @@ main <- function(){ #Get sig hits sig.vars = lapply(ase.res, function(sample.vars, alpha, pval.col){sample.vars = sample.vars[which(sample.vars[, pval.col] <= alpha), ]; return(sample.vars);}, alpha = alpha, pval.col = pval.col) - sig.vars.mat = get.sig.mat(ase.res, alpha, pval.col = pval.col) + + #Filter on alternative allele direction + alt.vars = lapply(sig.vars, alt.filter) + alt.sig.vars = unique(unlist(alt.vars)) + length(alt.sig.vars) #21092 + + #Dump + save(sig.vars, alt.sig.vars, file = sig.ase.res.file) ### @@ -147,14 +153,7 @@ main <- function(){ #number of uniq sig.vars length(unique(unlist(lapply(sig.vars, '[[', 'site')))) #45952 - - - - #((see varstats.R for inspiration, or its enough with the stuff in multiplesnsps2gene.R)) - #gene annot - #DE - ### #Gene based analysis. From 9897832980419c69f52ccbb84b447a1fea2916bd Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Fri, 12 Oct 2012 10:34:57 +0200 Subject: [PATCH 5/7] changed a hardcoded path --- ase/fdr.R | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/ase/fdr.R b/ase/fdr.R index e6ee29c..4bc1fc1 100644 --- a/ase/fdr.R +++ b/ase/fdr.R @@ -5,23 +5,21 @@ sys = 'local' if(sys == 'kalk'){ } if(sys == 'local'){ - simdir = '~/Dropbox/postdoc/projects/ase/data/sim/synt' - snpdir = '~/Documents/postdoc/ase/sim/snps' + syntdir = '/proj/b2012046/edsgard/ase/sim/data/synt/ase' ase.method.resdir = '~/Dropbox/postdoc/projects/ase/res/sim' - ase.method.datadir = '~/Dropbox/postdoc/projects/ase/data/sim/ase' - annovardir = '~/Dropbox/postdoc/projects/ase/data/sim/synt/annovar' } -exprdir = file.path(simdir, 'expr') -asedir = file.path(simdir, 'ase') +#Variant-based analysis +synt.sig.ase.res.file = file.path(syntdir, 'ase.noninduced.vars.RData') +sig.ase.res.file = file.path(ase.method.resdir, 'sig.ase.res.RData') + +#Gene-based analysis +synt.genes.pass.file = file.path(syntdir, 'multisnp.genes.pass.RData') +method.genes.pass.file = file.path(ase.method.resdir, 'multisnp.genes.pass.RData') + +#Cond-dep ASE +induced.ase.sig.file = file.path(syntdir, 'induced.ase.sig.RData') -ase.res.file = file.path(asedir, 'ase.res.RData') -annovar.file = file.path(annovardir, 'synt.snps.geneannot.RData') -vars.annot.list.file = file.path(asedir, 'ase.res.annot.RData') -genes.pass.file = file.path(asedir, 'multisnp.genes.pass.RData') -nonuniq.genes.file = file.path(asedir, 'nonuniq.genes.RData') -induced.ase.res.file = file.path(asedir, 'induced.ase.res.RData') -induced.ase.sig.file = file.path(asedir, 'induced.ase.sig.RData') main <- function(){ @@ -35,19 +33,17 @@ main <- function(){ ### #Load "TRUTH" - ase.noninduced.vars.file = file.path(asedir, 'ase.noninduced.vars.RData') - load(ase.noninduced.vars.file) #sig.vars, alt.sig.vars + load(synt.sig.ase.res.file) #sig.vars, alt.sig.vars true.sig.vars = unique(unlist(lapply(sig.vars, '[[', 'chrpos'))) true.alt.sig.vars = alt.sig.vars #Load significant variants after application of ASE method - sig.ase.res.file = file.path(ase.method.resdir, 'sig.ase.res.RData') load(sig.ase.res.file) #sig.vars, alt.sig.vars sig.vars = unique(unlist(lapply(sig.vars, '[[', 'chrpos'))) #get fdr, alt allele direction filtered: NO fdr = get.fdr(sig.vars, true.sig.vars) - print(fdr) # + print(fdr) #tbd #get fdr, alt allele direction filtered: YES fdr = get.fdr(alt.sig.vars, true.alt.sig.vars) @@ -59,13 +55,12 @@ main <- function(){ ### #Load "TRUTH" - load(genes.pass.file) + load(synt.genes.pass.file) #genes.pass, gene2nsamples true.genes.pass = genes.pass true.gene2nsamples = gene2nsamples #Load significant variants after application of ASE method - genes.pass.file = file.path(ase.method.datadir, 'multisnp.genes.pass.RData') - load(genes.pass.file) #genes.pass + load(method.genes.pass.file) #genes.pass, gene2nsamples #get fdr, n.samples = 2 fdr = get.fdr(genes.pass, true.genes.pass) @@ -86,7 +81,7 @@ main <- function(){ true.genes.pass = genes.pass true.sig.vars.nsamples1 = sig.vars.nsamples1 true.genes.pass.nsamples1 = genes.pass.nsamples1 - induced.ase.sig.file = file.path(ase.method.datadir, 'induced.ase.sig.RData') + induced.ase.sig.file = file.path(ase.method.resdir, 'induced.ase.sig.RData') load(induced.ase.sig.file) #sig.vars, genes.pass ## From f631ea0304726152a8e88294345c4f9cf01c7b4d Mon Sep 17 00:00:00 2001 From: Daniel Edsgard Date: Fri, 12 Oct 2012 14:16:51 +0200 Subject: [PATCH 6/7] Included handling of "chr" prefix --- ase/fdr.R | 20 ++++++++++++-------- ase/snv.ase.R | 6 +++--- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/ase/fdr.R b/ase/fdr.R index 4bc1fc1..421266b 100644 --- a/ase/fdr.R +++ b/ase/fdr.R @@ -1,12 +1,14 @@ - +sys = 'kalk' sys = 'local' if(sys == 'kalk'){ + syntdir = '/proj/b2012046/edsgard/ase/sim/data/synt/ase' + ase.method.resdir = 'FILLMEIN' } if(sys == 'local'){ - syntdir = '/proj/b2012046/edsgard/ase/sim/data/synt/ase' - ase.method.resdir = '~/Dropbox/postdoc/projects/ase/res/sim' + syntdir = '~/Dropbox/postdoc/projects/ase/data/sim/synt/ase' + ase.method.resdir = '~/Documents/postdoc/ase/asebenchmark/novercontrol/data/gsnap/ase' } #Variant-based analysis @@ -39,15 +41,17 @@ main <- function(){ #Load significant variants after application of ASE method load(sig.ase.res.file) #sig.vars, alt.sig.vars - sig.vars = unique(unlist(lapply(sig.vars, '[[', 'chrpos'))) - + sig.vars = unique(unlist(lapply(sig.vars, rownames))) + sig.vars = gsub('^chr', '', sig.vars) + alt.sig.vars = gsub('^chr', '', alt.sig.vars) + #get fdr, alt allele direction filtered: NO fdr = get.fdr(sig.vars, true.sig.vars) - print(fdr) #tbd + print(fdr) #7.4% #get fdr, alt allele direction filtered: YES fdr = get.fdr(alt.sig.vars, true.alt.sig.vars) - print(fdr) #29%, 6049, 21092 + print(fdr) #6.3% ### @@ -61,7 +65,7 @@ main <- function(){ #Load significant variants after application of ASE method load(method.genes.pass.file) #genes.pass, gene2nsamples - + #get fdr, n.samples = 2 fdr = get.fdr(genes.pass, true.genes.pass) print(fdr) #39%, 146, 373 diff --git a/ase/snv.ase.R b/ase/snv.ase.R index 9b58abd..1eedcad 100644 --- a/ase/snv.ase.R +++ b/ase/snv.ase.R @@ -32,7 +32,7 @@ dbsnp.file = file.path(dbsnpdir, 'dbsnp.v135.common.pos.txt') bc.file = file.path(basecount.datadir, 'bc.RData') bc.dbsnp.file = file.path(basecount.datadir, 'bc.dbsnp.RData') ase.res.file = file.path(ase.datadir, 'ase.res.RData') -sig.ase.res.file = file.path(resdir, 'sig.ase.res.RData') +sig.ase.res.file = file.path(ase.datadir, 'sig.ase.res.RData') main <- function(){ @@ -121,7 +121,7 @@ main <- function(){ #Dump save(ase.res, file = ase.res.file) - + ### #Sig hits stats @@ -136,7 +136,7 @@ main <- function(){ #Filter on alternative allele direction alt.vars = lapply(sig.vars, alt.filter) alt.sig.vars = unique(unlist(alt.vars)) - length(alt.sig.vars) #21092 + length(alt.sig.vars) #17832 #Dump save(sig.vars, alt.sig.vars, file = sig.ase.res.file) From 96ee7c815a16a5c4332f9cbc509166f7c1e2a2b4 Mon Sep 17 00:00:00 2001 From: edajeda Date: Mon, 29 Oct 2012 15:04:05 +0100 Subject: [PATCH 7/7] Added sensitivity calculation to fdr.R --- ase/fdr.R | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/ase/fdr.R b/ase/fdr.R index 421266b..cd4bfbc 100644 --- a/ase/fdr.R +++ b/ase/fdr.R @@ -1,4 +1,3 @@ - sys = 'kalk' sys = 'local' @@ -53,6 +52,9 @@ main <- function(){ fdr = get.fdr(alt.sig.vars, true.alt.sig.vars) print(fdr) #6.3% + #get sensitivity + sens = get.sense(alt.sig.vars, true.alt.sig.vars) + print(sens) ### #1.2 Genes @@ -178,6 +180,17 @@ get.fdr <- function(feat.pass, true.feat.pass){ return(fdr) } +get.sens <- function(feat.pass, true.feat.pass){ + tp = intersect(feat.pass, true.feat.pass) + fn = setdiff(true.feat.pass, feat.pass) + n.tp = length(tp) + n.fn = length(fn) + sens = n.tp / (n.tp + n.fn) + + sens = c(sens, n.tp, n.fn) + return(sens) +} + filter.altallele <- function(vars, frac.thr = 0.5, gene.col = 'gene annot', frac.col = 'frac', pval.col = 'mod.padj', min.alt = 1, alpha = 0.05){ sig.alt.vars = vars[which(vars[, pval.col] < alpha & vars[, frac.col] > frac.thr), ]