From 902bbd23ae9b23a54f02ce086348ead76b37daca Mon Sep 17 00:00:00 2001 From: Jeremiah Wala Date: Tue, 7 Jul 2015 12:09:44 -0400 Subject: [PATCH] added BamQCPlot.R --- R/BamQCPlot.R | 124 ++++++++++++++++++++++++++++++++++++ src/VariantBamWalker.h | 3 - src/autom4te.cache/output.0 | 1 + src/autom4te.cache/output.1 | 1 + src/autom4te.cache/traces.0 | 36 +++++------ src/autom4te.cache/traces.1 | 32 +++++----- src/configure | 1 + src/configure.ac | 1 + src/variant.cpp | 21 +++++- 9 files changed, 180 insertions(+), 40 deletions(-) create mode 100644 R/BamQCPlot.R diff --git a/R/BamQCPlot.R b/R/BamQCPlot.R new file mode 100644 index 0000000..4ed5c04 --- /dev/null +++ b/R/BamQCPlot.R @@ -0,0 +1,124 @@ +library(optparse) +#RLIBDIR = '/cga/meyerson/home/marcin/Software/R/x86_64-unknown-linux-gnu-library/3.1/' +#GIT.HOME = '/cga/meyerson/home/marcin/DB/git' + +option_list = list( + make_option(c("-i", "--input"), type = "character", default = "qcreport.txt", help = "Input txt file from a snowman preprocess qcreport.txt"), + make_option(c("-o", "--output"), type = "character", default = "qcreport.pdf", help = "Output pdf to generate") +) + +parseobj = OptionParser(option_list=option_list) +opt = parse_args(parseobj) + +if (is.null(opt$input)) + stop(print_help(parseobj)) + +print(getwd()) +if (!file.exists(opt$input)) + stop(paste("Input file does not exist", opt$input)) + +print(opt) + +require(ggplot2) +require(reshape2) +require(gridExtra) + +## read the table +con <- file(opt$input, open = "r") + +df.mapq <- df.nm <- df.isize <- df.as <- df.xp <- df.len <- df.phred <- df.clip <- data.frame() +rg <- list() + +## function for parsing histogram strings +function(ll, id) { + ss <- strsplit(ll[id],',')[[1]] + df <- data.frame(matrix(unlist(strsplit(ss, '_')), nrow=length(ss), byrow=T)) + colnames(df) <- c("Min", "Max", "Count") + return (df) +} + +while (length(line <- readLines(con, n = 1, warn = FALSE)) > 0) { + + ll <- strsplit(line, '\t')[[1]] + + thisrg = ll[1] + rg[[thisrg]] <- data.frame(readgroup = thisrg, ReadCount = as.numeric(ll[2]), Supplementary = as.numeric(ll[3]), + Unmapped = as.numeric(ll[4]), MateUnmapped = as.numeric(ll[5]), + QCFailed = as.numeric(ll[6]), Duplicate = as.numeric(ll[7])) + + + ss <- strsplit(ll[8],',')[[1]] + df.mapq <- data.frame(matrix(unlist(strsplit(ss, '_')), nrow=length(ss), byrow=T)) + colnames(df.mapq) <- c("Min", "Max", "Count") + + ## found a new one + #if (grepl("ReadGroup", line)) { + # thisrg = gsub("READGROUP:BI:(.*)", "\\1", line) + # rg[[thisrg]] <- data.frame(readgroup = thisrg) + #} else if (grepl("total", line)) { + # rg[[thisrg]]$total = as.numeric(gsub("total,([0-9]+)", "\\1", line)) + #} else if (grepl("unmap", line)) { + # rg[[thisrg]]$unmap = as.numeric(gsub("unmap,([0-9]+)", "\\1", line)) + #} else if (grepl("qcfail", line)) { + # rg[[thisrg]]$qcfail = as.numeric(gsub("qcfail,([0-9]+)", "\\1", line)) + #} else if (grepl("duplicate", line)) { + # rg[[thisrg]]$duplicate = as.numeric(gsub("duplicate,([0-9]+)", "\\1", line)) + #} else if (grepl("supplementary", line)) { + # rg[[thisrg]]$supp = as.numeric(gsub("supplementary,([0-9]+)", "\\1", line)) + #} else if (grepl("mapq", line)) { + # df.mapq <- rbind(df.mapq, as.numeric(strsplit(gsub("mapq,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("nm", line)) { + # df.nm <- rbind(df.nm, as.numeric(strsplit(gsub("nm,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("isize", line)) { + # df.isize <- rbind(df.isize, as.numeric(strsplit(gsub("isize,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("as", line)) { + # df.as <- rbind(df.as, as.numeric(strsplit(gsub("as,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("xp", line)) { + # df.xp <- rbind(df.xp, as.numeric(strsplit(gsub("xp,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("clip", line)) { + # df.clip <- rbind(df.clip, as.numeric(strsplit(gsub("clip,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("len", line)) { + # df.len <- rbind(df.len, as.numeric(strsplit(gsub("len,([0-9]+)", "\\1", line), ",")[[1]])) + #} else if (grepl("phred", line)) { + # df.phred <- rbind(df.phred, as.numeric(strsplit(gsub("phred,([0-9]+)", "\\1", line), ",")[[1]])) + #} else { + # stop(paste("Failed to read file at line:", line)) + #} + +} + +close(con) + +colnames(df.mapq) <- seq(from=0,to=60) +colnames(df.nm) <- seq(from=0,to=ncol(df.nm)-1) +colnames(df.isize) <- seq(from=0,to=ncol(df.isize)-1) +colnames(df.xp) <- seq(from=0, to=ncol(df.xp)-1) +colnames(df.as) <- seq(from=0, to=ncol(df.as)-1) +colnames(df.len) <- seq(from=0, to=ncol(df.len)-1) +colnames(df.phred) <- seq(from=0, to=ncol(df.phred)-1) +colnames(df.clip) <- seq(from=0, to=ncol(df.clip)-1) +readg <- sapply(rg, function(x) x$readgroup) + +df.mapq$readgroup <- readg +df.nm$readgroup <- readg +df.isize$readgroup <- readg +df.xp$readgroup <- readg +df.as$readgroup <- readg +df.phred$readgroup <- readg +df.len$readgroup <- readg +df.clip$readgroup <- readg + +g.mapq <- ggplot(df.mapq.m <- melt(df.mapq, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(30,61)) + ylab('Reads') + xlab('Mapping Quality') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.nm <- ggplot(df.nm.m <- melt(df.nm, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,50)) + ylab('Reads') + xlab('NM Tag') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.isize <- ggplot(df.isize.m <- melt(df.isize, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(10,2001)) + ylab('Reads') + xlab('InsertSize') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.xp <- ggplot(df.xp.m <- melt(df.xp, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,100)) + ylab('Reads') + xlab('XP Tag') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.as <- ggplot(df.as.m <- melt(df.as, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,100)) + ylab('Reads') + xlab('AS Tag') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.len <- ggplot(df.len.m <- melt(df.len, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(20,102)) + ylab('Reads') + xlab('Read Length') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.clip <- ggplot(df.clip.m <- melt(df.clip, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,85)) + ylab('Reads') + xlab('Clipped bases') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) +g.phred <- ggplot(df.phred.m <- melt(df.phred, id='readgroup'), aes(x=as.numeric(variable), y=pmax(log(value,10),0), group=readgroup, color=readgroup)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,43)) + ylab('Reads') + xlab('Mean read Phred quality') + scale_y_continuous(breaks=seq(0,8), labels=parse(text=paste('10', seq(0,8), sep='^'))) + + +pdf(opt$output, width=15, height=8) +print(grid.arrange(g.mapq, g.nm, g.isize, g.len, g.clip, g.phred, ncol=3)) +dev.off() + diff --git a/src/VariantBamWalker.h b/src/VariantBamWalker.h index 1843bb2..f36d7f5 100644 --- a/src/VariantBamWalker.h +++ b/src/VariantBamWalker.h @@ -18,11 +18,8 @@ class VariantBamWalker: public SnowTools::BamWalker void TrackSeenRead(SnowTools::BamRead &r); void printMessage(const SnowTools::ReadCount &rc_main, const SnowTools::BamRead &r) const; - - private: SnowTools::BamStats m_stats; - }; #endif diff --git a/src/autom4te.cache/output.0 b/src/autom4te.cache/output.0 index 1bf1ca2..657f0ee 100644 --- a/src/autom4te.cache/output.0 +++ b/src/autom4te.cache/output.0 @@ -4412,6 +4412,7 @@ fi # Set compiler flags. AM_CXXFLAGS="-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas" +##AC_SUBST(AM_CXXFLAGS, "-g -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas") CXXFLAGS="-O3" CFLAGS="-O3" diff --git a/src/autom4te.cache/output.1 b/src/autom4te.cache/output.1 index 1bf1ca2..657f0ee 100644 --- a/src/autom4te.cache/output.1 +++ b/src/autom4te.cache/output.1 @@ -4412,6 +4412,7 @@ fi # Set compiler flags. AM_CXXFLAGS="-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas" +##AC_SUBST(AM_CXXFLAGS, "-g -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas") CXXFLAGS="-O3" CFLAGS="-O3" diff --git a/src/autom4te.cache/traces.0 b/src/autom4te.cache/traces.0 index 84a5949..854aec3 100644 --- a/src/autom4te.cache/traces.0 +++ b/src/autom4te.cache/traces.0 @@ -685,25 +685,25 @@ m4trace:configure.ac:20: -1- AC_DEFINE_TRACE_LITERAL([HAVE_CLOCK_GETTIME]) m4trace:configure.ac:20: -1- AH_OUTPUT([HAVE_CLOCK_GETTIME], [/* clock_getttime found */ #undef HAVE_CLOCK_GETTIME]) m4trace:configure.ac:107: -1- AC_SUBST([AM_CXXFLAGS], ["-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas"]) -m4trace:configure.ac:108: -1- AC_SUBST([CXXFLAGS], ["-O3"]) -m4trace:configure.ac:109: -1- AC_SUBST([CFLAGS], ["-O3"]) -m4trace:configure.ac:110: -1- AC_SUBST([CPPFLAGS], ["$CPPFLAGS $hts_include $snowtools_include $boost_include"]) -m4trace:configure.ac:111: -1- AC_SUBST([LDFLAGS], ["$snowtools_ldflags $hts_ldflags $LDFLAGS -L."]) -m4trace:configure.ac:113: -1- AC_SUBST([LIBS], ["$LIBS $snowtools_libs -lhts -lpthread"]) -m4trace:configure.ac:116: -1- AC_CHECK_HEADERS([htslib/hts.h], [], [AC_MSG_ERROR([The HTS library must be installed (https://github.com/samtools/htslib). You can specify its path with the --with-hts=PATH option])]) -m4trace:configure.ac:116: -1- AH_OUTPUT([HAVE_HTSLIB_HTS_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:109: -1- AC_SUBST([CXXFLAGS], ["-O3"]) +m4trace:configure.ac:110: -1- AC_SUBST([CFLAGS], ["-O3"]) +m4trace:configure.ac:111: -1- AC_SUBST([CPPFLAGS], ["$CPPFLAGS $hts_include $snowtools_include $boost_include"]) +m4trace:configure.ac:112: -1- AC_SUBST([LDFLAGS], ["$snowtools_ldflags $hts_ldflags $LDFLAGS -L."]) +m4trace:configure.ac:114: -1- AC_SUBST([LIBS], ["$LIBS $snowtools_libs -lhts -lpthread"]) +m4trace:configure.ac:117: -1- AC_CHECK_HEADERS([htslib/hts.h], [], [AC_MSG_ERROR([The HTS library must be installed (https://github.com/samtools/htslib). You can specify its path with the --with-hts=PATH option])]) +m4trace:configure.ac:117: -1- AH_OUTPUT([HAVE_HTSLIB_HTS_H], [/* Define to 1 if you have the header file. */ #undef HAVE_HTSLIB_HTS_H]) -m4trace:configure.ac:119: -1- AC_CHECK_HEADERS([boost/icl/interval_set.hpp], [], [AC_MSG_WARN([The Boost library must be installed for a few of the interval operations. NOT REQUIRED for VariantBam. Specify its path with the --with-boostxs=PATH option])]) -m4trace:configure.ac:119: -1- AH_OUTPUT([HAVE_BOOST_ICL_INTERVAL_SET_HPP], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:120: -1- AC_CHECK_HEADERS([boost/icl/interval_set.hpp], [], [AC_MSG_WARN([The Boost library must be installed for a few of the interval operations. NOT REQUIRED for VariantBam. Specify its path with the --with-boostxs=PATH option])]) +m4trace:configure.ac:120: -1- AH_OUTPUT([HAVE_BOOST_ICL_INTERVAL_SET_HPP], [/* Define to 1 if you have the header file. */ #undef HAVE_BOOST_ICL_INTERVAL_SET_HPP]) -m4trace:configure.ac:122: -1- AC_CHECK_HEADERS([SnowTools/SnowTools.h], [], [AC_MSG_ERROR([The SnowTools library must be installed. You can specify its path with the --with-snowtools=PATH option])]) -m4trace:configure.ac:122: -1- AH_OUTPUT([HAVE_SNOWTOOLS_SNOWTOOLS_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:123: -1- AC_CHECK_HEADERS([SnowTools/SnowTools.h], [], [AC_MSG_ERROR([The SnowTools library must be installed. You can specify its path with the --with-snowtools=PATH option])]) +m4trace:configure.ac:123: -1- AH_OUTPUT([HAVE_SNOWTOOLS_SNOWTOOLS_H], [/* Define to 1 if you have the header file. */ #undef HAVE_SNOWTOOLS_SNOWTOOLS_H]) -m4trace:configure.ac:125: -1- AC_CHECK_HEADERS([ahocorasick/ahocorasick.h], [], [AC_MSG_WARN([The Aho-Corasick library must be installed (http://sourceforge.net/projects/multifast/), to do string matching rules. Can specify --with-aho=PATH option])]) -m4trace:configure.ac:125: -1- AH_OUTPUT([HAVE_AHOCORASICK_AHOCORASICK_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:126: -1- AC_CHECK_HEADERS([ahocorasick/ahocorasick.h], [], [AC_MSG_WARN([The Aho-Corasick library must be installed (http://sourceforge.net/projects/multifast/), to do string matching rules. Can specify --with-aho=PATH option])]) +m4trace:configure.ac:126: -1- AH_OUTPUT([HAVE_AHOCORASICK_AHOCORASICK_H], [/* Define to 1 if you have the header file. */ #undef HAVE_AHOCORASICK_AHOCORASICK_H]) -m4trace:configure.ac:127: -1- AC_CONFIG_FILES([Makefile]) -m4trace:configure.ac:129: -1- AC_SUBST([LIB@&t@OBJS], [$ac_libobjs]) -m4trace:configure.ac:129: -1- AC_SUBST([LTLIBOBJS], [$ac_ltlibobjs]) -m4trace:configure.ac:129: -1- _AC_AM_CONFIG_HEADER_HOOK([$ac_file]) -m4trace:configure.ac:129: -1- _AM_OUTPUT_DEPENDENCY_COMMANDS +m4trace:configure.ac:128: -1- AC_CONFIG_FILES([Makefile]) +m4trace:configure.ac:130: -1- AC_SUBST([LIB@&t@OBJS], [$ac_libobjs]) +m4trace:configure.ac:130: -1- AC_SUBST([LTLIBOBJS], [$ac_ltlibobjs]) +m4trace:configure.ac:130: -1- _AC_AM_CONFIG_HEADER_HOOK([$ac_file]) +m4trace:configure.ac:130: -1- _AM_OUTPUT_DEPENDENCY_COMMANDS diff --git a/src/autom4te.cache/traces.1 b/src/autom4te.cache/traces.1 index 0984a73..3a53c9d 100644 --- a/src/autom4te.cache/traces.1 +++ b/src/autom4te.cache/traces.1 @@ -161,23 +161,23 @@ m4trace:configure.ac:20: -1- AC_DEFINE_TRACE_LITERAL([HAVE_CLOCK_GETTIME]) m4trace:configure.ac:20: -1- AH_OUTPUT([HAVE_CLOCK_GETTIME], [/* clock_getttime found */ #undef HAVE_CLOCK_GETTIME]) m4trace:configure.ac:107: -1- AC_SUBST([AM_CXXFLAGS], ["-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas"]) -m4trace:configure.ac:108: -1- AC_SUBST([CXXFLAGS], ["-O3"]) -m4trace:configure.ac:109: -1- AC_SUBST([CFLAGS], ["-O3"]) -m4trace:configure.ac:110: -1- AC_SUBST([CPPFLAGS], ["$CPPFLAGS $hts_include $snowtools_include $boost_include"]) -m4trace:configure.ac:111: -1- AC_SUBST([LDFLAGS], ["$snowtools_ldflags $hts_ldflags $LDFLAGS -L."]) -m4trace:configure.ac:113: -1- AC_SUBST([LIBS], ["$LIBS $snowtools_libs -lhts -lpthread"]) -m4trace:configure.ac:116: -1- AC_CHECK_HEADERS([htslib/hts.h], [], [AC_MSG_ERROR([The HTS library must be installed (https://github.com/samtools/htslib). You can specify its path with the --with-hts=PATH option])]) -m4trace:configure.ac:116: -1- AH_OUTPUT([HAVE_HTSLIB_HTS_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:109: -1- AC_SUBST([CXXFLAGS], ["-O3"]) +m4trace:configure.ac:110: -1- AC_SUBST([CFLAGS], ["-O3"]) +m4trace:configure.ac:111: -1- AC_SUBST([CPPFLAGS], ["$CPPFLAGS $hts_include $snowtools_include $boost_include"]) +m4trace:configure.ac:112: -1- AC_SUBST([LDFLAGS], ["$snowtools_ldflags $hts_ldflags $LDFLAGS -L."]) +m4trace:configure.ac:114: -1- AC_SUBST([LIBS], ["$LIBS $snowtools_libs -lhts -lpthread"]) +m4trace:configure.ac:117: -1- AC_CHECK_HEADERS([htslib/hts.h], [], [AC_MSG_ERROR([The HTS library must be installed (https://github.com/samtools/htslib). You can specify its path with the --with-hts=PATH option])]) +m4trace:configure.ac:117: -1- AH_OUTPUT([HAVE_HTSLIB_HTS_H], [/* Define to 1 if you have the header file. */ #undef HAVE_HTSLIB_HTS_H]) -m4trace:configure.ac:119: -1- AC_CHECK_HEADERS([boost/icl/interval_set.hpp], [], [AC_MSG_WARN([The Boost library must be installed for a few of the interval operations. NOT REQUIRED for VariantBam. Specify its path with the --with-boostxs=PATH option])]) -m4trace:configure.ac:119: -1- AH_OUTPUT([HAVE_BOOST_ICL_INTERVAL_SET_HPP], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:120: -1- AC_CHECK_HEADERS([boost/icl/interval_set.hpp], [], [AC_MSG_WARN([The Boost library must be installed for a few of the interval operations. NOT REQUIRED for VariantBam. Specify its path with the --with-boostxs=PATH option])]) +m4trace:configure.ac:120: -1- AH_OUTPUT([HAVE_BOOST_ICL_INTERVAL_SET_HPP], [/* Define to 1 if you have the header file. */ #undef HAVE_BOOST_ICL_INTERVAL_SET_HPP]) -m4trace:configure.ac:122: -1- AC_CHECK_HEADERS([SnowTools/SnowTools.h], [], [AC_MSG_ERROR([The SnowTools library must be installed. You can specify its path with the --with-snowtools=PATH option])]) -m4trace:configure.ac:122: -1- AH_OUTPUT([HAVE_SNOWTOOLS_SNOWTOOLS_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:123: -1- AC_CHECK_HEADERS([SnowTools/SnowTools.h], [], [AC_MSG_ERROR([The SnowTools library must be installed. You can specify its path with the --with-snowtools=PATH option])]) +m4trace:configure.ac:123: -1- AH_OUTPUT([HAVE_SNOWTOOLS_SNOWTOOLS_H], [/* Define to 1 if you have the header file. */ #undef HAVE_SNOWTOOLS_SNOWTOOLS_H]) -m4trace:configure.ac:125: -1- AC_CHECK_HEADERS([ahocorasick/ahocorasick.h], [], [AC_MSG_WARN([The Aho-Corasick library must be installed (http://sourceforge.net/projects/multifast/), to do string matching rules. Can specify --with-aho=PATH option])]) -m4trace:configure.ac:125: -1- AH_OUTPUT([HAVE_AHOCORASICK_AHOCORASICK_H], [/* Define to 1 if you have the header file. */ +m4trace:configure.ac:126: -1- AC_CHECK_HEADERS([ahocorasick/ahocorasick.h], [], [AC_MSG_WARN([The Aho-Corasick library must be installed (http://sourceforge.net/projects/multifast/), to do string matching rules. Can specify --with-aho=PATH option])]) +m4trace:configure.ac:126: -1- AH_OUTPUT([HAVE_AHOCORASICK_AHOCORASICK_H], [/* Define to 1 if you have the header file. */ #undef HAVE_AHOCORASICK_AHOCORASICK_H]) -m4trace:configure.ac:127: -1- AC_CONFIG_FILES([Makefile]) -m4trace:configure.ac:129: -1- AC_SUBST([LIB@&t@OBJS], [$ac_libobjs]) -m4trace:configure.ac:129: -1- AC_SUBST([LTLIBOBJS], [$ac_ltlibobjs]) +m4trace:configure.ac:128: -1- AC_CONFIG_FILES([Makefile]) +m4trace:configure.ac:130: -1- AC_SUBST([LIB@&t@OBJS], [$ac_libobjs]) +m4trace:configure.ac:130: -1- AC_SUBST([LTLIBOBJS], [$ac_ltlibobjs]) diff --git a/src/configure b/src/configure index a0e3829..4f49d10 100755 --- a/src/configure +++ b/src/configure @@ -4412,6 +4412,7 @@ fi # Set compiler flags. AM_CXXFLAGS="-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas" +##AC_SUBST(AM_CXXFLAGS, "-g -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas") CXXFLAGS="-O3" CFLAGS="-O3" diff --git a/src/configure.ac b/src/configure.ac index 5fa269b..f5b12cd 100644 --- a/src/configure.ac +++ b/src/configure.ac @@ -105,6 +105,7 @@ fi # Set compiler flags. AC_SUBST(AM_CXXFLAGS, "-g -std=c++11 -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas") +##AC_SUBST(AM_CXXFLAGS, "-g -Wall -Wextra $fail_on_warning -Wno-unknown-pragmas") AC_SUBST(CXXFLAGS, "-O3") AC_SUBST(CFLAGS, "-O3") AC_SUBST(CPPFLAGS, "$CPPFLAGS $hts_include $snowtools_include $boost_include") diff --git a/src/variant.cpp b/src/variant.cpp index ac75961..05fc509 100644 --- a/src/variant.cpp +++ b/src/variant.cpp @@ -53,13 +53,14 @@ namespace opt { static std::string tag_list = ""; static std::string counts_file = ""; static bool counts_only = false; + static std::string bam_qcfile = ""; } enum { OPT_HELP }; -static const char* shortopts = "hvqji:o:r:k:g:Cf:s:ST:l:c:x:"; +static const char* shortopts = "hvji:o:r:k:g:Cf:s:ST:l:c:x:q:Q:"; static const struct option longopts[] = { { "help", no_argument, NULL, OPT_HELP }, { "linked-region", required_argument, NULL, 'l' }, @@ -73,7 +74,8 @@ static const struct option longopts[] = { { "include-header", no_argument, NULL, 'h' }, { "input", required_argument, NULL, 'i' }, { "output-bam", required_argument, NULL, 'o' }, - { "qc-only", no_argument, NULL, 'q' }, + { "qc-file", no_argument, NULL, 'q' }, + { "qc-only-file", no_argument, NULL, 'Q' }, { "rules", required_argument, NULL, 'r' }, { "region", required_argument, NULL, 'g' }, { "region-with-mates", required_argument, NULL, 'c' }, @@ -165,6 +167,9 @@ int main(int argc, char** argv) { if (!opt::counts_only) walk.OpenWriteBam(opt::out); + // if counts or qc only, dont write output + //walk.fop = nullptr; + // print out some info if (opt::verbose) std::cerr << walk << std::endl; @@ -178,6 +183,14 @@ int main(int argc, char** argv) { std::cerr << "...starting filtering" << std::endl; walk.writeVariantBam(); + // dump the stats file + if (opt::bam_qcfile.length()) { + std::ofstream ofs; + ofs.open(opt::bam_qcfile); + ofs << walk.m_stats; + ofs.close(); + } + // display the rule counts walk.MiniRulesToFile(opt::counts_file); @@ -247,6 +260,8 @@ void parseVarOptions(int argc, char** argv) { break; case 'c': arg >> opt::counts_file; break; case 'x': arg >> opt::counts_file; opt::counts_only = true; break; + case 'q': arg >> opt::bam_qcfile; break; + case 'Q': arg >> opt::bam_qcfile; opt::counts_only = true; break; case 'r': { std::string tmp; @@ -302,7 +317,7 @@ void parseVarOptions(int argc, char** argv) { if (opt::bam == "") die = true; - if (opt::out == "") + if (opt::out == "" && !opt::counts_only) opt::to_stdout = true; // dont stop the run for bad bams for quality checking only