diff --git a/.travis.yml b/.travis.yml index 4b712c7c..cf835c18 100644 --- a/.travis.yml +++ b/.travis.yml @@ -26,12 +26,12 @@ script: - cd $TRAVIS_BUILD_DIR # generate autotools's files - bash autotools-init.sh - # run configure + # run configure (without boost checks) - ./configure --prefix=$HOME/IntaRNA --with-vrna=$HOME/miniconda/envs/build-IntaRNA --with-boost=no --without-zlib # compile documentation # - make doxygen-doc # compile, test and install IntaRNA - make -j 2 && make tests -j 2 && make install ##### check IntaRNA build ##### - # run IntaRNA with help output + # run installed IntaRNA with help output - $HOME/IntaRNA/bin/IntaRNA -h diff --git a/ChangeLog b/ChangeLog index db2b1b99..c037bc21 100644 --- a/ChangeLog +++ b/ChangeLog @@ -10,9 +10,59 @@ # changes in development version since last release ################################################################################ + ################################################################################ ################################################################################ + +################################################################################ +### version 3.1.4 +################################################################################ + +# IntaRNA +- bugfix generation and tracing of seeds with bulges and no GU ends +- bugfix seed-extension prediction for seeds with bulges +- noLP for seeds with bulges enabled + +# R +- `IntaRNA_CSV_p-value.R` script to estimate p-values based on energy values +- `IntaRNA_plotRegions.R` = renaming of former `plotRegions.R` + +################################################################################ + +200130 Martin Raden + * IntaRNA/SeedHandlerMfe : + * bugfix generation and tracing of seeds with bulges and no GU ends + * IntaRNA/PredictorMfe*SeedExtension* : + * bugfix enumeration of seeds with bulges + * bin/CommandLineParseing : + * error msgs rephrased + + noLP for seeds with bulges enabled + + setup noLP for seed constraints via outNoLP + * IntaRNA/SeedConstraint : + + isLpAllowed : whether or not lps are allowed in seeds + * IntaRNA/SeedHandlerMfe : + + support for noLP constraint + * test/SeedHandlerMfe : + + test with lp + + test no lp (boundary) + + test no lp (internal) + * test/* + * adaptation to SeedConstraint constructor changes + +200121 Martin Raden + + R/IntaRNA_CSV_p-value.R : former addPvalues2csv.R + + R/IntaRNA_plotRegions.R : former plotRegions.R + - R/addPvalues2csv.R : renamed + - R/plotRegions.R : renamed + * README.md : adapted to renamings + + R/Makefile.am : install R scripts + +191115 Martin Raden + + R/addPvalues2csv.R + * R/README.md : + + docu of addPvalue2csv.R + ################################################################################ ### version 3.1.3 ################################################################################ diff --git a/Makefile.am b/Makefile.am index b9703ef5..43dd1734 100644 --- a/Makefile.am +++ b/Makefile.am @@ -5,7 +5,7 @@ export GCC_COLORS ?= "error=01;31:warning=01;35:note=01;36:caret=01;32:locus=01:quote=01" # sub directories to check for Makefiles -SUBDIRS = src python perl tests doc . +SUBDIRS = src python perl R tests doc . # list of all personalities to be installed #PERSONALITIES = `grep -P "^\\s+case\\s+IntaRNA\\S+\\s*:\\s*return" $(abs_top_srcdir)/src/bin/CommandLineParsing.h | sed "s/^\\s*case\\s\\+\\(IntaRNA\\S*\\)\\s\\+:\\s\\+return.*/\\1/g"` diff --git a/R/IntaRNA_CSV_p-value.R b/R/IntaRNA_CSV_p-value.R new file mode 100644 index 00000000..23d45326 --- /dev/null +++ b/R/IntaRNA_CSV_p-value.R @@ -0,0 +1,150 @@ +#!/usr/bin/env Rscript + +#################################################################### +# Computes p-values and false discovery rates (fdr following Benjamin+Hochberg) +# by fitting a GEV on the energy values computed by IntaRNA. +# Note, such p-value estimates are only useful for genome-wide predictions. +# +# arguments: [ = ] [ = E] +# +# 1 = ";"-separated CSV output of IntaRNA +# 2 = file name to write the extended CSV output to (2 new columns) +# 3 = the column name that holds the energy values to be fitted +# +# example call: +# +# Rscript --vanilla IntaRNA_CSV_p-value.R predictions.csv +# +# This script is part of the IntaRNA source code package. See +# respective licence and documentation for further information. +# +# https://github.com/BackofenLab/IntaRNA +# +#################################################################### + + +#################################################################### +# get command line arguments +#################################################################### + +args = commandArgs(trailingOnly=TRUE) +# check and parse +if (length(args)<1) { stop("call with [ = ] [ = E]", call.=FALSE) } + +# get input file = IntaRNA csv output +intarnaOutputFile = args[1]; +if (!file.exists(intarnaOutputFile )) { stop("intarna-csv-output file '", intarnaOutputFile, "' does not exist!", call.=FALSE) } + +# get output file +outFile = intarnaOutputFile; +if (length(args)>=2) { + outFile = args[2] +} + +# set column to get energies from +colNameE = "E" +# get column name from argument if present +if (length(args)>=3) { + colNameE = args[3] +} + +# column delimiter used in CSV input / output +csvColSep = ";" + +# number of digits of p-values +pValPrec = 7 + +#################################################################### +# fits a generalized extreme value distribution to the given energy data +# adopted from 'gev' function of 'evir' library +# @param energy the IntaRNA energy values to fit (a vector) +# @return the fitting parameters xi, mu, and sigma +gevFitting <- function (energy) +#################################################################### +{ + n.all <- NA + energy <- as.numeric(-energy) + n <- length(energy) + sigma0 <- sqrt(6 * var(energy))/pi + mu0 <- mean(energy) - 0.57722 * sigma0 + xi0 <- 0.1 + theta <- c(xi0, sigma0, mu0) + negloglik <- function(theta, tmp) { + y <- 1 + (theta[1] * (tmp - theta[3]))/theta[2] + if ((theta[2] < 0) || (min(y) < 0)) + out <- 1e+06 + else { + term1 <- length(tmp) * logb(theta[2]) + term2 <- sum((1 + 1/theta[1]) * logb(y)) + term3 <- sum(y^(-1/theta[1])) + out <- term1 + term2 + term3 + } + out + } + # compute fit + fit <- optim(theta, negloglik, hessian = TRUE, tmp = energy) + if (fit$convergence) + warning("gev fit optimization may not have succeeded") + + return( list( xi=fit$par[1], sigma=fit$par[2], mu=fit$par[3] ) ) +} + + +#################################################################### +# computes p-values for the given energy values and GEV distribution +# adopted from 'pgev' function of 'evir' library +# @param energy IntaRNA energy values +# @param gev GEV parameters +# @return p-values for each energy value +gevPvalue <- function (energy, gev=list( xi = 1, mu = 0, sigma = 1) ) +#################################################################### +{ + return ( 1 - exp( - (1 + (gev$xi * ((-energy) - gev$mu))/gev$sigma)^(-1 /gev$xi))) +} + + + +#################################################################### +# parse IntaRNA CSV +#################################################################### + +d = read.csv( intarnaOutputFile, sep=csvColSep ) + +# check if energy column present +if (!is.element(colNameE, colnames(d))) { + stop("'",colNameE,"' is not among the column names of '",intarnaOutputFile,"'", call.=FALSE); +} +# check if unique +if (sum(colnames(d) == colNameE)>1) { + stop("column name '",colNameE,"' occurs more than once in '",intarnaOutputFile,"'", call.=FALSE); +} + + +#################################################################### +# fit p-values +#################################################################### + +# get energies to fit +E = d[,colnames(d) == colNameE] + + +# fit negated energies +gevfit <- gevFitting(E) # fitten + +# get rounded pValue +pVal <- round( gevPvalue( E, gevfit ), digits=pValPrec ) +# get rounded fdr +fdr <- round( p.adjust(pVal, method="BH"), digits=pValPrec ) + +#################################################################### +# write output +#################################################################### + +o = cbind( d, pVal, fdr ) +colnames(o)[ncol(o)-1] = "p-value" +colnames(o)[ncol(o)] = "fdr" + +write.table( o, outFile, sep=csvColSep, row.names=FALSE, col.names = TRUE, quote=FALSE ) + + +#################################################################EOF diff --git a/R/plotRegions.R b/R/IntaRNA_plotRegions.R similarity index 56% rename from R/plotRegions.R rename to R/IntaRNA_plotRegions.R index 3cb2b55e..99048fd2 100644 --- a/R/plotRegions.R +++ b/R/IntaRNA_plotRegions.R @@ -1,151 +1,194 @@ -#!/usr/bin/env Rscript - -#################################################################### -# Visualization of sequence regions covered by RNA-RNA interactions -# predicted by IntaRNA. -# -# arguments: <1|2> -# 1 = ";"-separated CSV output of IntaRNA -# 2 <1|2> = suffix of "start,end,id" CSV cols to plot -# 3 = file name of the otuput figure suffixed -# by one of ".pdf",".png",".svg",".eps",".ps",".jpeg",".tiff" -# -# example call: -# -# Rscript --slave -f plotRegions.R --args pred.csv 1 regions.png -# -# This script is part of the IntaRNA source code package. See -# respective licence and documentation for further information. -# -# https://github.com/BackofenLab/IntaRNA -# -#################################################################### -# check and load dependencies -#################################################################### - -options(warn=-1) -suppressPackageStartupMessages(require(ggplot2)) -suppressPackageStartupMessages(require(ggalt)) -suppressPackageStartupMessages(require(cowplot)) # cowplot starts with a note -options(warn=0) - -theme_set(theme_cowplot()) - -#################################################################### -# get command line arguments -#################################################################### - -args = commandArgs(trailingOnly=TRUE) -# check and parse -if (length(args)!=3) { stop("call with <1|2> ", call.=FALSE) } - -intarnaOutputFile = args[1]; -if (!file.exists(intarnaOutputFile )) { stop("intarna-csv-output file '", intarnaOutputFile, "' does not exist!", call.=FALSE) } - -seqNr = args[2]; -if (seqNr != "1" && seqNr != "2") { stop("second call argument as to be '1' or '2' to specify which regions to plot"); } -# compile column names -id = paste("id",seqNr,sep=""); -start = paste("start",seqNr,sep=""); -end = paste("end",seqNr,sep=""); - -outFile = args[3]; -fileExtensions = c(".pdf",".png",".svg",".eps",".ps",".jpeg",".tiff"); -outFileExtOk = FALSE; -for( ext in fileExtensions ) { outFileExtOk = outFileExtOk || endsWith(outFile,ext); } -if ( !outFileExtOk ) {stop(" has to have one of the following file extensions ",paste(fileExtensions,sep=" "), call.=FALSE);} - -# if set to some x-position, this will trigger the plotting of a vertical line at that location -xVline = NA; -#xVline = 250; # DEBUG TEST VALUE - - -#################################################################### -# parse IntaRNA output -#################################################################### - -d = read.csv2( intarnaOutputFile ) -# check if all columns present -for( x in c(id,start,end)) { - if (!is.element(x, colnames(d))) { - stop("'",id,"' is not among the column names of '",intarnaOutputFile,"'", call.=FALSE); - } -} - -#################################################################### -# create count plot -#################################################################### - -allPos = c(); -for( i in 1:nrow(d) ) { - allPos = c( allPos, d[i,start]:d[i,end] ); -} -allPos = as.data.frame(allPos,ncol=1) -#allPos # DEBUG OUT - -coveragePlot = - ggplot( allPos, aes(x=allPos, stat(count))) + - geom_density() + - ylab("coverage") + - xlab("position") + - scale_y_continuous(position = "right", expand=expand_scale(mult = c(0, .02))) + - scale_x_continuous(expand = c(0, 0)) + - theme( axis.title.x=element_blank() ) - -if ( ! is.na(xVline)) { - coveragePlot = coveragePlot + - geom_vline(aes(xintercept=xVline)); -} - -#################################################################### -# create region plot -#################################################################### - -dRegion = data.frame() -dRegion[1:nrow(d),1:3] = d[,c(start,end,id)] -dRegion[1:nrow(d),4] = factor(sprintf("%08d",nrow(d):1)) -colnames(dRegion) = c("start","end","id","idx"); -#dRegion # DEBUG OUT - -yLabelScale = 0.6 # if you have to alter for more/less sequence IDs per inch; see below - -regionPlot = - ggplot(dRegion, aes(x=start,xend=end,y=idx)) + - geom_dumbbell(color="dodgerblue", size=2) + - xlab("position") + - scale_x_continuous( expand = c(0, 0) ) + - ylab("") + - scale_y_discrete(position = "right", breaks=dRegion$idx, labels=dRegion$id) + - geom_vline(aes(xintercept=min(allPos))) + - theme(panel.grid.major.y=element_line(size=0.7,color="lightgray") - , axis.text.y=element_text(size=rel(yLabelScale)) - ) - -if ( ! is.na(xVline)) { - regionPlot = regionPlot + - geom_vline(aes(xintercept=xVline)); -} - -#################################################################### -# plot to file -#################################################################### - -plotWidth = 6 -plotHeightDensity = 2 -plotHeight = plotHeightDensity + max(2,nrow(d)/9) -plotHeightDensityRel = plotHeightDensity / plotHeight - - -plot_grid( coveragePlot, regionPlot - , nrow=2, ncol=1 - , align = "hv" - , axis = "r" - , rel_heights= c( plotHeightDensityRel, 1.0-plotHeightDensityRel) -) - -ggsave( outFile - , width= plotWidth - , height= plotHeight -); - -#############################################################EOF +#!/usr/bin/env Rscript + +#################################################################### +# Visualization of sequence regions covered by RNA-RNA interactions +# predicted by IntaRNA. +# +# arguments: <1|2> +# +# 1 = ";"-separated CSV output of IntaRNA +# 2 <1|2|paramFile> = suffix of "start,end,id" CSV cols to plot +# or paramFile containing '='-separated assignments for +# id, start, end, title, xvline, xmin, xmax +# 3 = file name of the output figure suffixed +# by one of ".pdf",".png",".svg",".eps",".ps",".jpeg",".tiff" +# +# example call: +# +# Rscript --vanilla plotRegions.R pred.csv 1 regions.png +# +# This script is part of the IntaRNA source code package. See +# respective licence and documentation for further information. +# +# https://github.com/BackofenLab/IntaRNA +# +#################################################################### + + +#################################################################### +# check and load dependencies +#################################################################### + +options(warn=-1) +suppressPackageStartupMessages(library(ggplot2)) +suppressPackageStartupMessages(library(ggalt)) +suppressPackageStartupMessages(library(cowplot)) # cowplot starts with a note +options(warn=0) + +theme_set(theme_cowplot()) + +#################################################################### +# get command line arguments +#################################################################### + +args = commandArgs(trailingOnly=TRUE) +# check and parse +if (length(args)!=3) { stop("call with <1|2|paramFile> ", call.=FALSE) } + +intarnaOutputFile = args[1]; +if (!file.exists(intarnaOutputFile )) { stop("intarna-csv-output file '", intarnaOutputFile, "' does not exist!", call.=FALSE) } + +# init columns to plot +id = NA +start = NA +end = NA +title = NULL +# if set to some x-position, this will trigger the plotting of a vertical line at that location +xVline = NA; +xmin = NA; +xmax = NA; +rowsMax = 200; +# set columns to plot +seqNr = args[2]; +if (seqNr == "1" || seqNr == "2") { + id = paste("id",seqNr,sep=""); + start = paste("start",seqNr,sep=""); + end = paste("end",seqNr,sep=""); +} else { + if (!file.exists(seqNr )) { stop("second call argument as to be '1' or '2' to specify which regions to plot or the name of the parameter file to parse"); } + p = read.table(seqNr, header=FALSE, row.names = 1, sep="=", quote = "", strip.white=TRUE, blank.lines.skip=TRUE, comment.char="#") + # set parsed data + id = as.character(p["id",1]) + start = as.character(p["start",1]) + end = as.character(p["end",1]) + title = as.character(p["title",1]) + xVline = as.numeric(as.character(p["xvline",1])) + xmin = as.numeric(as.character(p["xmin",1])) + xmax = as.numeric(as.character(p["xmax",1])) + rowsMax = as.numeric(as.character(p["rows",1])) +} + +outFile = args[3]; +fileExtensions = c(".pdf",".png",".svg",".eps",".ps",".jpeg",".tiff"); +outFileExtOk = FALSE; +for( ext in fileExtensions ) { outFileExtOk = outFileExtOk || endsWith(outFile,ext); } +if ( !outFileExtOk ) {stop(" has to have one of the following file extensions ",paste(fileExtensions,sep=" "), call.=FALSE);} + +#################################################################### +# parse IntaRNA output +#################################################################### + +d = read.csv2( intarnaOutputFile ) +# check if all columns present +for( x in c(id,start,end)) { + if (!is.element(x, colnames(d))) { + stop("'",id,"' is not among the column names of '",intarnaOutputFile,"'", call.=FALSE); + } +} + +# reduce to rows of interest +rowsMax = min(rowsMax, nrow(d)) +d = d[1:rowsMax,] + +#################################################################### +# create count plot +#################################################################### + +allPos = c(); +for( i in 1:nrow(d) ) { + allPos = c( allPos, d[i,start]:d[i,end] ); +} +allPos = as.data.frame(allPos,ncol=1) +#allPos # DEBUG OUT + +if ( is.na(xmin) ) { + xmin = min(allPos) +} +if ( is.na(xmax) ) { + xmax = max(allPos) +} + +coveragePlot = + ggplot( allPos, aes(x=allPos, stat(count))) + + geom_density() + + ylab("coverage") + + xlab( ifelse( is.null(title) , "position" , title ) ) + + scale_y_continuous(position = "right", expand=expand_scale(mult = c(0, .02))) + + scale_x_continuous(expand = c(0, 0), limits=c(xmin,xmax)); + +if (!is.null(title)) { + coveragePlot = coveragePlot + theme( axis.title.x=element_text(hjust = 0.5, face = "bold")); +} else { + coveragePlot = coveragePlot + theme( axis.title.x=element_blank()); +} + +# plot vertical line if requested +if ( ! is.na(xVline) && xVline >= xmin && xVline <= xmax ) { + coveragePlot = coveragePlot + + geom_vline(aes(xintercept=xVline)); +} + +#################################################################### +# create region plot +#################################################################### + +dRegion = data.frame() +dRegion[1:nrow(d),1:3] = d[,c(start,end,id)] +dRegion[1:nrow(d),4] = factor(sprintf("%08d",nrow(d):1)) +colnames(dRegion) = c("start","end","id","idx"); +#dRegion # DEBUG OUT + +yLabelScale = 0.6 # if you have to alter for more/less sequence IDs per inch; see below + +regionPlot = + ggplot(dRegion, aes(x=start,xend=end,y=idx)) + + geom_dumbbell(color="dodgerblue", size=2) + + xlab("position") + + scale_x_continuous( expand = c(0, 0), limits=c(xmin,xmax) ) + + ylab("") + + scale_y_discrete(position = "right", breaks=dRegion$idx, labels=dRegion$id) + + geom_vline(aes(xintercept=(xmin))) + + theme(panel.grid.major.y=element_line(size=0.7,color="lightgray") + , axis.text.y=element_text(size=rel(yLabelScale)) + #, plot.title = element_blank() + ) + +if ( ! is.na(xVline)) { + regionPlot = regionPlot + + geom_vline(aes(xintercept=xVline)); +} + +#################################################################### +# plot to file +#################################################################### + +plotWidth = 6 +plotHeightDensity = 2 +plotHeight = plotHeightDensity + max(2,nrow(d)/9) +plotHeightDensityRel = plotHeightDensity / plotHeight + +options(warn=-1) # disable warnings of unneeded rows +plot_grid( coveragePlot, regionPlot + , nrow=2, ncol=1 + , align = "hv" + , axis = "r" + , rel_heights= c( plotHeightDensityRel, 1.0-plotHeightDensityRel) +) +options(warn=0) # reenable warnings + +ggsave( outFile + , width= plotWidth + , height= plotHeight +); + +#############################################################EOF diff --git a/R/Makefile.am b/R/Makefile.am new file mode 100644 index 00000000..f0562006 --- /dev/null +++ b/R/Makefile.am @@ -0,0 +1,8 @@ + +# scripts to be installed in $(bindir) and distributed + +dist_bin_SCRIPTS = \ + IntaRNA_CSV_p-value.R \ + IntaRNA_plotRegions.R + + diff --git a/R/README.md b/R/README.md index 0690fa97..3d718cfb 100644 --- a/R/README.md +++ b/R/README.md @@ -1,12 +1,42 @@ # Auxiliary R scripts of the IntaRNA package +- [IntaRNA_CSV_p-value.R](#intarnacsvpvaluer) +- [IntaRNA_plotRegions.R](#intrarnaplotregionsr) -# `plotRegions.R` - Visualization of RRI-covered regions +# `IntaRNA_CSV_p-value.R` +### p-value estimates of IntaRNA energies from present energy distribution + +For genome-wide target sRNA predictions, we assume the set of energy values +predicted for all target sufficiently diverse to be used as a background energy +model for minimum energies of putative target sequences. Thus, we can fit a +generalized extreme value (GEV) distribution to the data that is subsequently +used to estimate p-values for each energy. + +The `IntaRNA_CSV_p-value.R` script takes an IntaRNA CSV output file (assumed to be +sufficiently large and sane enough for GEV fitting) to compute respective +p-value estimates. The output consists of the input table extended with a +`p-value` column. If no output file is given or in- and output file names are +equal, the input file is overwritten! +You can (optionally) specify the column name for which p-values are to be +estimated. Example calls are given below. + +```bash +# overwriting the input file with p-value-extended table +Rscript --vanilla IntaRNA_CSV_p-value.R IntaRNA-output.csv +# creating a new output file for p-value extension +Rscript --vanilla IntaRNA_CSV_p-value.R IntaRNA-output.csv IntaRNA-output-with-pValue.csv +# computing p-values for normalized energies (has to be present in file IN.csv) +Rscript --vanilla IntaRNA_CSV_p-value.R IN.csv IN-pValue.csv E_norm +``` + + +# `IntaRNA_plotRegions.R` +### Visualization of RRI-covered regions To visualize sequences' regions covered by RNA-RNA interactions predicted by -IntaRNA, you can use `plotRegions.R` by providing the following arguments (in +IntaRNA, you can use `IntaRNA_plotRegions.R` by providing the following arguments (in the given order) 1. CSV-IntaRNA output file (semicolon separated) covering the columns `start,end,id` @@ -17,7 +47,7 @@ the given order) An example is given below, when calling ```bash -Rscript --vanilla plotRegions.R pred.csv 1 regions.png +Rscript --vanilla IntaRNA_plotRegions.R pred.csv 1 plotRegions.example.png ``` with `pred.csv` containing @@ -31,5 +61,5 @@ b0005;281;295;query;5;22 ``` will produce the output -![plotRegions.R example](plotRegions.example.png) +![IntaRNA_plotRegions.R example](plotRegions.example.png) diff --git a/configure.ac b/configure.ac index 7d2cbd24..41871bab 100644 --- a/configure.ac +++ b/configure.ac @@ -1,7 +1,7 @@ AC_PREREQ([2.65]) # 5 argument version only available with aclocal >= 2.64 -AC_INIT([IntaRNA], [3.1.3], [], [intaRNA], [http://www.bioinf.uni-freiburg.de] ) +AC_INIT([IntaRNA], [3.1.4], [], [intaRNA], [http://www.bioinf.uni-freiburg.de] ) # minimal required version of the boost library BOOST_REQUIRED_VERSION=1.50.0 @@ -59,7 +59,7 @@ AX_CXX_COMPILE_STDCXX([11], [noext], [mandatory]) # check if python is available AM_PATH_PYTHON([$PYTHON_REQUIRED_VERSION],, [:]) -AM_CONDITIONAL([HAVE_PYTHON], [test "$PYTHON" != :]) +AM_CONDITIONAL([HAVE_PYTHON], [test "$PYTHON" != ":"]) AM_COND_IF([HAVE_PYTHON], [HAVE_PYTHON_yesno=yes], [HAVE_PYTHON_yesno=no]) ############################################################################### @@ -291,11 +291,13 @@ AM_CONDITIONAL([INSTALL_INTARNAPVALUE],[test x"$intarnapvalue" = x"yes"]) ############################################################################### AX_BOOST_BASE([$BOOST_REQUIRED_VERSION], [FOUND_BOOST=1;], [FOUND_BOOST=0;]) -AX_BOOST_SYSTEM -AX_BOOST_FILESYSTEM -AX_BOOST_PROGRAM_OPTIONS -AX_BOOST_REGEX -AX_BOOST_IOSTREAMS +AS_IF([test "x$want_boost" != "xno" && test "x$FOUND_BOOST" = "x1"], [ + AX_BOOST_SYSTEM + AX_BOOST_FILESYSTEM + AX_BOOST_PROGRAM_OPTIONS + AX_BOOST_REGEX + AX_BOOST_IOSTREAMS +]) AX_CHECK_ZLIB([], []) @@ -313,7 +315,7 @@ AC_HEADER_STDC ########################################################################## # FOUND_BOOST is only defined if want_boost is "yes" -AS_IF([test $want_boost = "no" || test $FOUND_BOOST != 1], [ +AS_IF([test $want_boost != "no" && test $FOUND_BOOST != 1], [ AC_MSG_NOTICE([]) AC_MSG_NOTICE([The Boost Library was not found!]) AC_MSG_NOTICE([ -> If installed in a non-standard path, please use '--with-boost=PREFIX'.]) @@ -430,6 +432,7 @@ AC_CONFIG_FILES([src/IntaRNA/Makefile]) AC_CONFIG_FILES([src/IntaRNA/intarna_config.h]) AC_CONFIG_FILES([src/bin/Makefile]) AC_CONFIG_FILES([perl/Makefile]) +AC_CONFIG_FILES([R/Makefile]) AC_CONFIG_FILES([python/Makefile]) AC_CONFIG_FILES([python/bin/Makefile]) AC_CONFIG_FILES([python/intarnapvalue/Makefile]) diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp index 4e124e0f..58c2035f 100644 --- a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp @@ -73,8 +73,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) ) + , 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) ) { E_type seedE = seedHandler.getSeedE(si1, si2); const size_t sl1 = seedHandler.getSeedLength1(si1, si2); diff --git a/src/IntaRNA/PredictorMfe2dSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dSeedExtension.cpp index b98278ee..a178b895 100644 --- a/src/IntaRNA/PredictorMfe2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfe2dSeedExtension.cpp @@ -74,8 +74,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) ) + , 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) ) { // get energy and boundaries of seed E_type seedE = seedHandler.getSeedE(si1, si2); @@ -327,8 +327,8 @@ traceBack( Interaction & interaction ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , i1,j1+1-seedHandler.getConstraint().getBasePairs() - , i2,j2+1-seedHandler.getConstraint().getBasePairs() ) ) + , i1,j1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , i2,j2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2() ) ) { E_type seedE = seedHandler.getSeedE(si1, si2); const size_t sl1 = seedHandler.getSeedLength1(si1, si2); diff --git a/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp b/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp index 6e933154..864e652a 100644 --- a/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp +++ b/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp @@ -78,8 +78,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , 0, interaction_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, interaction_size2+1-seedHandler.getConstraint().getBasePairs()) ) + , 0, interaction_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , 0, interaction_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) ) { E_type seedE = seedHandler.getSeedE(si1, si2); size_t sl1 = seedHandler.getSeedLength1(si1, si2); @@ -373,8 +373,8 @@ traceBack( Interaction & interaction ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , i1,j1+1-seedHandler.getConstraint().getBasePairs() - , i2,j2+1-seedHandler.getConstraint().getBasePairs() ) ) + , i1,j1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , i2,j2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2() ) ) { E_type seedE = seedHandler.getSeedE(si1, si2); diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp index bc1fe839..21523ea6 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp @@ -75,8 +75,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) ) + , 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) ) { const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) ); const size_t sl1 = seedHandler.getSeedLength1(si1, si2); diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index cc7aa459..4e5da9bb 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -76,8 +76,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) ) + , 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1() + , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) ) { // get Z and boundaries of seed const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) ); diff --git a/src/IntaRNA/SeedConstraint.h b/src/IntaRNA/SeedConstraint.h index 04532809..1d7bd4da 100644 --- a/src/IntaRNA/SeedConstraint.h +++ b/src/IntaRNA/SeedConstraint.h @@ -39,6 +39,7 @@ class SeedConstraint { * @param explicitSeeds the encodings of explicit seed interactions to be used * @param noGUallowed whether or not GU base pairs are allowed within seeds * @param noGUendAllowed whether or not GU base pairs are allowed at the ends of seeds + * @param noLP whether or not lonely base pairs are allowed */ SeedConstraint( const size_t bp , const size_t maxUnpairedOverall @@ -52,6 +53,7 @@ class SeedConstraint { , const std::string & explicitSeeds , const bool noGUallowed , const bool noGUendAllowed + , const bool noLP ); virtual ~SeedConstraint(); @@ -148,6 +150,13 @@ class SeedConstraint { bool isGUendAllowed() const; + /** + * Whether or not lonely base pairs are allowed + * @return true if lonely base pairs are allowed; false otherwise + */ + bool + isLpAllowed() const; + /** * Index ranges in seq1 to be searched for seeds or empty if all indices * are to be considered. @@ -240,6 +249,9 @@ class SeedConstraint { //! whether or not GU base pairs are allowed at seed ends bool bpGUendAllowed; + //! whether or not lonely base pairs are allowed + bool lpAllowed; + }; @@ -261,6 +273,7 @@ SeedConstraint::SeedConstraint( , const std::string & explicitSeeds , const bool noGUallowed , const bool noGUendAllowed + , const bool noLP ) : bp(bp_) @@ -275,6 +288,7 @@ SeedConstraint::SeedConstraint( , explicitSeeds(explicitSeeds) , bpGUallowed(!noGUallowed) , bpGUendAllowed(!noGUendAllowed) + , lpAllowed(!noLP) { if (bp < 2) throw std::runtime_error("SeedConstraint() : base pair number ("+toString(bp)+") < 2"); } @@ -440,6 +454,15 @@ isGUendAllowed() const { ///////////////////////////////////////////////////////////////////////////// +inline +bool +SeedConstraint:: +isLpAllowed() const { + return lpAllowed; +} + +///////////////////////////////////////////////////////////////////////////// + inline std::ostream& operator<<(std::ostream& out, const SeedConstraint& c) diff --git a/src/IntaRNA/SeedHandlerMfe.cpp b/src/IntaRNA/SeedHandlerMfe.cpp index a0308161..5c276ecd 100644 --- a/src/IntaRNA/SeedHandlerMfe.cpp +++ b/src/IntaRNA/SeedHandlerMfe.cpp @@ -37,6 +37,11 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size size_t i1, i2, bpIn, u1, u2, j1, j2, u1p, u2p, k1,k2, u1best, u2best; E_type curE, bestE; + // determine whether or not lonely base pairs are allowed or if we have to + // ensure a stacking to the right of the left boundary (i1,i2) + const size_t noLpShift = seedConstraint.isLpAllowed() ? 0 : 1; + E_type iStackE = E_type(0); + size_t seedCountNotInf = 0, seedCount = 0; // fill for all start indices @@ -50,8 +55,8 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size // init according to no seed interaction seed(i1-offset1,i2-offset2) = SeedMatrix::value_type( E_INF, 0 ); - // skip non-complementary or infeasible left seed boundaries - if (!isFeasibleSeedBasePair(i1,i2,true)) { + // skip non-complementary or infeasible seed base pair + if (!isFeasibleSeedBasePair(i1,i2)) { continue; // go to next seedE index } @@ -75,6 +80,17 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size // check if this index range is to be considered for seed search bool validSeedSite = isFeasibleSeedBasePair(j1,j2,true); + // if no LP : check for direct right-stacking of i + if (noLpShift > 0) { + // check if feasible extension + if (isFeasibleSeedBasePair(i1+1,i2+1)) { + // get stacking energy + iStackE = energy.getE_interLeft(i1,i1+1,i2,i2+1); + } else { + validSeedSite = false; + } + } + // init current seed energy curE = E_INF; @@ -83,31 +99,47 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size // base case: only left and right base pair present if (bpIn==0) { - // energy for stacking/bulge/interior depending on u1/u2 - curE = energy.getE_interLeft(i1,j1,i2,j2); + // if lonely bps are allowed or no bulge + if (noLpShift == 0 || (u1==0 && u2==0)) { + // energy for stacking/bulge/interior depending on u1/u2 + curE = energy.getE_interLeft(i1,j1,i2,j2); + } } else { - // split seed recursively into all possible leading interior loops - // i1 .. i1+u1p+1 .. j1 - // i2 .. i2+u2p+1 .. j2 - for (u1p=1+std::min(u1,energy.getMaxInternalLoopSize1()); u1p-- > 0;) { - for (u2p=1+std::min(u2,energy.getMaxInternalLoopSize2()); u2p-- > 0;) { - - k1 = i1+u1p+1; - k2 = i2+u2p+1; - // check if split pair is complementary - // and recursed entry is < E_INF - if (! (isFeasibleSeedBasePair(k1,k2) && E_isNotINF( getSeedE( k1-offset1, k2-offset2, bpIn-1, u1-u1p, u2-u2p ) ) ) ) { - continue; // not complementary -> skip - } - - // update mfe for split at k1,k2 - curE = std::min( curE, - energy.getE_interLeft(i1,k1,i2,k2) - + getSeedE( k1-offset1, k2-offset2, bpIn-1, u1-u1p, u2-u2p ) - ); - } // u2p - } // u1p + + // explicitly check direct stacking extension in noLP mode + if (noLpShift > 0 && E_isNotINF( getSeedE( i1+1-offset1, i2+1-offset2, bpIn-1, u1, u2 ) )) { + curE = std::min( curE, iStackE + getSeedE( i1+1-offset1, i2+1-offset2, bpIn-1, u1, u2 ) ); + } + + // if enough interior base pairs left + if (bpIn >= 1+noLpShift) { + // split seed recursively into all possible leading interior loops + // i1 .. i1+u1p+1 .. j1 + // i2 .. i2+u2p+1 .. j2 + for (u1p=1+std::min(u1,energy.getMaxInternalLoopSize1()); u1p-- > 0;) { + for (u2p=1+std::min(u2,energy.getMaxInternalLoopSize2()); u2p-- > 0;) { + + // skip stacked extension for noLP since already covered above + if (u1p+u2p < noLpShift) { continue; } + + k1 = i1+u1p+1+noLpShift; + k2 = i2+u2p+1+noLpShift; + // check if split pair is complementary + // and recursed entry is < E_INF + if (! (isFeasibleSeedBasePair(k1,k2) && E_isNotINF( getSeedE( k1-offset1, k2-offset2, bpIn-1-noLpShift, u1-u1p, u2-u2p ) ) ) ) { + continue; // not complementary -> skip + } + + // update mfe for split at k1,k2 + curE = std::min( curE, + iStackE + + energy.getE_interLeft(i1+noLpShift,k1,i2+noLpShift,k2) + + getSeedE( k1-offset1, k2-offset2, bpIn-1-noLpShift, u1-u1p, u2-u2p ) + ); + } // u2p + } // u1p + } } // more than two base pairs } // (j1,j2) complementary @@ -119,9 +151,10 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size } // u1 // check if full base pair number reached - if (bpIn+1==seedE_rec.shape()[2]) { + // check if feasible left seed boundaries + if ( bpIn+1==seedE_rec.shape()[2] && isFeasibleSeedBasePair(i1,i2,true) ) { - // find best unpaired combination in seed seed for i1,i2,bp + // find best unpaired combination in seed for i1,i2,bp u1best = 0; u2best = 0; bestE = E_INF; @@ -137,16 +170,21 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size // skip if ED boundary exceeded if (energy.getED1(i1,j1) >= seedConstraint.getMaxED() - || energy.getED2(i2,j2) >= seedConstraint.getMaxED() ) + || energy.getED2(i2,j2) >= seedConstraint.getMaxED() ) { continue; } // get overall interaction energy - curE = energy.getE( i1, j1, i2, j2, getSeedE( i1-offset1, i2-offset2, bpIn, u1, u2 ) ) + energy.getE_init(); + E_type curEhyb = getSeedE( i1-offset1, i2-offset2, bpIn, u1, u2 ) + energy.getE_init(); + curE = energy.getE( i1, j1, i2, j2, curEhyb ); + // check hybrid energy bound including Einit // check if better than what is known so far - if ( curE < bestE ) { + if ( curEhyb <= seedConstraint.getMaxEhybrid() + && curE <= seedConstraint.getMaxE() + && curE < bestE ) + { bestE = curE; u1best = u1; u2best = u2; @@ -156,18 +194,10 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size // reduce bestE to hybridization energy only (init+loops) if (E_isNotINF( bestE )) { - // overwrite all seeds with too high energy -> infeasible start interactions - if (bestE >= seedConstraint.getMaxE() - // check hybridization energy bound (incl E_init) - || (getSeedE( i1-offset1, i2-offset2, bpIn, u1best, u2best ) + energy.getE_init()) >= seedConstraint.getMaxEhybrid()) - { - bestE = E_INF; - } else { - // get seed's hybridization loop energies only - bestE = getSeedE( i1-offset1, i2-offset2, bpIn, u1best, u2best ); - // count true seed - seedCountNotInf++; - } + // get seed's hybridization loop energies only + bestE = getSeedE( i1-offset1, i2-offset2, bpIn, u1best, u2best ); + // count true seed + seedCountNotInf++; } // store best (mfe) seed for all u1/u2 @@ -217,10 +247,17 @@ traceBackSeed( Interaction & interaction // get energy of provided seed E_type curE = getSeedE(i1_,i2_,bpInbetween,u1_,u2_); + // determine whether or not lonely base pairs are allowed or if we have to + // ensure a stacking to the right of the left boundary (i1,i2) + const size_t noLpShift = seedConstraint.isLpAllowed() ? 0 : 1; + E_type iStackE = E_type(0); + // trace seed // trace each seed base pair (excluding right most) for( size_t bpIn=1+bpInbetween; bpIn-- > 0; ) { + + // base case: only left and right base pair present if (bpIn==0) { // add left base pair if not left seed boundary @@ -229,6 +266,24 @@ traceBackSeed( Interaction & interaction } } else { + + // if no LP : check for direct right-stacking of i + if (noLpShift > 0) { + // check if feasible extension + assert(isFeasibleSeedBasePair(i1+1,i2+1)); + // get stacking energy + iStackE = energy.getE_interLeft(i1+offset1,i1+1+offset1,i2+offset2,i2+1+offset2); + // noLP : check stacking of i + if ( E_equal( curE, iStackE + getSeedE( i1+1, i2+1, bpIn-1, u1max, u2max )) ) { + i1++; + i2++; + curE = getSeedE( i1+1, i2+1, bpIn-1, u1max, u2max ); + continue; + } + // sanity check for noLP mode + assert( bpIn >= 1+noLpShift ); + } + // split seed recursively into all possible leading interior loops // i1 .. i1+u1p+1 .. j1 // i2 .. i2+u2p+1 .. j2 @@ -236,26 +291,33 @@ traceBackSeed( Interaction & interaction for (u1=1+u1max; traceNotFound && u1-- > 0;) { for (u2=1+u2max; traceNotFound && u2-- > 0;) { // check if overall number of unpaired is not exceeded - if (u1+u2 > uMax) { + // or skip stacked extension since covered above + if (u1+u2 > uMax || u1+u2 < noLpShift) { continue; } - k1 = i1+u1+1; - k2 = i2+u2+1; + k1 = i1+u1+1+noLpShift; + k2 = i2+u2+1+noLpShift; // check if valid trace - if ( E_isNotINF( getSeedE( k1, k2, bpIn-1, u1max-u1, u2max-u2 ) ) ) { + if ( isFeasibleSeedBasePair(k1+offset1, k2+offset2) && E_isNotINF( getSeedE( k1, k2, bpIn-1, u1max-u1, u2max-u2 ) ) ) { // check if correct trace - if ( E_equal( curE, energy.getE_interLeft(i1+offset1,k1+offset1,i2+offset2,k2+offset2) - + getSeedE( k1, k2, bpIn-1, u1max-u1, u2max-u2 )) ) + if ( E_equal( curE, iStackE + + energy.getE_interLeft(i1+noLpShift+offset1,k1+offset1,i2+noLpShift+offset2,k2+offset2) + + getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 )) ) { // store left base pair if not left seed boundary if (i1 != i1_) { interaction.basePairs.push_back( energy.getBasePair(i1+offset1,i2+offset2) ); } + if (noLpShift > 0) { + interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift+offset1,i2+noLpShift+offset2) ); + // reflect additional base pair + bpIn--; + } // store next energy value to trace - curE = getSeedE( k1, k2, bpIn-1, u1max-u1, u2max-u2 ); + curE = getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 ); // reset for next trace step i1 = k1; i2 = k2; diff --git a/src/IntaRNA/SeedHandlerMfe.h b/src/IntaRNA/SeedHandlerMfe.h index 488b76d0..0f770176 100644 --- a/src/IntaRNA/SeedHandlerMfe.h +++ b/src/IntaRNA/SeedHandlerMfe.h @@ -334,7 +334,6 @@ E_type SeedHandlerMfe:: getSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2 ) const { -// return seedE_rec[i1][i2][bpInbetween][u1][u2]; return seedE_rec( SeedIndex({{ (SeedRecMatrix::index) i1 , (SeedRecMatrix::index) i2 @@ -350,7 +349,6 @@ void SeedHandlerMfe:: setSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2, const E_type E ) { -// seedE_rec[i1][i2][bpInbetween][u1][u2] = E; seedE_rec( SeedIndex({{ (SeedRecMatrix::index) i1 , (SeedRecMatrix::index) i2 diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index e853d06f..9d56491f 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -384,7 +384,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ->default_value(qIdxPos0.def) ->notifier(boost::bind(&CommandLineParsing::validate_numberArgument,this,qIdxPos0,_1)) , std::string("index of first (5') sequence position of all query sequences" - " (arg in range ["+toString(qAccW.min)+","+toString(qAccW.max)+"];").c_str()) + " (arg in range ["+toString(qIdxPos0.min)+","+toString(qIdxPos0.max)+"];").c_str()) ("qSet" , value(&(qSetString)) ->notifier(boost::bind(&CommandLineParsing::validate_qSet,this,_1)) @@ -485,7 +485,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ->default_value(tIdxPos0.def) ->notifier(boost::bind(&CommandLineParsing::validate_numberArgument,this,tIdxPos0,_1)) , std::string("index of first (5') sequence position of all target sequences" - " (arg in range ["+toString(qAccW.min)+","+toString(qAccW.max)+"];").c_str()) + " (arg in range ["+toString(tIdxPos0.min)+","+toString(tIdxPos0.max)+"];").c_str()) ("tSet" , value(&(tSetString)) ->notifier(boost::bind(&CommandLineParsing::validate_tSet,this,_1)) @@ -2210,7 +2210,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'B': { switch ( mode.val ) { case 'H' : return new PredictorMfe2dHelixBlockHeuristic( energy, output, predTracker, getHelixConstraint(energy)); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // single-site mfe interactions (contain only interior loops) @@ -2218,7 +2218,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const switch ( mode.val ) { case 'H' : return new PredictorMfe2dHeuristic( energy, output, predTracker ); case 'M' : return new PredictorMfe2d( energy, output, predTracker ); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // single-site mfe ensemble interactions (contain only interior loops) @@ -2226,16 +2226,16 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const switch ( mode.val ) { case 'H' : return new PredictorMfeEns2dHeuristic( energy, output, predTracker ); case 'M' : return new PredictorMfeEns2d( energy, output, predTracker ); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // multi-site mfe interactions (contain interior and multi-loops loops) case 'M' : { switch ( mode.val ) { - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); + default : INTARNA_NOT_IMPLEMENTED("no model "+toString(model.val)+" available for mode "+toString(mode.val)); } } else { // seed-constrained predictors @@ -2243,7 +2243,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'B' : { switch ( mode.val ) { case 'H' : return new PredictorMfe2dHelixBlockHeuristicSeed(energy, output, predTracker, getHelixConstraint(energy), getSeedHandler(energy)); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // single-site mfe interactions (contain only interior loops) @@ -2252,7 +2252,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'H' : return new PredictorMfe2dHeuristicSeed( energy, output, predTracker, getSeedHandler( energy ) ); case 'M' : return new PredictorMfe2dSeed( energy, output, predTracker, getSeedHandler( energy ) ); case 'S' : return new PredictorMfeSeedOnly( energy, output, predTracker, getSeedHandler( energy ) ); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // single-site min energy interactions via seed extension (contain only interior loops) @@ -2271,16 +2271,16 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'H' : return new PredictorMfeEns2dHeuristicSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'M' : return new PredictorMfeEns2dSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'S' : return new PredictorMfeEnsSeedOnly( energy, output, predTracker, getSeedHandler(energy) ); - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; // multi-site mfe interactions (contain interior and multi-loops loops) case 'M' : { switch ( mode.val ) { - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented for model "+toString(model.val)); + default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; - default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); + default : INTARNA_NOT_IMPLEMENTED("no model "+toString(model.val)+" available for mode "+toString(mode.val)); } } } @@ -2382,6 +2382,7 @@ getSeedConstraint( const InteractionEnergy & energy ) const , seedTQ , seedNoGU , seedNoGUend + , outNoLP ); } return *seedConstraint; @@ -2402,9 +2403,6 @@ getSeedHandler( const InteractionEnergy & energy ) const } else { // check if we have to allow for bulges in seed if (seedConstr.getMaxUnpaired1()+seedConstr.getMaxUnpaired2()+seedConstr.getMaxUnpairedOverall() > 0) { - if (outNoLP) { - INTARNA_NOT_IMPLEMENTED("outNoLP not yet implemented for seeds with bulges"); - } // create new seed handler using mfe computation return new SeedHandlerMfe( energy, seedConstr ); } else { diff --git a/tests/HelixHandlerNoBulgeMaxSeedIdxOffset_test.cpp b/tests/HelixHandlerNoBulgeMaxSeedIdxOffset_test.cpp index a621a41a..bb0222c5 100644 --- a/tests/HelixHandlerNoBulgeMaxSeedIdxOffset_test.cpp +++ b/tests/HelixHandlerNoBulgeMaxSeedIdxOffset_test.cpp @@ -30,7 +30,7 @@ TEST_CASE( "HelixHandlerIdxOffset for NoBulgeMaxSeed", "[HelixHandlerIdxOffset]" // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); HelixHandler *hhS = new HelixHandlerNoBulgeMax(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -128,7 +128,7 @@ TEST_CASE( "HelixHandlerIdxOffset for NoBulgeMaxSeed", "[HelixHandlerIdxOffset]" // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); HelixHandler *hhS = new HelixHandlerNoBulgeMax(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -192,7 +192,7 @@ TEST_CASE( "HelixHandlerIdxOffset for NoBulgeMaxSeed", "[HelixHandlerIdxOffset]" // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); HelixHandler *hhS = new HelixHandlerNoBulgeMax(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -268,7 +268,7 @@ TEST_CASE( "HelixHandlerIdxOffset for NoBulgeMaxSeed", "[HelixHandlerIdxOffset]" // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 2, 1, 1, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); HelixHandler *hhS = new HelixHandlerNoBulgeMax(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -364,7 +364,7 @@ TEST_CASE( "HelixHandlerIdxOffset for NoBulgeMaxSeed", "[HelixHandlerIdxOffset]" // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); HelixHandler *hhS = new HelixHandlerNoBulgeMax(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); diff --git a/tests/HelixHandlerNoBulgeMaxSeed_test.cpp b/tests/HelixHandlerNoBulgeMaxSeed_test.cpp index b2be598d..b2712bf5 100644 --- a/tests/HelixHandlerNoBulgeMaxSeed_test.cpp +++ b/tests/HelixHandlerNoBulgeMaxSeed_test.cpp @@ -30,7 +30,7 @@ TEST_CASE( "HelixHandlerNoBulgeMaxSeed", "[HelixHandlerNoBulgeMax]" ) { , AccessibilityDisabled::ED_UPPER_BOUND, 0 , IndexRangeList("") , IndexRangeList("") - , "", false, false ); + , "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerNoBulgeMax hhS(energy, hC); @@ -118,7 +118,7 @@ TEST_CASE( "HelixHandlerNoBulgeMaxSeed", "[HelixHandlerNoBulgeMax]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); SeedHandlerMfe sH(energy, sC); HelixHandlerNoBulgeMax hhS(energy, hC); @@ -172,7 +172,7 @@ TEST_CASE( "HelixHandlerNoBulgeMaxSeed", "[HelixHandlerNoBulgeMax]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); SeedHandlerMfe sH(energy, sC); HelixHandlerNoBulgeMax hhS(energy, hC); @@ -239,7 +239,7 @@ TEST_CASE( "HelixHandlerNoBulgeMaxSeed", "[HelixHandlerNoBulgeMax]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ / seedNoGU SeedConstraint sC(3, 2, 1, 1, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false); + "", false, false, false); SeedHandlerMfe sH(energy, sC); HelixHandlerNoBulgeMax hhS(energy, hC); diff --git a/tests/HelixHandlerUnpairedSeedIdxOffset_test.cpp b/tests/HelixHandlerUnpairedSeedIdxOffset_test.cpp index 8384ba4b..88a1af5e 100644 --- a/tests/HelixHandlerUnpairedSeedIdxOffset_test.cpp +++ b/tests/HelixHandlerUnpairedSeedIdxOffset_test.cpp @@ -27,7 +27,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { HelixConstraint hC(2, 10, 2, 999, 0, false); // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -67,7 +67,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -164,7 +164,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -227,7 +227,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -302,7 +302,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 2, 1, 1, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -397,7 +397,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -494,7 +494,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -582,7 +582,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -685,7 +685,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -751,7 +751,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -833,7 +833,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); @@ -936,7 +936,7 @@ TEST_CASE( "HelixSeed for Unpaired with offset", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 1, 1, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); HelixHandler *hhU = new HelixHandlerUnpaired(energy, hC); SeedHandler *sH = new SeedHandlerMfe(energy, sC); diff --git a/tests/HelixHandlerUnpairedSeed_test.cpp b/tests/HelixHandlerUnpairedSeed_test.cpp index 7966e1b7..9e18ff85 100644 --- a/tests/HelixHandlerUnpairedSeed_test.cpp +++ b/tests/HelixHandlerUnpairedSeed_test.cpp @@ -27,7 +27,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -116,7 +116,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -170,7 +170,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -237,7 +237,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 2, 1, 1, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -324,7 +324,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -404,7 +404,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -499,7 +499,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -557,7 +557,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -632,7 +632,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -728,7 +728,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 1, 1, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); @@ -829,7 +829,7 @@ TEST_CASE( "HelixHandlerUnpairedSeed", "[HelixHandlerUnpaired]" ) { // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); HelixHandlerUnpaired hhU(energy, hC); diff --git a/tests/PredictorMfe2dHelixBlockHeuristicSeed_test.cpp b/tests/PredictorMfe2dHelixBlockHeuristicSeed_test.cpp index cbc7797d..262af548 100644 --- a/tests/PredictorMfe2dHelixBlockHeuristicSeed_test.cpp +++ b/tests/PredictorMfe2dHelixBlockHeuristicSeed_test.cpp @@ -30,7 +30,7 @@ TEST_CASE( "PredictorMfe2dHelixBlockHeuristicSeed", "[PredictorMfe2dHelixBlockHe HelixConstraint hc(2, 4, 0, 999, 0, false); // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); OutputConstraint outC(1, OutputConstraint::OVERLAP_SEQ2, 0, 100); OutputHandlerInteractionList out(outC,1); @@ -67,7 +67,7 @@ TEST_CASE( "PredictorMfe2dHelixBlockHeuristicSeed", "[PredictorMfe2dHelixBlockHe HelixConstraint hc(2, 4, 0, 999, 0, false); // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ SeedConstraint sC(3, 0, 0, 0, 0, AccessibilityDisabled::ED_UPPER_BOUND, 0, IndexRangeList(""), IndexRangeList(""), - "", false, false ); + "", false, false, false ); SeedHandlerMfe sH(energy, sC); OutputConstraint outC(1, OutputConstraint::OVERLAP_SEQ2, 0, 100); diff --git a/tests/SeedHandlerExplicit_test.cpp b/tests/SeedHandlerExplicit_test.cpp index 944845d4..5ce02c4c 100644 --- a/tests/SeedHandlerExplicit_test.cpp +++ b/tests/SeedHandlerExplicit_test.cpp @@ -125,7 +125,7 @@ TEST_CASE( "SeedHandlerExplicit", "[SeedHandlerExplicit]" ) { , IndexRangeList("") , IndexRangeList("") , "1||&3||,2|.|&1||" - , false, false); + , false, false, false); // create instance (triggers parsing) SeedHandlerExplicit sh( energy, sC ); diff --git a/tests/SeedHandlerIdxOffset_test.cpp b/tests/SeedHandlerIdxOffset_test.cpp index 814c9e53..434618eb 100644 --- a/tests/SeedHandlerIdxOffset_test.cpp +++ b/tests/SeedHandlerIdxOffset_test.cpp @@ -31,7 +31,7 @@ TEST_CASE( "SeedHandlerMfe with offset", "[SeedHandlerIdxOffset]") { , IndexRangeList("") , IndexRangeList("") , "" - , false, false ); + , false, false, false ); SeedHandlerNoBulge shOrig(energy, sC); SeedHandlerIdxOffset shIO( &shOrig, false ); diff --git a/tests/SeedHandlerMfe_test.cpp b/tests/SeedHandlerMfe_test.cpp index 6bcee5e7..6f5487cc 100644 --- a/tests/SeedHandlerMfe_test.cpp +++ b/tests/SeedHandlerMfe_test.cpp @@ -29,7 +29,7 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { , IndexRangeList("") , IndexRangeList("") , "" - , false, false + , false, false, false ); SeedHandlerMfe sHM(energy, sC); @@ -39,4 +39,104 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { ////////////////////////////////////////////////////////////////////////////////////////////////////////// REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) > 0); } + + SECTION("SeedHandlerMfe: with lp", "[SeedHandlerMfe]") { + RnaSequence r1("r1", "GGCGCGG"); + RnaSequence r2("r2", "CCCCC"); + AccessibilityDisabled acc1(r1, 0, NULL); + AccessibilityDisabled acc2(r2, 0, NULL); + ReverseAccessibility racc(acc2); + InteractionEnergyBasePair energy(acc1, racc); + + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(4,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, false + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 4); // getting only one seed per si, thus 2x2 + } + + SECTION("SeedHandlerMfe: no lp - 1", "[SeedHandlerMfe]") { + RnaSequence r1("r1", "GCGCGG"); + RnaSequence r2("r2", "CCCCC"); + AccessibilityDisabled acc1(r1, 0, NULL); + AccessibilityDisabled acc2(r2, 0, NULL); + ReverseAccessibility racc(acc2); + InteractionEnergyBasePair energy(acc1, racc); + + { + size_t bp = 2; + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(bp,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, true + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 4); + } + + for (size_t bp = 3; bp <= 4; bp++) { + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(bp,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, true + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 0); + } + } + + SECTION("SeedHandlerMfe: no lp - 2", "[SeedHandlerMfe]") { + RnaSequence r1("r1", "GGCGCGG"); + RnaSequence r2("r2", "CCCCC"); + AccessibilityDisabled acc1(r1, 0, NULL); + AccessibilityDisabled acc2(r2, 0, NULL); + ReverseAccessibility racc(acc2); + InteractionEnergyBasePair energy(acc1, racc); + + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(4,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, true + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 2); + } } diff --git a/tests/SeedHandlerNoBulge_test.cpp b/tests/SeedHandlerNoBulge_test.cpp index e3663514..2e2fa14f 100644 --- a/tests/SeedHandlerNoBulge_test.cpp +++ b/tests/SeedHandlerNoBulge_test.cpp @@ -32,7 +32,7 @@ TEST_CASE( "SeedHandlerNoBulge", "[SeedHandlerNoBulge]" ) { // LOG(DEBUG) <<"r2 : "< 0 ); } { // should find some seed - SeedConstraint sConstr(3,0,0,0,Ekcal_2_E(-2),10,0,IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,Ekcal_2_E(-2),10,0,IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) > 0 ); } { // should find no seed - SeedConstraint sConstr(3,0,0,0,Ekcal_2_E(-3),10,0,IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,Ekcal_2_E(-3),10,0,IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) == 0 ); } @@ -117,19 +117,19 @@ TEST_CASE( "SeedHandlerNoBulge", "[SeedHandlerNoBulge]" ) { { // should find some seeds - SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-1),IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-1),IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) > 0 ); } { // should find some seeds - SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-2),IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-2),IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) > 0 ); } { // should find no seed - SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-3),IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-3),IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) == 0 ); } @@ -144,7 +144,7 @@ TEST_CASE( "SeedHandlerNoBulge", "[SeedHandlerNoBulge]" ) { InteractionEnergyBasePair energy( acc1, rAcc2 ); // should find some seeds - SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-1),IndexRangeList(),IndexRangeList(),"", false, false); + SeedConstraint sConstr(3,0,0,0,0,10,Ekcal_2_E(-1),IndexRangeList(),IndexRangeList(),"", false, false, false); SeedHandlerNoBulge sh(energy,sConstr); REQUIRE( sh.fillSeed(0,energy.size1()-1,0,energy.size2()-1) > 0 );