diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 18668a64..d7f73eab 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -4,7 +4,7 @@ on: push: branches: [master] pull_request: - branches: [master] + branches: [master, feat/294-python-intarna-fuse] defaults: run: @@ -105,6 +105,8 @@ jobs: make -j 2 make install make tests -j 2 + python3 -m pip install -r python/requirements-output-tools.txt + PATH="$RUNNER_TEMP/intarna-install/bin:$PATH" make -C python check-python-tools "$RUNNER_TEMP/intarna-install/bin/IntaRNA" -h - name: Compile every installed public header independently diff --git a/.github/workflows/python-output-tools.yml b/.github/workflows/python-output-tools.yml new file mode 100644 index 00000000..9ee186b1 --- /dev/null +++ b/.github/workflows/python-output-tools.yml @@ -0,0 +1,35 @@ +name: Python companion tools + +on: + push: + branches: [master] + pull_request: + +jobs: + companion-tools: + runs-on: ubuntu-24.04 + strategy: + matrix: + include: + - python: '3.10' + dependencies: numpy==1.23.5 scipy==1.9.3 matplotlib==3.6.3 + - python: '3.13' + dependencies: -r python/requirements-output-tools.txt + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python }} + - run: python -m pip install setuptools wheel ${{ matrix.dependencies }} + - run: python -m unittest discover -s python/tests -v + env: + MPLBACKEND: Agg + - run: python tests/testIntaRNAFuse.py + - name: Check installed commands away from the checkout + run: | + python -m venv --system-site-packages "$RUNNER_TEMP/companions" + "$RUNNER_TEMP/companions/bin/python" -m pip install --no-deps --no-build-isolation ./python + cd "$RUNNER_TEMP" + for tool in CopomuS.py IntaRNA-fuse.py IntaRNA_CSV_p-value.py IntaRNA_plotRegions.py IntaRNApvalue; do + "$RUNNER_TEMP/companions/bin/$tool" --help + done diff --git a/.gitignore b/.gitignore index ee819cc0..aa570ee3 100644 --- a/.gitignore +++ b/.gitignore @@ -82,3 +82,11 @@ python/intarnapvalue.* /IntaRNA.pc /IntaRNA.exe.stackdump /py-compile + +# Python companion-tool build products +__pycache__/ +*.py[cod] +*.egg-info/ +python/build/ +python/dist/ +python/bin/build/ diff --git a/AGENTS.md b/AGENTS.md index fa368888..4a900620 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -11,7 +11,7 @@ with the code and build files when their conventions change. - `src/bin/`: the executable and `CommandLineParsing` for CLI options and wiring. - `tests/`: Catch-based API tests; `tests/data/`: CLI parameters and reference results exercised by `tests/runIntaRNA.sh`. -- `python/`, `R/`: companion tools; consult their local documentation +- `python/`: companion tools; consult their local documentation when changing them. - [README.md](README.md): user-facing installation, options and library usage. [doc/refactor/](doc/refactor/) records earlier analyses and measurements; diff --git a/ChangeLog b/ChangeLog index 35ce0126..18db4a1f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,12 @@ ## Interface and handling +- fix IntaRNApvalue both-strand sampling and parameter consistency, bound retries, + validate FASTA/fits and expose empirical p-values (PR #298) + +- fix CopomuS explicit mutations, GC filters, signed indices, zero scores and + shared profiles; validate inputs and preserve output on failure (PR #298) + - orient probability SVGs from the upper-left corner and highlight every MFE pair, splitting seed/MFE cell borders into orange and green (PR #259 review) @@ -68,6 +74,9 @@ ## Technical changes and Optimizations +- share Python companion-tool command execution, input validation and atomic + output handling after rebasing onto IntaRNA-fuse (PR #298 review) + - reuse MFE seed-extension scratch matrices and bound predecessor scans in exact and heuristic prediction; evaluate physical energies directly (PR #276 review) @@ -210,6 +219,17 @@ and accessibility models, and the output can include distinct structures, energies, restricted partition sums, and trackers. +# Auxiliary Python output tools + +- handle symlinked temporary directories in companion executable-lookup tests, + including macOS /var aliases (PR #298) + +- align companion tools on Python 3.10+, declare pip dependencies, install shared + modules, and exercise source/installed/native commands in CI (PR #298) + +- replace both R scripts with Python GEV/FDR estimation and region plotting; + document optional dependencies, migration and regression tests (issue #295) + # IntaRNA_plotRegions.R - fix deprecation warning @@ -226,6 +246,92 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261024 Alexander Mitrofanov + * python/tests/test_common_tools.py + * Compare executable lookup against canonical fixture paths so macOS /var + aliases do not cause false failures after the helper resolves symlinks. + + Exercise sibling, source-tree fallback, missing and frozen executable + lookup through an explicit directory symlink on every tested platform. + Reproduce and verify the macOS CI failure using an aliased TMPDIR (PR #298). + +261024 Alexander Mitrofanov + * python/setup.py, python/pyproject.toml, python/MANIFEST.in, python/*/Makefile.am, + python/Makefile.am, configure.ac, python/intarna_tools.py + * Declare Python 3.10+ and NumPy/SciPy dependencies, offer optional plotting, + and install shared helpers alongside standalone scripts and Python modules. + * Replace legacy setup.py installation with an opt-in pip hook honoring DESTDIR + without fetching dependencies or uninstalling packages in another prefix. + * python/bin/build.py, python/intarnapvalue/test.py, python/README.md, + python/bin/README.md, python/DEVELOPMENT.md, python/requirements-output-tools.txt + * Make executable builds independent of the working directory and fail visibly; + refresh documentation, source-test entry point and audit/performance notes. + * Preserve direct invocation of the shuffle utility with argparse validation, + optional sample counts and local random seeds. + * python/copomus/measures.py, python/CopomuS.py, python/tests/*.py, + .github/workflows/{build,python-output-tools}.yml, .gitignore + * Validate shared profiles and thresholds; test malformed/NA profile data, + installed lookup, atomic writes, native sampling and mutation workflows. + * Exercise all companion tools and installation in Python/C++ CI, including + the dependent PR base; exclude Python caches/build outputs (PR #298). + +261024 Alexander Mitrofanov + * python/intarnapvalue/intarna_pvalue.py, python/intarnapvalue/__main__.py + * Independently shuffle the correct target and query in pairwise batches; + apply the same parameter file to observed and shuffled predictions. + * Bound favorable-interaction resampling with --max-attempts, fail immediately + on command errors, and validate sample IDs/energies before publication. + * Normalize wrapped single-record FASTA and DNA, validate counts, preserve + zero energies, isolate RNG state and keep score ordering deterministic. + * Evaluate fitted CDFs instead of numerical integration; retain historical + raw-energy models, validate fits and expose empirical/none modes (PR #298). + * python/intarnapvalue/dinucl_shuffle.py, python/tests/test_intarnapvalue.py + * Handle singleton and invalid sequences; remove quadratic front-deletion + from Euler-path construction and support instance-local RNGs. + + Test composition, sampling, bounded errors, constraints and CDF equivalence. + * python/intarnapvalue/README.md, python/Makefile.am + * Document real options, conditional sampling, dependencies and API behavior. + +261024 Alexander Mitrofanov + * python/CopomuS.py, python/copomus/*.py, python/tests/test_copomus.py + * Validate sequences/FASTA, mutation positions and bases, delimiters and counts; + repair the explicit-selector call signature and CG/GC filter mismatch. + * Preserve zero energies and quoted CSV identifiers; fix signed-index conversion, + neighboring pairs across zero, and minimum-delta ordering at/below zero. + * Share profile contents across accessibility measures, clean private temporary + files, and reuse the wildtype prediction for all candidates (PR #298). + * Report errors on stderr and replace output only after successful computation. + * python/copomus/README.md, python/Makefile.am + * Correct option names, defaults, Python requirements and mutation/ranking docs; + distribute the focused regressions. + +261024 Alexander Mitrofanov + * python/intarna_tools.py, python/IntaRNA-fuse.py, python/IntaRNA_CSV_p-value.py, + python/IntaRNA_plotRegions.py, python/Makefile.am, tests/testIntaRNAFuse.py, + python/tests/test_output_tools.py + + Share dependency-free executable lookup, checked subprocess calls, RNA/FASTA + validation, CSV validation and atomic text publication between Python tools. + * Rebase onto feat/294-python-intarna-fuse; retain both migrations, their + installation rules and regression suites (PR #298). + +261023 Alexander Mitrofanov + * python/IntaRNA_CSV_p-value.py, python/IntaRNA_plotRegions.py + + Replace the R tools with standalone Python scripts preserving positional + arguments, negated-energy GEV tails, seven-digit rounding/BH correction, + ordered region rows, inclusive coverage and plotting parameter files. + + Validate input and fitting; publish CSV atomically, preserve original CSV + fields, and support all previous figure formats without a display. + * python/tests/test_output_tools.py, .github/workflows/python-output-tools.yml + + Test analytical tails, known GEV samples, BH, CSV errors/round trips, + plotting parameters, coverage and all seven output formats in Python CI. + * python/README.md, python/requirements-output-tools.txt, python/examples/*, + python/plotRegions.example.png, README.md, AGENTS.md + + Document Python dependencies, migration, statistical and rendering + differences, parameter defaults and examples; replace README R snippets. + * Makefile.am, configure.ac, python/Makefile.am, R/* + * Install/distribute the Python replacements, documentation and fixtures. + - Remove the R scripts and their build directory. + * Resolves https://github.com/BackofenLab/IntaRNA/issues/295 + 261023 Alexander Mitrofanov * python/IntaRNA-fuse.py, tests/testIntaRNAFuse.py, tests/runIntaRNAFuse.sh + Port the non-overlapping-site compatibility check and joint-accessibility diff --git a/Makefile.am b/Makefile.am index 1387c875..a44f9d79 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 R tests doc . +SUBDIRS = src python tests doc . # list of all personalities to be installed PERSONALITIES=@PERSONALITIES@ diff --git a/R/IntaRNA_CSV_p-value.R b/R/IntaRNA_CSV_p-value.R deleted file mode 100644 index 23d45326..00000000 --- a/R/IntaRNA_CSV_p-value.R +++ /dev/null @@ -1,150 +0,0 @@ -#!/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/IntaRNA_plotRegions.R b/R/IntaRNA_plotRegions.R deleted file mode 100644 index a805f5ea..00000000 --- a/R/IntaRNA_plotRegions.R +++ /dev/null @@ -1,194 +0,0 @@ -#!/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=expansion(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(linewidth=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 deleted file mode 100644 index f0562006..00000000 --- a/R/Makefile.am +++ /dev/null @@ -1,8 +0,0 @@ - -# 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 deleted file mode 100644 index 41162344..00000000 --- a/R/README.md +++ /dev/null @@ -1,65 +0,0 @@ - -# Auxiliary R scripts of the IntaRNA package - -- [IntaRNA_CSV_p-value.R](#intarna_csv_p-valuer) -- [IntaRNA_plotRegions.R](#intarna_plotregionsr) - - -# `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 `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` - with suffix `1` or `2` to plot target or query regions, respectively -2. `1` or `2` to select whether to plot target or query regions -3. output file name with a file-format-specific suffix from `.pdf`, `.png`, - `.svg`, `.eps`, `.ps`, `.jpeg`, `.tiff` - -An example is given below, when calling -```bash -Rscript --vanilla IntaRNA_plotRegions.R pred.csv 1 plotRegions.example.png -``` - -with `pred.csv` containing -``` -id1;start1;end1;id2;start2;end2 -b0001;266;273;query;116;123 -b0002;204;231;query;85;111 -b0003;229;262;query;96;125 -b0004;265;300;query;10;38 -b0005;281;295;query;5;22 -``` - -will produce the output -![IntaRNA_plotRegions.R example](plotRegions.example.png) - diff --git a/R/plotRegions.example.png b/R/plotRegions.example.png deleted file mode 100644 index 5f6f080b..00000000 Binary files a/R/plotRegions.example.png and /dev/null differ diff --git a/README.md b/README.md index 4af23ea2..e5aed90e 100644 --- a/README.md +++ b/README.md @@ -126,10 +126,9 @@ The following topics are covered by this documentation: - [Scaling factors for partition function computation for accessibility estimation](#accPfScale) - [Read/write accessibility from/to file or stream](#accFromFile) - [Library for integration in external tools](#lib) -- [Auxiliary R scripts for output visualization etc.](/R) -- [Auxiliary python scripts for IntaRNA-based pipelines](/python) +- [Auxiliary Python tools for output analysis, visualization and pipelines](python/README.md) - [IntaRNA-fuse.py](python/README.md#intarna-fusepy) combines compatible - interaction sites using joint accessibility costs (requires Python 3 and + interaction sites using joint accessibility costs (requires Python 3.10+ and RNAfold; installed by `make install`). @@ -1963,44 +1962,40 @@ pairs using `--out=pMinE:MYPAIRMINE.csv`. This generates a CSV file (`;`separate holding for each index pair the minimal energy of any interaction covering this index combination or `NA` if no covers it at all. -This information can be visualized with your preferred program. In the following, -the provided R call is used to generate a heatmap visualization of the -RNA-RNA interaction possibilities. - -```R -# read data, skip first column, and replace NA and E>0 values with 0 -d <- read.table("MYPAIRMINE.csv",header=T,sep=";"); -d <- d[,2:ncol(d)]; -d[is.na(d)] = 0; -d[d>0] = 0; -# plot -image( 1:nrow(d), 1:ncol(d), as.matrix(d), col = heat.colors(100), xlab="index in sequence 1", ylab="index in sequence 2"); -box(); +This information can be visualized with Python, NumPy and Matplotlib (see the +[optional Python dependencies](python/README.md#dependencies-and-installation)). +The following example uses the coordinate labels from the CSV and masks missing +or positive energies as zero: + +```python +import csv +import numpy as np +import matplotlib.pyplot as plt + +with open("MYPAIRMINE.csv", newline="") as handle: + rows = list(csv.reader(handle, delimiter=";")) +# Labels have the form nucleotide_position; split at the final underscore. +query = [int(label.rsplit("_", 1)[1]) for label in rows[0][1:]] +target = [int(row[0].rsplit("_", 1)[1]) for row in rows[1:]] +energy = np.array([[float(value) if value != "NA" else 0 for value in row[1:]] + for row in rows[1:]]) +energy = np.minimum(energy, 0) +figure, axis = plt.subplots() +heatmap = axis.pcolormesh(target, query, energy.T, shading="nearest", cmap="hot") +axis.set_xlabel("index in sequence 1") +axis.set_ylabel("index in sequence 2") +figure.colorbar(heatmap, ax=axis, label="energy [kcal/mol]") +figure.savefig("pair-minE.png", dpi=150, bbox_inches="tight") ``` The following plot (for the [minimal energy profile](#profileMinE) example from -above) reveals, that the alternative stable (*E*<0) interactions all involve the -mfe-site in the second sequence and are thus less likely to occure. +above) reveals that the alternative stable (*E*<0) interactions all involve the +MFE site in the second sequence and are thus less likely to occur. ![Minimal interaction energy index pair information](./doc/figures/pair-minE.png) -A more sophisticated plot can be done using `ggplot2` from the `tidyverse` package. - -```R -library(tidyverse) - -read_delim("pairMinE-t1q1.csv", delim=";") %>% - replace(is.na(.), 0) %>% - rename( target = minE ) %>% - pivot_longer( cols=-target, names_to = "query", values_to = "E" ) %>% - separate( query, into = c("qSeq","qIdx"), sep = "_") %>% - separate( target, into = c("tSeq","tIdx"), sep = "_") %>% - mutate(across(c(tIdx,qIdx,E),as.numeric)) %>% - mutate(E = ifelse( E>0, 0, E)) %>% - ggplot(aes(x=tIdx,y=qIdx, fill=E)) + - geom_tile() + - scale_fill_gradient(low="blue",high="red") -``` +For interval and coverage plots from standard interaction CSV output, use +[IntaRNA_plotRegions.py](python/README.md#intarna_plotregionspy). [![up](doc/figures/icon-up.28.png) back to overview](#overview) diff --git a/configure.ac b/configure.ac index fd470280..db1bf57b 100644 --- a/configure.ac +++ b/configure.ac @@ -8,7 +8,7 @@ BOOST_REQUIRED_VERSION=1.50.0 # minimal Vienna RNA version providing the logger locking API VRNA_REQUIRED_VERSION=2.7.0 # minimal required version of python interpreter (to install intarnapvalue) -PYTHON_REQUIRED_VERSION=3.6 +PYTHON_REQUIRED_VERSION=3.10 AC_CANONICAL_HOST @@ -570,7 +570,6 @@ AC_CONFIG_FILES([src/Makefile]) AC_CONFIG_FILES([src/IntaRNA/Makefile]) AC_CONFIG_FILES([src/IntaRNA/intarna_config.h]) AC_CONFIG_FILES([src/bin/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/python/CopomuS.py b/python/CopomuS.py index c5b9d7f3..d2460fc9 100644 --- a/python/CopomuS.py +++ b/python/CopomuS.py @@ -5,16 +5,15 @@ __intarna_min__ = '3.1.2' -__python_min__ = '3.7.0' -import os +import math import sys -from platform import python_version import csv -from tempfile import gettempdir -from argparse import ArgumentParser, FileType, SUPPRESS +from tempfile import TemporaryDirectory +from pathlib import Path +from argparse import ArgumentParser, SUPPRESS from typing import List from copomus.mutation import Mutation from copomus.IntaRNA import IntaRNA @@ -23,9 +22,7 @@ from copomus.mutation_generators import get_generator -if IntaRNA.version_tuple(python_version()) < IntaRNA.version_tuple(__python_min__): - sys.stderr.write(f'CopomuS requires python >= {__python_min__}\n') - sys.exit(1) +from intarna_tools import atomic_text_output, nonnegative_integer, read_sequence class CopomuS: @@ -45,9 +42,9 @@ def __init__(self): self.alpha, self.beta = 0.0, 0.0 self.threads = 0 - def main(self): + def main(self, argv=None): # PARSE COMMAND LINE ARGS - self.parse_cmd_args() + self.parse_cmd_args(argv) # CHECK INTARNA BIN AND VERSION IntaRNA.check_binary(__intarna_min__) @@ -64,8 +61,7 @@ def main(self): bp = func(bp, self.query, self.target, self.qidxpos0, self.tidxpos0) if not bp: - sys.stderr.write('There are no base pair candidates left after filtering!\n') - sys.exit(1) + raise ValueError('There are no base pair candidates left after filtering') # CREATE MUTATIONS if not self.mutation_encoding: @@ -89,7 +85,7 @@ def main(self): # OUTPUT self.output_csv(mutations) - def parse_cmd_args(self): + def parse_cmd_args(self, argv=None): defaults = { 'measures': ['mfeCover', 'E', 'minDeltaE'], 'selector': 'mfe', @@ -101,8 +97,8 @@ def parse_cmd_args(self): } parser = ArgumentParser(description='Checks different measures for rating mutations') - parser.add_argument('-q', '--query', dest='query', required=True, help='The query sequence.') - parser.add_argument('-t', '--target', dest='target', required=True, help='The target sequence.') + parser.add_argument('-q', '--query', dest='query', required=True, help='Query RNA/DNA sequence or single-record FASTA file.') + parser.add_argument('-t', '--target', dest='target', required=True, help='Target RNA/DNA sequence or single-record FASTA file.') parser.add_argument('--qIdxPos0', dest='qidxpos0', default=1, type=int, help='The starting index for the query. (Default: 1)') parser.add_argument('--tIdxPos0', dest='tidxpos0', default=1, type=int, @@ -120,14 +116,14 @@ def parse_cmd_args(self): default=[], help='Filters candidate base pairs, can be used multiple times. ' f"(Default: {defaults['filters']})", - choices=['GU', 'AU', 'CG', 'lp', 'lpMfe', 'he', 'heMfe']) + choices=['GU', 'AU', 'CG', 'GC', 'lp', 'lpMfe', 'he', 'heMfe']) parser.add_argument('-g', '--generator', dest='mutation_generator', default=defaults['generator'], choices=['flip', 'any'], help='Defines the method used for generating mutated sequences. ' f"(Default: {defaults['generator']})") parser.add_argument('--mutationEncoding', dest='mutation_encoding', help='Allows direct candidate selection by specifying a mutation encoding. ' 'Overwrites options -s, -f, and -g') - parser.add_argument('-o', '--output', dest='output', nargs='?', type=FileType('w'), default=sys.stdout, + parser.add_argument('-o', '--output', dest='output', default=None, help='Which file the output should be written to. (Default: STDOUT)') parser.add_argument('-d', '--delimiter', dest='delimiter', default='\t', help='Defines the delimiter used to separate columns in the output, default tab. ' @@ -135,60 +131,59 @@ def parse_cmd_args(self): parser.add_argument('-p', '--parameterFile', dest='param_file', default='', help='Optional parameter file for IntaRNA to provide further parameters and ' 'prediction constraints.') - parser.add_argument('--threads', dest='threads', default=defaults['threads'], type=int, + parser.add_argument('--threads', dest='threads', default=defaults['threads'], type=nonnegative_integer, help='Threads used for IntaRNA call') parser.add_argument('--alpha', dest='alpha', type=float, default=defaults['alpha'], help=SUPPRESS) parser.add_argument('--beta', dest='beta', type=float, default=defaults['beta'], help=SUPPRESS) - args = parser.parse_args() - - mm = [] # filter out any duplicate entries for measures - args.measures = defaults['measures'] if not args.measures else args.measures - for m in args.measures: - if m not in mm: - mm.append(m) - args.measures = mm - - ff = [] # filter out any duplicate entries for filters - args.candidate_filters = defaults['filters'] if not args.candidate_filters else args.candidate_filters - for f in args.candidate_filters: - if f not in ff: - ff.append(f) - args.candidate_filters = ff - - args.query = args.query.upper().replace('T', 'U') # capitalize sequences and replace T -> U - args.target = args.target.upper().replace('T', 'U') - - if args.param_file and not os.path.exists(args.param_file): - sys.stderr.write(f'Cannot find provided parameter file: {args.param_file}\n') - sys.exit(1) - - for key, value in args.__dict__.items(): - self.__setattr__(key, value) + args = parser.parse_args(argv) + args.measures = list(dict.fromkeys(args.measures or defaults['measures'])) + args.candidate_filters = list(dict.fromkeys(args.candidate_filters)) + try: + args.query = read_sequence(args.query) + args.target = read_sequence(args.target) + if not all(math.isfinite(value) for value in (args.alpha, args.beta)): + raise ValueError('alpha and beta must be finite') + if args.param_file and not Path(args.param_file).is_file(): + raise ValueError(f'Cannot find provided parameter file: {args.param_file}') + if len(args.delimiter) != 1 or args.delimiter in '\r\n\0': + raise ValueError('delimiter must be one non-newline character') + except (OSError, ValueError) as error: + parser.error(str(error)) + for key, value in vars(args).items(): + setattr(self, key, value) def run_tests(self, mutations: List[Mutation]): - for mut in mutations: - outcsvcols = list(set([j for i in [mm.outcsvcols for mm in mut.measures] for j in i])) # combine lists - outcsvcols = ','.join(outcsvcols) # join columns - outcsvcols = 'id1' if not outcsvcols else outcsvcols # if no col -> select id1 - - out_pipes = [j for i in [m.out_pipes for m in mut.measures] for j in i] # combine lists - - # create args for IntaRNA call - out_param = [f'--out={x}:{os.path.join(gettempdir(), f"CopomuS_{x}")}_{os.getpid()}.temp' for x in out_pipes] - - for mode in ['ww', 'wm', 'mw', 'mm']: # go through each mutation mode - q = mut.qw if mode[0] == 'w' else mut.qm # get actual sequences - t = mut.tw if mode[1] == 'w' else mut.tm - - i = IntaRNA() # call IntaRNA, we have 4 calls per base pair, one for each mutation mode - data = i.run(q, t, mut.qidxpos0, mut.tidxpos0, outcsvcols, self.threads, 1, self.param_file, out_param) - - for m in mut.measures: # add results to each measure - m.calculate_score(data, mode, mut, self.param_file) - - for m in mut.measures: - m.post_compute() + # Evaluate wildtype once: its sequences and requested outputs are identical + # for every candidate. Profile contents are immutable and shared by measures. + wildtype = None + with TemporaryDirectory(prefix='CopomuS-') as directory: + for mut in mutations: + columns = list(dict.fromkeys(c for m in mut.measures for c in m.outcsvcols)) + pipes = list(dict.fromkeys(p for m in mut.measures for p in m.out_pipes)) + paths = {pipe: Path(directory) / pipe for pipe in pipes} + out_params = [f'--out={pipe}:{path}' for pipe, path in paths.items()] + for mode in ('ww', 'wm', 'mw', 'mm'): + if mode == 'ww' and wildtype is not None: + data, profiles = wildtype + else: + query = mut.qw if mode[0] == 'w' else mut.qm + target = mut.tw if mode[1] == 'w' else mut.tm + # Remove prior outputs so a missing profile cannot reuse stale data. + for path in paths.values(): + path.unlink(missing_ok=True) + data = IntaRNA().run(query, target, mut.qidxpos0, mut.tidxpos0, + ','.join(columns) or 'id1', self.threads, + 1, self.param_file, out_params) + profiles = {pipe: path.read_text(encoding='utf-8') + for pipe, path in paths.items()} + if mode == 'ww': + wildtype = data, profiles + for measure in mut.measures: + measure.profiles = profiles + measure.calculate_score(data, mode, mut, self.param_file) + for measure in mut.measures: + measure.post_compute() def calculate_ranks(self, mutations: List[Mutation]): # Calculate total ranking order @@ -215,8 +210,16 @@ def calculate_ranks(self, mutations: List[Mutation]): mutations[mut_id].ranks[m_id] = rank def output_csv(self, mutations: List[Mutation]): + # Open a destination only after all predictions and ranking succeeded. + if self.output: + with atomic_text_output(self.output) as handle: + self.write_csv(mutations, handle) + else: + self.write_csv(mutations, sys.stdout) + + def write_csv(self, mutations: List[Mutation], handle): # CREATE CSV WRITER - writer = csv.writer(self.output, delimiter=self.delimiter) + writer = csv.writer(handle, delimiter=self.delimiter, lineterminator="\n") # WRITE CSV headers = [j for i in @@ -237,6 +240,14 @@ def output_csv(self, mutations: List[Mutation]): writer.writerow([str(m), m.ranks['total'], m.query_n, m.target_n, m.bp_ww, m.bp_mm] + values) # write csv +def main(argv=None): + try: + CopomuS().main(argv) + except (OSError, ValueError, RuntimeError, csv.Error) as error: + print(f'Error: {error}', file=sys.stderr) + return 1 + return 0 + + if __name__ == '__main__': - c = CopomuS() - c.main() + raise SystemExit(main()) diff --git a/python/DEVELOPMENT.md b/python/DEVELOPMENT.md new file mode 100644 index 00000000..9e079875 --- /dev/null +++ b/python/DEVELOPMENT.md @@ -0,0 +1,91 @@ +# Companion-tool audit (PR #298) + +The branch builds on `feat/294-python-intarna-fuse` (PR #296). The R-to-Python +migration, shared infrastructure, CopomuS fixes, IntaRNApvalue fixes, and +packaging/validation are separate commits. + +## Survey and common conventions + +| Area | Review outcome | +| --- | --- | +| IntaRNA-fuse | Retain its fusion model, decimal energies and literal-sequence interface; share checked subprocess calls and executable lookup. | +| CSV p-values / region plots | Retain positional interfaces and the R-compatible formulas; share strict CSV header/row validation and atomic CSV publication. | +| CopomuS | Share RNA/FASTA validation and subprocess handling; fix explicit mutation dispatch, CG/GC filtering, signed indexing, zero/NA handling, ranking, and profile ownership. | +| IntaRNApvalue | Correct both-strand sampling, consistent constraints and failure handling; share FASTA and executable handling; expose empirical estimates and preserve fitted models. | +| Dinucleotide shuffling | Normalize at entry, handle length one, use an optional RNG, and consume adjacency lists with deques. | +| Packaging / executable builder | Declare Python 3.10 and real scientific dependencies; install shared modules with both pip and Autotools; replace shell-based PyInstaller calls and cleanup tied to the caller's directory. | + +`intarna_tools.py` contains only standard-library I/O and validation. Scientific +calculations stay in their tools: the genome-wide CSV fit uses negated energies, +whereas IntaRNApvalue retains its historical raw-energy fitted distributions. +Combining these fits as one helper would change one tool's statistical model. +Likewise, fusion requires literal sequences for constraints, while CopomuS and +IntaRNApvalue can read a single FASTA record. + +All entry points have side-effect-free imports and `main(argv)` interfaces. +Existing options remain; `GC` is an alias of CopomuS's advertised `CG` filter. +IntaRNApvalue adds `--max-attempts` and `-d empirical` (`none` alias). Sampled +p-values remain conditional on a reported interaction. Errors never enter +numeric/CSV stdout, and CopomuS/output-analysis CSV files are replaced only after +successful computation. See each tool's README for interpretation and examples. + +## Correctness and performance checks + +The Python unittest suite covers malformed/quoted CSV, zero/NA values, signed +coordinate round trips, explicit mutation bounds, all plot formats, independent +shuffle pairing, composition preservation, retry limits, subprocess errors, +parameter forwarding, and old-integration/new-CDF equivalence. Native checks +exercise both/single-strand sampling, deterministic output across thread counts, +profile combinations, explicit mutations and failure limits. The fusion suite +retains its independent synthetic energy oracle and optional native check. + +CopomuS reuses the one wildtype prediction and immutable profile contents. +For N mutation candidates, measurement calls decrease from 4N to 3N+1, with +candidate-selection calls unchanged. A two-candidate regression verifies seven +calls and complete temporary-file cleanup. Filter membership uses sets without +changing candidate order. No native predictor code or defaults are changed. + +Local comparison on Linux, Python 3.12.7, the same GCC 14 release IntaRNA binary, +ViennaRNA 2.7.2, and one IntaRNA thread (three runs each): + +```bash +python3 CopomuS.py \ + -q GCUGAAAAACAUAACCCAUAAAAUGCUAGCUGUACCAGGAACCA \ + -t GGUUUCUUCGCCUCUGCGUUCACCAAAGUGUUCACCC --threads 1 -m E +``` + +Comparing pre-review commit `b4a9c99` with the revised sources, all columns, +mutation order, ranks and energy values match after normalizing line endings. +Wall times measured by `/usr/bin/time -f '%e %M'` were 2.90/2.94/2.97 seconds +before and 2.22/2.29/2.23 seconds after. Peak RSS ranges were 17,704–17,712 KiB +and 17,584–17,712 KiB. These are workload-specific measurements made during a +concurrent build, not a general throughput guarantee. + +For the shuffle loop, use `s = 'AGGAUGGGGGACCUUCUCGUUCCU' * 4000`, reset +`random.seed(7)` before each call, and time `dinucl_shuffle(s)` with +`time.perf_counter()`. Three runs took 0.239–0.243 seconds before and 0.093–0.095 +seconds after, with byte-identical shuffles. Deque consumption removes repeated +front deletion from Python lists; mononucleotide/dinucleotide counts and endpoint +preservation are tested separately. Both versions retain adjacency storage +linear in sequence length. + +## Reproducing validation + +```bash +python3 -m pip install -r python/requirements-output-tools.txt +python3 -m unittest discover -s python/tests -v +python3 tests/testIntaRNAFuse.py +bash autotools-init.sh +./configure # supply your compiler and dependency prefixes as needed +make tests -j2 +make -C python check-python-tools +make install +make dist +``` + +Place a native IntaRNA on PATH to enable the optional native Python tests and +RNAfold on PATH for the fusion integration test. CI covers Python 3.10 with +minimum dependency versions, Python 3.13 with current dependencies, installed +entry points, and native execution in the C++ build matrix. PyInstaller itself +is optional; tests check its invocation and failure propagation without building +a platform-specific frozen executable. diff --git a/python/IntaRNA-fuse.py b/python/IntaRNA-fuse.py index abfc24d8..25e0844c 100755 --- a/python/IntaRNA-fuse.py +++ b/python/IntaRNA-fuse.py @@ -11,11 +11,12 @@ import io from pathlib import Path import re -import shutil -import subprocess import sys +from intarna_tools import find_program, positive_integer, run_command + + ESSENTIAL_COLUMNS = ('E', 'start1', 'end1', 'start2', 'end2', 'ED1', 'ED2', 'E_init') DEFAULT_COLUMNS = ESSENTIAL_COLUMNS + ('hybridDB',) @@ -27,13 +28,6 @@ def sequence(value): return value -def positive_integer(value): - number = int(value) - if number < 1: - raise argparse.ArgumentTypeError('expected a positive integer') - return number - - def parse_args(argv=None): parser = argparse.ArgumentParser( description='Report non-overlapping IntaRNA sites and compatible fused interactions.', @@ -55,35 +49,6 @@ def parse_args(argv=None): return args -def find_program(name): - """Prefer installed siblings, then the source-tree binary, then PATH.""" - script_dir = Path(__file__).resolve().parent - candidates = [script_dir / name] - if name == 'IntaRNA': - candidates.append(script_dir.parent / 'src' / 'bin' / name) - for candidate in candidates: - found = shutil.which(str(candidate)) - if found: - return found - found = shutil.which(name) - if not found: - raise RuntimeError('{} executable not found beside the script or on PATH'.format(name)) - return found - - -def run_command(command, input_text=None): - """Pass sequences and constraints without shell interpretation; check failures.""" - result = subprocess.run(command, input=input_text, universal_newlines=True, - stdout=subprocess.PIPE, stderr=subprocess.PIPE) - if result.returncode: - detail = (result.stderr or result.stdout).strip() - raise RuntimeError('{} exited with status {}: {}'.format( - Path(command[0]).name, result.returncode, detail)) - if result.stderr: - sys.stderr.write(result.stderr) - return result.stdout - - def energy(value): """Read a finite energy without introducing binary rounding in the sum.""" try: diff --git a/python/IntaRNA_CSV_p-value.py b/python/IntaRNA_CSV_p-value.py new file mode 100755 index 00000000..c1957c09 --- /dev/null +++ b/python/IntaRNA_CSV_p-value.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +"""Append genome-wide GEV p-values and Benjamini-Hochberg FDRs to IntaRNA CSV.""" + +import argparse +import csv +import math + + +from intarna_tools import atomic_text_output, csv_header, csv_rows + + +def fit_gev(energies): + """Fit negated energies; return SciPy's (shape=-xi, location, scale).""" + import numpy as np + from scipy.optimize import minimize + from scipy.stats import genextreme + + # The background must contain enough distinct, finite observations to fit. + values = np.asarray(energies, dtype=float) + if not np.all(np.isfinite(values)) or np.unique(values).size < 3: + raise ValueError("GEV fitting requires at least three distinct finite energies") + data = -values + scale = float(np.std(data, ddof=1) * math.sqrt(6) / math.pi) + location = float(np.mean(data) - 0.57722 * scale) + if not math.isfinite(scale) or scale <= 0 or not math.isfinite(location): + raise ValueError("energy range is too large for GEV fitting") + + # Retain the R script's initial xi=0.1 and moment-based scale/location. + # SciPy uses the opposite shape sign and handles the Gumbel limit (xi=0). + def optimizer(function, initial, args=(), disp=0): + result = minimize(function, initial, args=args, method="Nelder-Mead", + options={"maxiter": 10000, "maxfev": 20000, + "xatol": 1e-8, "fatol": 1e-8}) + if not result.success: + raise ValueError("GEV fit did not converge: " + result.message) + return result.x + + params = genextreme.fit(data, -0.1, loc=location, scale=scale, optimizer=optimizer) + if (not np.all(np.isfinite(params)) or params[2] <= 0 + or not np.all(np.isfinite(genextreme.logpdf(data, *params)))): + raise ValueError("GEV fit is not finite or does not cover all energies") + return params + + +def gev_pvalues(energies, params): + """Return P(E <= observed E), using the upper tail of the negated energies.""" + import numpy as np + from scipy.stats import genextreme + + # sf avoids cancellation for significant interactions and handles support ends. + values = genextreme.sf(-np.asarray(energies, dtype=float), *params) + if not np.all(np.isfinite(values)): + raise ValueError("GEV fit produced nonfinite p-values") + return values + + +def adjust_bh(pvalues): + """Adjust p-values by Benjamini-Hochberg, retaining the original row order.""" + result = [0.0] * len(pvalues) + previous = 1.0 + order = sorted(range(len(pvalues)), key=lambda i: pvalues[i]) + for rank in range(len(order), 0, -1): + index = order[rank - 1] + previous = min(previous, pvalues[index] * len(order) / rank) + result[index] = previous + return result + + +def extend_csv(source, destination, energy_column="E"): + """Validate and fit the entire table before atomically publishing the output.""" + with open(source, newline="", encoding="utf-8-sig") as handle: + reader = csv.reader(handle, delimiter=";", strict=True) + header = csv_header(reader) + if energy_column not in header: + raise ValueError("energy column {!r} is missing".format(energy_column)) + if "p-value" in header or "fdr" in header: + raise ValueError("CSV already contains a 'p-value' or 'fdr' column") + column = header.index(energy_column) + rows, energies = [], [] + for row in csv_rows(reader, header): + try: + energy = float(row[column]) + except ValueError: + raise ValueError("invalid energy on CSV line {}".format(reader.line_num)) from None + if not math.isfinite(energy): + raise ValueError("nonfinite energy on CSV line {}".format(reader.line_num)) + rows.append(row) + energies.append(energy) + + # Preserve the historical seven-decimal rounding *before* BH adjustment. + pvalues = [round(float(p), 7) for p in gev_pvalues(energies, fit_gev(energies))] + fdrs = adjust_bh(pvalues) + with atomic_text_output(destination) as handle: + writer = csv.writer(handle, delimiter=";", lineterminator="\n") + writer.writerow(header + ["p-value", "fdr"]) + for row, pvalue, fdr in zip(rows, pvalues, fdrs): + writer.writerow(row + [format(pvalue, ".7f"), format(fdr, ".7f")]) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__, epilog=( + "Use a sufficiently large, diverse genome-wide prediction set as the background. " + "If output is omitted, input is overwritten after a successful fit. " + "Requires Python 3.10+, NumPy and SciPy (see python/requirements-output-tools.txt).")) + parser.add_argument("input", help="semicolon-separated IntaRNA CSV") + parser.add_argument("output", nargs="?", help="output CSV (default: overwrite input)") + parser.add_argument("energy_column", nargs="?", default="E", help="energy column (default: E)") + args = parser.parse_args(argv) + try: + extend_csv(args.input, args.output or args.input, args.energy_column) + except ImportError as error: + parser.exit(1, "Missing Python dependency: {}. Install numpy and scipy.\n".format(error)) + except (OSError, ValueError, RuntimeError, csv.Error) as error: + parser.exit(1, "Error: {}\n".format(error)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/python/IntaRNA_plotRegions.py b/python/IntaRNA_plotRegions.py new file mode 100755 index 00000000..2cb59bdb --- /dev/null +++ b/python/IntaRNA_plotRegions.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python3 +"""Plot IntaRNA interaction regions and their smoothed nucleotide coverage.""" + +import argparse +import csv +import math +from pathlib import Path + + +from intarna_tools import csv_header, csv_rows + + +FORMATS = {".pdf", ".png", ".svg", ".eps", ".ps", ".jpeg", ".tiff"} + + +def read_parameters(selector): + """Read a strand selector or the original key=value plotting parameters.""" + parameters = {"title": None, "xvline": None, "xmin": None, "xmax": None, "rows": 200} + if selector in ("1", "2"): + parameters.update({key: key + selector for key in ("id", "start", "end")}) + return parameters + + # Comments and optional values retain the original parameter-file syntax. + seen = set() + with open(selector, encoding="utf-8-sig") as handle: + for line_number, line in enumerate(handle, 1): + line = line.split("#", 1)[0].strip() + if not line: + continue + key, separator, value = line.partition("=") + key, value = key.strip(), value.strip() + if not separator or key not in {"id", "start", "end", *parameters}: + raise ValueError("unknown or malformed parameter on line {}".format(line_number)) + if key in seen: + raise ValueError("duplicate parameter: " + key) + seen.add(key) + parameters[key] = value + for key in ("id", "start", "end"): + if not parameters.get(key): + raise ValueError("parameter file must define " + key) + for key in ("xmin", "xmax", "xvline"): + value = parameters[key] + if value is None or value in ("", "NA"): + parameters[key] = None + else: + parameters[key] = float(value) + if not math.isfinite(parameters[key]): + raise ValueError(key + " must be finite") + try: + parameters["rows"] = int(parameters["rows"]) + except ValueError: + raise ValueError("rows must be a positive integer") from None + if parameters["rows"] < 1: + raise ValueError("rows must be a positive integer") + return parameters + + +def read_regions(source, parameters): + """Read the selected rows in input order, keeping inclusive signed positions.""" + regions = [] + with open(source, newline="", encoding="utf-8-sig") as handle: + reader = csv.reader(handle, delimiter=";", strict=True) + header = csv_header(reader) + for key in ("id", "start", "end"): + if parameters[key] not in header: + raise ValueError("missing CSV column: " + parameters[key]) + columns = [header.index(parameters[key]) for key in ("id", "start", "end")] + for row in csv_rows(reader, header): + try: + start, end = (int(row[index]) for index in columns[1:]) + except ValueError: + raise ValueError("positions must be integers on CSV line {}".format(reader.line_num)) from None + regions.append((row[columns[0]], start, end)) + if len(regions) == parameters["rows"]: + break + if not regions: + raise ValueError("CSV contains no interaction regions") + return regions + + +def coverage_density(regions, xmin, xmax): + """Gaussian count density of inclusive covered positions inside the x limits.""" + import numpy as np + from scipy.stats import norm + + # ggplot's original scale limits filter positions before estimating density. + positions = [np.arange(max(min(start, end), math.ceil(xmin)), + min(max(start, end), math.floor(xmax)) + 1, dtype=float) + for _, start, end in regions] + positions = np.concatenate(positions) + x = np.linspace(xmin, xmax, 512) + if positions.size == 0: + return x, np.zeros_like(x) + + # R's default bw.nrd0 bandwidth, including its zero-spread fallback. + spread = float(np.std(positions, ddof=1)) if positions.size > 1 else 0.0 + iqr = float(np.subtract(*np.percentile(positions, [75, 25]))) + bandwidth = 0.9 * (min(spread, iqr / 1.34) or spread or abs(positions[0]) or 1.0) + bandwidth *= positions.size ** -0.2 + # Evaluate in chunks so the number of plotted bases does not multiply memory by 512. + density = np.zeros_like(x) + for offset in range(0, len(positions), 4096): + distances = (x[:, None] - positions[None, offset:offset + 4096]) / bandwidth + density += norm.pdf(distances).sum(axis=1) / bandwidth + return x, density + + +def create_figure(regions, parameters): + """Create aligned coverage and region panels without requiring a display.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + # Preserve the input ordering and independent rows even for duplicate IDs. + xmin = parameters["xmin"] + xmax = parameters["xmax"] + if xmin is None: + xmin = min(min(start, end) for _, start, end in regions) + if xmax is None: + xmax = max(max(start, end) for _, start, end in regions) + if xmin > xmax: + raise ValueError("xmin must not exceed xmax") + if xmin == xmax: + xmin, xmax = xmin - 0.5, xmax + 0.5 + x, density = coverage_density(regions, xmin, xmax) + region_height = max(2, len(regions) / 9) + figure, (coverage, intervals) = plt.subplots( + 2, 1, figsize=(6, 2 + region_height), sharex=True, + gridspec_kw={"height_ratios": [2, region_height]}, layout="constrained") + coverage.plot(x, density, color="black", linewidth=1) + coverage.set_ylabel("coverage") + coverage.set_ylim(bottom=0) + if parameters["title"]: + coverage.set_title(parameters["title"], fontweight="bold") + + for index, (_, start, end) in enumerate(regions): + intervals.plot([start, end], [index, index], color="dodgerblue", + linewidth=2, marker="o", markersize=3) + intervals.set_yticks(range(len(regions)), labels=[region[0] for region in regions], fontsize=7) + intervals.set_ylim(len(regions) - 0.5, -0.5) + intervals.grid(axis="y", color="lightgray", linewidth=0.7) + intervals.set_axisbelow(True) + intervals.set_xlabel("position") + intervals.axvline(xmin, color="black", linewidth=0.7) + for axis in (coverage, intervals): + axis.set_xlim(xmin, xmax) + axis.yaxis.tick_right() + axis.yaxis.set_label_position("right") + axis.spines["top"].set_visible(False) + axis.spines["left"].set_visible(False) + line = parameters["xvline"] + if line is not None and xmin <= line <= xmax: + axis.axvline(line, color="black", linewidth=0.7) + return figure + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__, epilog=( + "Parameter files accept id, start, end, title, xvline, xmin, xmax and rows (default: 200). " + "Requires Python 3.10+, NumPy, SciPy and Matplotlib; see python/requirements-output-tools.txt.")) + parser.add_argument("input", help="semicolon-separated IntaRNA CSV") + parser.add_argument("selector", help="1 (target), 2 (query), or a key=value parameter file") + parser.add_argument("output", help="plot filename: " + ", ".join(sorted(FORMATS))) + args = parser.parse_args(argv) + try: + if Path(args.output).suffix.lower() not in FORMATS: + raise ValueError("unsupported output extension; choose " + ", ".join(sorted(FORMATS))) + parameters = read_parameters(args.selector) + regions = read_regions(args.input, parameters) + figure = create_figure(regions, parameters) + import matplotlib.pyplot as plt + try: + figure.savefig(args.output, dpi=150) + finally: + plt.close(figure) + except ImportError as error: + parser.exit(1, "Missing Python dependency: {}. Install numpy, scipy and matplotlib.\n".format(error)) + except (OSError, ValueError, RuntimeError, csv.Error) as error: + parser.exit(1, "Error: {}\n".format(error)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/python/MANIFEST.in b/python/MANIFEST.in new file mode 100644 index 00000000..5d01a021 --- /dev/null +++ b/python/MANIFEST.in @@ -0,0 +1,6 @@ +include README.md DEVELOPMENT.md requirements-output-tools.txt +recursive-include copomus *.md *.jpg +recursive-include intarnapvalue *.md +recursive-include examples *.csv *.params +recursive-include tests *.py +include plotRegions.example.png diff --git a/python/Makefile.am b/python/Makefile.am index 681ba63e..722eab9a 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -1,32 +1,33 @@ - -# sub directories to check for Makefiles +# Companion modules and optional executable builder. SUBDIRS = copomus intarnapvalue bin - -# scripts to be installed in $(bindir) directly - dist_bin_SCRIPTS = \ CopomuS.py \ - IntaRNA-fuse.py + IntaRNA-fuse.py \ + IntaRNA_CSV_p-value.py \ + IntaRNA_plotRegions.py -EXTRA_DIST = README.md - -# auxiliary scripts to be part of distribution +# Standalone scripts import this module from their installed directory. +dist_pythonhelpers_DATA = intarna_tools.py +pythonhelpersdir = $(bindir) noinst_PYTHON = setup.py -# install python modules if HAVE_PYTHON - if INSTALL_INTARNAPVALUE - -# ensure pip is installed -# create python directory if non-existing with current prefix +# Install modules using the selected interpreter without fetching dependencies. +# Users install NumPy/SciPy in that interpreter's environment as documented. install-exec-hook: - @ $(PYTHON) -m ensurepip --default-pip - @ $(MKDIR_P) -p ${pythondir} - @ export PYTHONPATH=$PYTHONPATH:${pythondir}; \ - $(PYTHON) setup.py install --prefix=${prefix} - + $(PYTHON) -m pip install --ignore-installed --no-deps --no-build-isolation \ + --prefix="$(DESTDIR)$(prefix)" "$(abs_srcdir)" endif - endif + +EXTRA_DIST = README.md DEVELOPMENT.md requirements-output-tools.txt pyproject.toml MANIFEST.in \ + plotRegions.example.png examples/regions.csv examples/regions.params \ + tests/test_output_tools.py tests/test_copomus.py tests/test_intarnapvalue.py \ + tests/test_common_tools.py tests/test_native_tools.py + +.PHONY: check-output-tools check-python-tools +check-output-tools: check-python-tools +check-python-tools: + $(PYTHON) -m unittest discover -s $(srcdir)/tests -v diff --git a/python/README.md b/python/README.md index bf69709e..f6d0c219 100644 --- a/python/README.md +++ b/python/README.md @@ -1,5 +1,165 @@ +# Auxiliary Python tools -# auxiliary tools +All companion tools require **Python 3.10+**. They use `argparse` for `--help`, +write results to stdout or the requested output file, and report failures on +stderr with a nonzero exit status. Existing positional arguments and option +names remain supported. CopomuS and IntaRNApvalue accept RNA/DNA literals or +one-record FASTA files; IntaRNA-fuse retains literal sequences (including N). + +| Tool | Purpose | Additional dependencies | +| --- | --- | --- | +| [IntaRNA-fuse.py](#intarna-fusepy) | Combine compatible interaction sites | IntaRNA, RNAfold | +| [CopomuS.py](copomus/README.md) | Rank compensatory point mutations | IntaRNA >= 3.1.2 | +| [IntaRNApvalue](intarnapvalue/README.md) | Significance from shuffled sequences | IntaRNA; NumPy/SciPy for fits | +| [IntaRNA_CSV_p-value.py](#intarna_csv_p-valuepy) | Genome-wide energy p-values/FDRs | NumPy, SciPy | +| [IntaRNA_plotRegions.py](#intarna_plotregionspy) | Plot interaction regions | NumPy, SciPy, Matplotlib | + +Run from the checkout or install all commands/modules with +`python3 -m pip install './python[plot]'` from the repository root (omit `[plot]` +if plotting is not needed). Native IntaRNA and RNAfold are installed separately. +`make install` installs the four `.py` scripts with their supporting modules; +`--enable-intarnapvalue` additionally installs the Python package using pip with +no dependency downloads. That optional hook requires pip/setuptools in the +selected Python environment. Install scientific dependencies separately for it. +For a custom prefix outside that interpreter's normal search path, add +`PREFIX/lib/pythonX.Y/site-packages` to `PYTHONPATH` when invoking IntaRNApvalue; +using `pip install` in a virtual environment avoids that manual path setup. +The optional hook does not uninstall an existing package in another prefix. + +`intarna_tools.py` shares checked command execution, executable lookup, input +validation and atomic output between the scripts. Keep it beside scripts when +copying them manually; keep `copomus/` beside `CopomuS.py` as well. Installed +Python packages include the helper automatically. Imports perform no predictions. + +## Output analysis and visualization + +The former R tools are now standalone Python scripts in this folder. Replace +`Rscript --vanilla IntaRNA_CSV_p-value.R` with +`python3 IntaRNA_CSV_p-value.py`, and `Rscript --vanilla IntaRNA_plotRegions.R` +with `python3 IntaRNA_plotRegions.py`. Positional arguments are unchanged. +R and the R packages ggplot2, ggalt and cowplot are no longer needed. + +### Dependencies and installation + +Both scripts require **Python 3.10 or newer** and **NumPy >= 1.23**, +**SciPy >= 1.9**; plotting additionally requires **Matplotlib >= 3.6**. +From the repository root, install these optional dependencies in your Python +environment: + +```bash +python3 -m pip install -r python/requirements-output-tools.txt +``` + +Run the scripts directly from the checkout as shown below. `make install` also +installs both executable `.py` scripts in the configured `bin` directory; their +`#!/usr/bin/env python3` interpreter must have the dependencies installed. +Autotools does not install these Python dependencies. They are optional for the +C++ IntaRNA predictor. The fitted CSV and shuffled-sequence workflows have separate statistical models, +although they are distributed together by the Python package. + +### `IntaRNA_CSV_p-value.py` + +Estimate p-values for a **genome-wide** target screen from its observed energy +distribution. A sufficiently large and diverse collection of target energies is +required to serve as a background. This is not a significance test calibrated +for a small selection of interactions or a single target/query pair. + +The script reads semicolon-separated IntaRNA CSV, fits a generalized extreme +value (GEV) distribution to the **negated** energies by maximum likelihood, and +appends `p-value` and `fdr` columns. More negative energies receive smaller +p-values. FDRs use the Benjamini-Hochberg correction over all input rows. As in +the former script, p-values are rounded to seven decimal places **before** the +correction, and FDRs are also written with seven decimal places. + +```bash +# Overwrite the input only after successful validation and fitting. +python3 python/IntaRNA_CSV_p-value.py IntaRNA-output.csv +# Write a separate output file. +python3 python/IntaRNA_CSV_p-value.py IntaRNA-output.csv IntaRNA-output-with-pValue.csv +# Fit an existing custom energy column. +python3 python/IntaRNA_CSV_p-value.py IN.csv IN-pValue.csv E_norm +``` + +The default energy column is `E`. All original field values, column names and +row order are retained; fields containing semicolons or quotes are CSV-quoted. +Empty tables, fewer than three distinct energies, missing/duplicate columns, +nonfinite or nonnumeric energies, malformed rows, existing `p-value`/`fdr` +columns and unsuccessful fits are errors. The distinct-value check is a numerical +minimum, not a guarantee of a statistically adequate background. Errors leave +existing output files untouched. + +SciPy uses the opposite GEV shape sign to the old R/evir convention +(`c = -xi`); the script accounts for this and evaluates the survival function +of the negated energies. This also handles the zero-shape Gumbel limit and +extreme tails. The statistical model and initial estimates are retained, but +SciPy's optimizer and convergence tolerances can produce slightly different +fitted values from R; outputs are not guaranteed to match digit for digit. +See [SciPy's GEV documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.genextreme.html). + +### `IntaRNA_plotRegions.py` + +Visualize the regions covered by IntaRNA interactions in two aligned panels: +smoothed nucleotide coverage above and one blue interval per interaction below. +The first CSV row appears at the top. Repeated IDs remain separate rows. Both +interval endpoints are included in the coverage calculation; signed coordinates +and intervals specified in reverse order are supported. + +Arguments, in order: + +1. Semicolon-separated IntaRNA CSV containing `id1,start1,end1` or + `id2,start2,end2`. +2. `1` for target regions, `2` for query regions, or a parameter-file path. +3. Output filename ending in `.pdf`, `.png`, `.svg`, `.eps`, `.ps`, `.jpeg` + or `.tiff`. + +```bash +python3 python/IntaRNA_plotRegions.py python/examples/regions.csv 1 regions.png +python3 python/IntaRNA_plotRegions.py python/examples/regions.csv 2 query-regions.svg +python3 python/IntaRNA_plotRegions.py python/examples/regions.csv python/examples/regions.params query-regions.pdf +``` + +The [example CSV](examples/regions.csv) produces: + +![Interaction regions plotted with Python](plotRegions.example.png) + +By default, the first **200** interactions are plotted. A parameter file accepts +one `key = value` assignment per line, whitespace, blank lines and `#` comments. +`id`, `start` and `end` are required CSV column names; all other keys are optional: + +| Key | Meaning | Default | +| --- | --- | --- | +| `id` | CSV column used for row labels | Required | +| `start` | CSV column for integer interval starts | Required | +| `end` | CSV column for integer interval ends | Required | +| `title` | Title above the coverage panel | No title | +| `xvline` | Vertical reference position in both panels | None | +| `xmin` | Left x-axis limit | Smallest covered position | +| `xmax` | Right x-axis limit | Largest covered position | +| `rows` | Positive maximum number of interactions | 200 | + +An example is provided in [examples/regions.params](examples/regions.params). +Blank or `NA` values for optional x positions use the defaults. The coverage +curve uses a Gaussian kernel with R's original `bw.nrd0` bandwidth rule, scaled +by the number of covered bases (a smoothed count, not a fraction). As in the R +plot, x limits filter coverage samples before estimating the curve. Matplotlib +replaces the R graphics backend, so fonts, margins and density interpolation +can differ. Plots are generated without a graphical display. + +### Tests + +After installing the optional dependencies, run from the repository root: + +```bash +python3 -m unittest discover -s python/tests -v +``` + +In a configured build, `make -C python check-python-tools` runs the same tests +(`check-output-tools` remains an alias). +They cover GEV formulas and known-distribution fitting, BH correction, CSV +preservation and errors, every plot format, mutation coordinates/ranking/profiles, +shuffled sampling and failure handling. Native smoke tests run when IntaRNA is +available. CI checks minimum/current scientific dependencies on Python 3.10/3.13, +installed command entry points, and native behavior in the C++ build matrix. ## [IntaRNA-fuse.py](IntaRNA-fuse.py) @@ -10,7 +170,7 @@ overlaps the original MFE by at least two nucleotides on **each** strand, it reports a fused interaction using the suboptimal site and that constrained MFE. This preserves the compatibility criterion of the former Perl script. -Requires Python 3.6 or newer (standard library only), IntaRNA, and ViennaRNA's +Requires Python 3.10 or newer (standard library only), IntaRNA, and ViennaRNA's `RNAfold`. `make install` installs `IntaRNA-fuse.py` alongside `IntaRNA`; no pip installation or `--enable-intarnapvalue` is needed. The script looks for executables beside itself and on `PATH`, and also checks `../src/bin/IntaRNA` @@ -67,12 +227,12 @@ produce a nonzero exit status and a diagnostic on stderr, with no result CSV. ## [IntaRNApvalue](intarnapvalue) -A tool for calculating p-values of IntaRNA energy scores. +Calculate p-values using energies from **shuffled sequences**. This is a +separate workflow from fitting an existing genome-wide CSV with +`IntaRNA_CSV_p-value.py`; see its own dependency and usage documentation. ## [CopomuS](copomus) -The Compensatory Point Mutation Selector helps to identify the most promising -base pairs from an IntaRNA mfe prediction to verify the interaction via -compensatory mutation experiments. - - +The Compensatory Point Mutation Selector helps identify promising base pairs +from an IntaRNA MFE prediction for validation through compensatory mutation +experiments. See its own dependency and usage documentation. diff --git a/python/bin/Makefile.am b/python/bin/Makefile.am index e2a46584..f10b522d 100644 --- a/python/bin/Makefile.am +++ b/python/bin/Makefile.am @@ -1,2 +1,4 @@ noinst_PYTHON = build.py + +EXTRA_DIST = README.md diff --git a/python/bin/README.md b/python/bin/README.md new file mode 100644 index 00000000..66c8b2dc --- /dev/null +++ b/python/bin/README.md @@ -0,0 +1,16 @@ +# Optional standalone IntaRNApvalue executable + +From the repository root, in a Python 3.10+ environment: + +```bash +python3 -m pip install ./python pyinstaller +python3 python/bin/build.py +python/bin/build/IntaRNApvalue --help +``` + +The builder uses the invoking interpreter, resolves source paths independently +of the working directory, and creates `python/bin/build/IntaRNApvalue` (with an +`.exe` suffix on Windows). Build metadata is removed automatically. An existing +executable with that name is replaced by PyInstaller. Build on each target OS; +no precompiled binaries are shipped. The resulting executable still needs the +native IntaRNA program beside it or on PATH. diff --git a/python/bin/build.py b/python/bin/build.py index 6ee071a3..c46ec090 100644 --- a/python/bin/build.py +++ b/python/bin/build.py @@ -1,25 +1,27 @@ #!/usr/bin/env python3 - # Copyright 2019 # Author: Fabio Gutmann -import os -import sys +"""Build IntaRNApvalue with PyInstaller without relying on the working directory.""" + +from pathlib import Path import subprocess -from shutil import rmtree, move +import sys +from tempfile import TemporaryDirectory -if __name__ == '__main__': - build_dir = os.path.abspath(os.path.join(__file__, os.pardir)) - main = os.path.abspath(os.path.join(__file__, os.pardir, os.pardir, 'intarnapvalue', '__main__.py')) +def main(): + directory = Path(__file__).resolve().parent + source = directory.parent + # Keep build metadata private; leave only the requested binary in bin/build. + with TemporaryDirectory(prefix='IntaRNApvalue-build-') as temporary: + subprocess.run([sys.executable, '-m', 'PyInstaller', '--noconfirm', '--onefile', + '--name', 'IntaRNApvalue', '--paths', str(source), + '--distpath', str(directory / 'build'), '--workpath', temporary, + '--specpath', temporary, str(source / 'intarnapvalue' / '__main__.py')], + check=True) + return 0 - subprocess.call( f'{sys.executable} -m PyInstaller -y -F "{main}"', shell=True ) - os.remove(os.path.join(build_dir, '__main__.spec')) - rmtree(os.path.join(build_dir, 'build')) - os.mkdir(os.path.join(build_dir, 'build')) - if os.name == 'posix': - move(os.path.join(build_dir, 'dist', '__main__'), os.path.join(build_dir, 'build', 'IntaRNApvalue')) - else: - move(os.path.join(build_dir, 'dist', '__main__.exe'), os.path.join(build_dir, 'build', 'IntaRNApvalue.exe')) - rmtree(os.path.join(build_dir, 'dist')) +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/python/copomus/IntaRNA.py b/python/copomus/IntaRNA.py index 4481a43c..8ad1d5d2 100644 --- a/python/copomus/IntaRNA.py +++ b/python/copomus/IntaRNA.py @@ -1,135 +1,81 @@ # Copyright 2019 # Author: Fabio Gutmann -import sys -from typing import Dict, Union, List -from subprocess import Popen, PIPE, run +"""Checked IntaRNA subprocess and CSV adapter used by CopomuS.""" + +import csv +import io +import math +import re +from intarna_tools import csv_header, csv_rows, find_program, run_command class IntaRNA: - """IntaRNA python3 wrapper""" - def __init__(self): - pass + """IntaRNA Python wrapper; no-interaction output is represented by False.""" @staticmethod - def version_tuple(versionstring: str): - """ - Decomposes a version string into a tuple of entries - :param versionstring: The string encoding of a version, i.e. dot-separated numbers - """ - return tuple(map(int, versionstring.split('.'))) + def version_tuple(versionstring): + """Read numeric release components, accepting development suffixes.""" + match = re.match(r'(\d+)\.(\d+)\.(\d+)', versionstring) + if not match: + raise ValueError(f'Unrecognized IntaRNA version: {versionstring}') + return tuple(map(int, match.groups())) @staticmethod - def check_binary(required_version: str): - """ - Checks weather or not the IntaRNA binary can be found in PATH and exits with return code 1 if not found. - :param required_version: The minimum required version for IntaRNA - """ - p = run('IntaRNA --version', shell=True, stdout=PIPE, stderr=PIPE, universal_newlines=True) - if p.returncode: - sys.stderr.write('Please add the IntaRNA binary to your PATH.\n') - sys.exit(1) - if IntaRNA.version_tuple(p.stdout.split('\n')[0].split()[1]) < IntaRNA.version_tuple(required_version): - sys.stderr.write(f'This tool requires IntaRNA version {required_version}, please upgrade\n') - sys.exit(1) + def check_binary(required_version): + """Check executable availability and the minimum supported release.""" + output = run_command([find_program('IntaRNA'), '--version']) + match = re.search(r'IntaRNA\s+(\d+\.\d+\.\d+\S*)', output) + if not match or IntaRNA.version_tuple(match[1]) < IntaRNA.version_tuple(required_version): + raise RuntimeError(f'This tool requires IntaRNA version {required_version} or newer') - def run(self, query: str, target: str, qidxpos0: int, tidxpos0: int, outcsvcols: str, threads: int, n: int = 1, - param_file: str = '', extra_params: list = '', raw_stdout: bool = False) -> \ - Union[Dict[str, any], List[Dict[str, any]], str]: - """ - Runs IntaRNA on a pair of sequences and returns a dictionary with the output - :param query: The query sequence - :param target: The target sequence - :param qidxpos0: The starting index for the query sequence - :param tidxpos0: The starting index for the target sequence - :param outcsvcols: The csvcols returned by IntaRNA - :param threads: Thread count used by IntaRNA - :param n: The number of suboptimal interactions returned - :param param_file: Optional path to a file with additional parameters - :param extra_params: Optional list of additional parameters - :param raw_stdout: Weather or not stdout should be returned in raw format, not as a dict - :return: Dictionary with column name as key and payload as value - """ - if not extra_params: - extra_params = [] - param_file = f'--parameterFile={param_file}' if param_file else '' - p = Popen(['IntaRNA', '-q', query, '-t', target, '--outMode=C', f'--outcsvcols={outcsvcols}', - f'--qIdxPos0={qidxpos0}', f'--tIdxPos0={tidxpos0}', f'--outNumber={n}', f'--threads={threads}', param_file] - + extra_params, stdout=PIPE, stderr=PIPE, universal_newlines=True) - - stdout, stderr = p.communicate() - if p.returncode: - sys.stderr.write(f'IntaRNA Error: {stdout} {stderr}\n') - sys.exit(p.returncode) - if not raw_stdout: - if n == 1: - d = self.output_to_dict(stdout) - if type(d) != bool: - return d[0] - else: - return d - return self.output_to_dict(stdout) - else: + def run(self, query, target, qidxpos0, tidxpos0, outcsvcols, threads, n=1, + param_file='', extra_params=None, raw_stdout=False): + """Return one typed CSV row, all suboptimal rows, or raw stdout.""" + command = [find_program('IntaRNA'), '-q', query, '-t', target, '--outMode=C', + f'--outCsvCols={outcsvcols}', f'--qIdxPos0={qidxpos0}', + f'--tIdxPos0={tidxpos0}', f'--outNumber={n}', f'--threads={threads}'] + if param_file: + command.append(f'--parameterFile={param_file}') + stdout = run_command(command + list(extra_params or [])) + if raw_stdout: return stdout + rows = self.output_to_dict(stdout) + return rows[0] if rows and n == 1 else rows @staticmethod - def output_to_dict(data: str) -> Union[List[Dict[str, any]], bool]: - """ - Transforms IntaRNA output into a dictionary - :param data: Raw IntaRNA output - :return: dictionary with colname as key and payload as value (casted into fitting type) - >>> d = IntaRNA.output_to_dict('start1;end1;id1;id2;E\\n34;42;target;query;-4.7\\n') - >>> d + def output_to_dict(data): + r"""Parse quoted CSV, preserving IDs and numeric zero; only NA/empty is missing. + + >>> IntaRNA.output_to_dict('start1;end1;id1;id2;E\n34;42;target;query;-4.7\n') [{'start1': 34, 'end1': 42, 'id1': 'target', 'id2': 'query', 'E': -4.7}] - >>> [type(x) for x in d[0].values()] - [, , , , ] """ - subopts = [] - - lines = data.rstrip().split('\n') - if len(lines) == 2: - headers, entries = lines - entries = [entries] - elif len(lines) > 2: - headers = lines[0] - entries = lines[1:] - else: - headers = lines[0] - entries = [''] - headers = headers.split(';') - - for values in entries: - d = {} - values = values.split(';') - if len(headers) != len(values): # no interaction between query and target - return False - - for i, h in enumerate(headers): + reader = csv.reader(io.StringIO(data), delimiter=';', strict=True) + header = csv_header(reader) + rows = [] + for fields in csv_rows(reader, header): + row = {} + for name, value in zip(header, fields): + if name in ('id1', 'id2', 'seq1', 'seq2', 'subseq1', 'subseq2'): + row[name] = value + continue + if value in ('', 'NA'): + row[name] = None + continue try: - value = int(values[i]) # try to cast to int + number = int(value) except ValueError: try: - value = float(values[i]) # try to cast to float + number = float(value) except ValueError: - value = values[i] # leave as string - d[h] = value - subopts.append(d) - - for subopt in subopts: - for k, v in subopt.items(): - if not v: - subopt[k] = None - - return subopts + number = value + if isinstance(number, float) and not math.isfinite(number): + raise ValueError(f'IntaRNA returned a nonfinite {name}') + row[name] = number + rows.append(row) + return rows or False @staticmethod - def to_fasta(sequences: Dict[str, str]) -> str: - """ - Transforms sequences into a FASTA string - :param sequences: A dictionary with sequence name as key and raw sequence as tuple - :return: A FASTA string containing all sequences - """ - fasta_str = '' - for name, seq in sequences.items(): - fasta_str += f'>{name}\n{seq}\n' - return fasta_str + def to_fasta(sequences): + """Format a mapping of record names to sequences.""" + return ''.join(f'>{name}\n{seq}\n' for name, seq in sequences.items()) diff --git a/python/copomus/Makefile.am b/python/copomus/Makefile.am index bd5cdd99..0dc9467e 100644 --- a/python/copomus/Makefile.am +++ b/python/copomus/Makefile.am @@ -11,3 +11,5 @@ dist_copomus_DATA = \ copomusdir = $(bindir)/copomus + +EXTRA_DIST = README.md CopomuS-logo.jpg diff --git a/python/copomus/README.md b/python/copomus/README.md index bc237923..87930a00 100644 --- a/python/copomus/README.md +++ b/python/copomus/README.md @@ -18,12 +18,12 @@ CopomuS is highly configurable for specific purposes, see below for an explanati - [Mutation Generators](#mutation-generators) ## Dependencies -- [IntaRNA](https://github.com/BackofenLab/IntaRNA) >= 3.1.2 (needs to be in PATH) -- [python](https://www.python.org/) >= 3.7.0 +- [IntaRNA](https://github.com/BackofenLab/IntaRNA) >= 3.1.2 (beside the installed scripts, in `src/bin`, or on PATH) +- [python](https://www.python.org/) >= 3.10 ## Usage ```console -./CopomoS.py -q -t +python3 python/CopomuS.py -q -t mutation rank qIndex tIndex bpWildtype bpMutated mfeCover_rank mfeCover_ww mfeCover_wm mfeCover_mw mfeCover_mm mfeCover_is_valid E_rank E_ww E_wm E_mw E_mm E_is_valid minDeltaE_rank minDeltaE_wm/ww minDeltaE_mw/ww minDeltaE_wm/mm minDeltaE_mw/mm minDeltaE_ww/mm minDeltaE_min U4A&A11U 1 4 11 UA AU 1 True False False True 1 1 -9.32 -5.43 -7.25 -8.46 1 3 3.89 2.07 3.03 1.21 -0.86 1.21 @@ -37,26 +37,46 @@ G15C&C2G 2 15 2 GC CG 1 True True False True 1 2 -9.32 -12.37 -8.12 -12.2 0 5 -3 | Flag | Value | Default | Description | | :------------------ |:---------------------- | :------ | :-------------------- | | -h, --help | | | Gives detailed command help. | -| -q, --query | sequence | | Query as a raw sequence. | -| -t, --target | sequence | | Target as a raw sequence. | -| -m, --measure | {E, ed1, ed2, maxED, mfeCover, mfeOverlap, spotProb, qAccessProfile, tAccessProfile, qEnergyProfile, tEnergyProfile} | E | The measure used to calculate the rank of a base pair. Can be used multiple times. In that case, base pairs are ranked in order of measurements specified. | +| -q, --query | sequence or FASTA | | Query as RNA/DNA or a single-record FASTA file. | +| -t, --target | sequence or FASTA | | Target as RNA/DNA or a single-record FASTA file. | +| -m, --measure | {E, minDeltaE, mfeCover, Eqi, Eti, EDqi, cEDqi} | mfeCover, E, minDeltaE | The measure used to calculate the rank of a base pair. Can be used multiple times. In that case, base pairs are ranked in order of measurements specified. | | --qIdxPos0 | integer | 1 | (Optional) starting index for the query. | | --tIdxPos0 | integer | 1 | (Optional) starting index for the target. | -| -c, --candidateSelection | {mfe, mfeSO} | mfe | (Optional) method by which base pair candidates are selected from the sequences. | -| -f, --candidateFilters | {GU, AU, CG, lp, lpMfe, he, heMfe} | None | (Optional) filter by which base pair candidates are filtered. Can be used multiple times to combine filters. | +| -s, --candidateSelection | {mfe, mfeSO} | mfe | (Optional) method by which base pair candidates are selected from the sequences. | +| -f, --candidateFilter | {GU, AU, CG, GC, lp, lpMfe, he, heMfe} | None | (Optional) filter by which base pair candidates are filtered. Can be used multiple times to combine filters. | | -g, --generator | {flip, any} | flip | (Optional) method used to generate mutated sequences by which base pairs are ranked. | -| --mutationEncoding | string (ex. G1C&U7G) | | (Optional) allows specific selection of candidates by specifying a mutation encoding. If this option is specified, -c, -f and -g will be ignored. | +| --mutationEncoding | string (ex. G1C&U7G) | | (Optional) allows specific selection of candidates by specifying a mutation encoding. If this option is specified, -s, -f and -g will be ignored. | | -o, --output | file name | STDOUT | (Optional) file to which the output should be written to. | | -d, --delimiter | character | "\t" (tab) | (Optional) delimiter used to separate csv output. | | -p, --parameterFile | file name | None | (optional) parameter file to be used in IntaRNA calls to further guide predictions. | -| -t, --threads | integer | 0 (all) | (Optional) thread count used by IntaRNA. | +| --threads | integer | 0 (all) | (Optional) thread count used by IntaRNA. | + +Input sequences are normalized to uppercase RNA. FASTA files must contain exactly +one nonempty record; wrapped lines are supported. Negative indexing skips zero, +as in IntaRNA; `--qIdxPos0=0` / `--tIdxPos0=0` explicitly use zero-based indexing. +Mutation encodings are query-first, e.g. `G1C&C2G`, and comma-separated for multiple +candidates. Positions must lie inside the sequences, wildtype bases must match, +and both bases must change. Quote encodings in the shell because they contain `&`. +`GC` and `CG` are aliases for the same filter. + +Output is tab-separated UTF-8 with LF line endings by default. Numeric zero is a +valid score; unavailable scores are `NA`. Errors go to stderr with nonzero status. +An output file is replaced only after all predictions and rankings succeed. +Temporary profiles are private to one invocation, shared between measures, and +cleaned up on success or failure. All candidates reuse one wildtype prediction. +For example, from the repository root: + +```bash +python3 python/CopomuS.py -q query.fa -t target.fa -m EDqi -m cEDqi -m Eqi -o mutations.tsv +python3 python/CopomuS.py -q GC -t GC --mutationEncoding 'G1C&C2G' -m E +``` #### Measures | Name | Description | | :---- | :-------------------------------------------| | E | (class) Interaction Energy between query and target. | -| minDeltaE | (sort) Energy difference between different combinations of mutation types (ww, wm, mw, mm). | +| minDeltaE | (sort) Minimum of wm−ww, mw−ww, wm−mm and mw−mm; larger is better, including zero/negative values. | | mfeCover | (class) Mutation position has to be in MFE. | | Eqi | (class) Energy profile on query at mutation position. | | Eti | (class) Energy profile on target at mutation position. | diff --git a/python/copomus/candidate_filters.py b/python/copomus/candidate_filters.py index 03899120..1a6dada9 100644 --- a/python/copomus/candidate_filters.py +++ b/python/copomus/candidate_filters.py @@ -2,7 +2,7 @@ # Author: Fabio Gutmann from typing import Tuple, List -from copomus.indexing import idx_to_array_index +from copomus.indexing import idx_to_array_index, add_sub_idx def _filter_bp_type(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos0: int, bp_type='') -> \ @@ -34,7 +34,7 @@ def _filter_bp_type(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, t return new_bps -def _neighbors_in_mfe(q_index: int, t_index: int, bps: List[Tuple[int, int]]) -> int: +def _neighbors_in_mfe(q_index: int, t_index: int, bps: List[Tuple[int, int]], qidxpos0=1, tidxpos0=1) -> int: """ Counts how many neighboring base pairs are in the MFE region :param q_index: Base pair index on the query @@ -48,7 +48,9 @@ def _neighbors_in_mfe(q_index: int, t_index: int, bps: List[Tuple[int, int]]) -> """ count = 0 for q_o, t_o in [(+1, -1), (-1, +1)]: # get neighbors by offset - if (q_index+q_o, t_index+t_o) in bps: + q_neighbor = q_index + q_o if qidxpos0 == 0 else add_sub_idx(q_index, q_o) + t_neighbor = t_index + t_o if tidxpos0 == 0 else add_sub_idx(t_index, t_o) + if (q_neighbor, t_neighbor) in bps: count += 1 return count @@ -99,9 +101,10 @@ def filter_gc(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos def filter_lp(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos0: int) -> List[Tuple[int, int]]: """Filters lonely base pairs that can not stack (both neighbors not in MFE and both neighbors cant pair)""" new_bps = [] + bp_set = set(bps) for q_index, t_index in bps: - if not _neighbors_in_mfe(q_index, t_index, bps) and \ + if not _neighbors_in_mfe(q_index, t_index, bp_set, qidxpos0, tidxpos0) and \ not _neighbors_can_pair(q_index, t_index, q, t, qidxpos0, tidxpos0): continue new_bps.append((q_index, t_index)) @@ -111,9 +114,10 @@ def filter_lp(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos def filter_lp_mfe(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos0: int) -> List[Tuple[int, int]]: """Filters lonely base pairs in MFE interaction (both neighbors not in MFE)""" new_bps = [] + bp_set = set(bps) for q_index, t_index in bps: - if not _neighbors_in_mfe(q_index, t_index, bps): + if not _neighbors_in_mfe(q_index, t_index, bp_set, qidxpos0, tidxpos0): continue new_bps.append((q_index, t_index)) return new_bps @@ -133,9 +137,10 @@ def filter_he(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos def filter_he_mfe(bps: List[Tuple[int, int]], q: str, t: str, qidxpos0: int, tidxpos0: int) -> List[Tuple[int, int]]: """Filters helix ends (One neighbor in MFE, the other one is not)""" new_bps = [] + bp_set = set(bps) for q_index, t_index in bps: - if _neighbors_in_mfe(q_index, t_index, bps) == 1: + if _neighbors_in_mfe(q_index, t_index, bp_set, qidxpos0, tidxpos0) == 1: continue new_bps.append((q_index, t_index)) return new_bps @@ -146,6 +151,7 @@ def get_filter(string: str) -> any: 'GU': filter_gu, 'AU': filter_au, 'GC': filter_gc, + 'CG': filter_gc, 'lp': filter_lp, 'lpMfe': filter_lp_mfe, 'he': filter_he, diff --git a/python/copomus/candidate_selectors.py b/python/copomus/candidate_selectors.py index fbcb8e1e..2f3eeca7 100644 --- a/python/copomus/candidate_selectors.py +++ b/python/copomus/candidate_selectors.py @@ -1,10 +1,9 @@ # Copyright 2019 # Author: Fabio Gutmann -import sys from typing import List, Tuple from copomus.IntaRNA import IntaRNA -from copomus.indexing import idx_to_array_index as idx_aidx +from copomus.mutation import parse_mutation_encoding def get_mfe_bps(query: str, target: str, qidxpos0: int, tidxpos0: int, param_file: str, threads: int, **kwargs) -> \ @@ -23,8 +22,7 @@ def get_mfe_bps(query: str, target: str, qidxpos0: int, tidxpos0: int, param_fil i = IntaRNA() data = i.run(query, target, qidxpos0, tidxpos0, 'bpList', threads, param_file=param_file) if not data or not data['bpList']: - print('No favorable interaction between the specified sequences!') - sys.exit(0) + return [] for t in data['bpList'].strip().split(':'): t_index, q_index = [int(x) for x in t.strip('(').strip(')').split(',')] candidates.append((q_index, t_index)) # note we do query first, then target @@ -47,46 +45,23 @@ def get_mfe_bps_so(query: str, target: str, qidxpos0: int, tidxpos0: int, param_ i = IntaRNA() data = i.run(query, target, qidxpos0, tidxpos0, 'start1,end1,start2,end2', threads, param_file=param_file) if not data: - print('No favorable interaction between the specified sequences!') - sys.exit(0) + return [] data = i.run(query, target, qidxpos0, tidxpos0, 'bpList', threads, 1000, param_file, ['--qRegion', f"{data['start2']}-{data['end2']}", '--tRegion', f"{data['start1']}-{data['end1']}"]) - for subopt in data: + seen = set() + for subopt in data or []: for t in subopt['bpList'].strip().split(':'): t_index, q_index = [int(x) for x in t.strip('(').strip(')').split(',')] - if (q_index, t_index) not in candidates: + if (q_index, t_index) not in seen: + seen.add((q_index, t_index)) candidates.append((q_index, t_index)) # note we do query first, then target return candidates -def get_mutation_enc(query: str, target: str, qidxpos0: int, tidxpos0: int, param_file: str, **kwargs) \ - -> List[Tuple[int, int]]: - """ - Only selects the base pairs specified by the mutation encoding - :param query: The query sequence - :param target: The target sequence - :param qidxpos0: The counting start index for the query - :param tidxpos0: The counting start index for the target - :param param_file: Parameter file used for IntaRNA call - :return: List of tuples of indexes in form (query_index, target_index) with each tuple representing a base pair - """ - try: - candidates = [] - enc = kwargs['enc'] - for mut in enc.split(','): # split block mutations into single mutations - q, t = mut.split('&') - q_index = int(q[1:len(q) - 1]) # mutation indices - t_index = int(t[1:len(t) - 1]) - - if query[idx_aidx(q_index, qidxpos0)] != q[0] or target[idx_aidx(t_index, tidxpos0)] != t[0]: - print('The encoding you specified does not match with the sequences you specified!') - sys.exit(1) - candidates.append((q_index, t_index)) - return candidates - - except (IndexError, ValueError, KeyError): - print('Error parsing the encoding you specified!') - sys.exit(1) +def get_mutation_enc(query, target, qidxpos0, tidxpos0, param_file, threads=0, **kwargs): + """Validate explicit mutations and return their query/target coordinates.""" + return [(q, t) for q, t, _ in parse_mutation_encoding( + kwargs['enc'], query, target, qidxpos0, tidxpos0)] def get_selector(string: str) -> any: diff --git a/python/copomus/indexing.py b/python/copomus/indexing.py index 6b08be81..e3e48c27 100644 --- a/python/copomus/indexing.py +++ b/python/copomus/indexing.py @@ -15,7 +15,9 @@ def idx_to_array_index(idx: int, idxpos0: int): >>> idx_to_array_index(1, -1) 1 """ - return idx - idxpos0 - (1 if (idxpos0 < 0 < idx or idx < 0 < idxpos0) else 0) + if idxpos0 < 0 and idx == 0: + raise ValueError('index zero is skipped for negative sequence offsets') + return idx - idxpos0 - (1 if idxpos0 < 0 < idx else 0) def array_index_to_idx(i: int, idxpos0: int): @@ -31,7 +33,7 @@ def array_index_to_idx(i: int, idxpos0: int): >>> array_index_to_idx(1, -1) 1 """ - return idxpos0 + i + (1 if (idxpos0 < 0 < i or i < 0 < idxpos0) else 0) + return idxpos0 + i + (1 if idxpos0 < 0 <= idxpos0 + i else 0) def add_sub_idx(i1: int, change: int): diff --git a/python/copomus/measures.py b/python/copomus/measures.py index 05f0341f..133359b4 100644 --- a/python/copomus/measures.py +++ b/python/copomus/measures.py @@ -1,8 +1,7 @@ # Copyright 2019 # Author: Fabio Gutmann -import os -from tempfile import gettempdir +import math from copomus.mutation import Mutation from copomus.indexing import idx_to_array_index @@ -20,18 +19,23 @@ def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_f def post_compute(self): pass - def read_temp_file(self): - data = [] - for pipe in self.out_pipes: - file_name = f'CopomuS_{pipe}_{os.getpid()}.temp' - with open(os.path.join(gettempdir(), file_name), 'r') as f: - data.append(f.read()) - os.remove(os.path.join(gettempdir(), file_name)) - return data + def profile_value(self, pipe, index): + """Read a finite profile energy or NA, rejecting malformed/missing rows.""" + header_lines, delimiter, column = (2, '\t', 1) if pipe == 'qAcc' else (1, ';', 2) + try: + value = self.profiles[pipe].splitlines()[index + header_lines].split(delimiter)[column] + except (KeyError, IndexError) as error: + raise ValueError(f'IntaRNA returned an incomplete {pipe} profile') from error + if value == 'NA': + return None + number = float(value) + if not math.isfinite(number): + raise ValueError(f'IntaRNA returned a nonfinite {pipe} profile energy') + return number # noinspection PyTypeChecker def is_valid(self, alpha: float, beta: float = None) -> bool: - beta = alpha if not beta else beta + beta = alpha if beta is None else beta res = (self.score['ww'] < 0 and self.score['mm'] < 0 and (self.score['mw'] >= 0 or @@ -92,14 +96,14 @@ def __init__(self): def post_compute(self): score = { - 'wm/ww': round(self.score['wm'] - self.score['ww'], 2) if self.score['wm'] and self.score['ww'] else '', - 'mw/ww': round(self.score['mw'] - self.score['ww'], 2) if self.score['mw'] and self.score['ww'] else '', - 'wm/mm': round(self.score['wm'] - self.score['mm'], 2) if self.score['wm'] and self.score['mm'] else '', - 'mw/mm': round(self.score['mw'] - self.score['mm'], 2) if self.score['mw'] and self.score['mm'] else '', - 'ww/mm': round(self.score['ww'] - self.score['mm'], 2) if self.score['ww'] and self.score['mm'] else '' + 'wm/ww': round(self.score['wm'] - self.score['ww'], 2) if self.score['wm'] is not None and self.score['ww'] is not None else None, + 'mw/ww': round(self.score['mw'] - self.score['ww'], 2) if self.score['mw'] is not None and self.score['ww'] is not None else None, + 'wm/mm': round(self.score['wm'] - self.score['mm'], 2) if self.score['wm'] is not None and self.score['mm'] is not None else None, + 'mw/mm': round(self.score['mw'] - self.score['mm'], 2) if self.score['mw'] is not None and self.score['mm'] is not None else None, + 'ww/mm': round(self.score['ww'] - self.score['mm'], 2) if self.score['ww'] is not None and self.score['mm'] is not None else None } self.score = score - entries = [v for k, v in self.score.items() if v != '' and k != 'ww/mm'] + entries = [v for k, v in self.score.items() if v is not None and k != 'ww/mm'] self.score['min'] = min(entries) if entries else None def __lt__(self, other): @@ -107,7 +111,7 @@ def __lt__(self, other): return False elif other.non_i(): return True - return self.score['min'] > 0 > other.score['min'] or 0 < other.score['min'] < self.score['min'] + return self.score['min'] > other.score['min'] def __eq__(self, other): if self.non_i() and not other.non_i() or not self.non_i() and other.non_i(): @@ -136,6 +140,7 @@ def in_mfe(self): def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_file: str = ''): if not data: + self.score[mut_type] = False return mfe_cover = (m.target_n in range(data['start1'], data['end1']+1)) and \ (m.query_n in range(data['start2'], data['end2']+1)) @@ -174,14 +179,8 @@ def post_compute(self): self.score['is_valid'] = int(self.is_valid(self.alpha, self.beta)) def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_file: str = ''): - data = self.read_temp_file()[0] - data = data.strip().split('\n') - value = data[idx_to_array_index(m.query_n, m.qidxpos0) + 1].split(';')[2] - if value == 'NA': - self.score[mut_type] = None - else: - # noinspection PyTypeChecker - self.score[mut_type] = float(value) + self.score[mut_type] = self.profile_value('qMinE', idx_to_array_index(m.query_n, m.qidxpos0)) + def __lt__(self, other): if self.non_i(): @@ -216,14 +215,8 @@ def post_compute(self): # noinspection PyTypeChecker def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_file: str = ''): - data = self.read_temp_file()[0] - data = data.strip().split('\n') - value = data[idx_to_array_index(m.target_n, m.tidxpos0) + 1].split(';')[2] - if value == 'NA': - self.score[mut_type] = None - else: - # noinspection PyTypeChecker - self.score[mut_type] = float(value) + self.score[mut_type] = self.profile_value('tMinE', idx_to_array_index(m.target_n, m.tidxpos0)) + def __lt__(self, other): if self.non_i(): @@ -255,12 +248,9 @@ def post_compute(self): self.score['is_valid'] = int(self.is_valid(0)) def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_file: str = ''): - data = self.read_temp_file()[0] - data = data.strip().split('\n') - # noinspection PyTypeChecker - self.score[mut_type] = float(data[idx_to_array_index(m.query_n, m.qidxpos0) + 2].split('\t')[1]) + self.score[mut_type] = self.profile_value('qAcc', idx_to_array_index(m.query_n, m.qidxpos0)) + - # noinspection PyTypeChecker def is_valid(self, alpha: float, beta: float = None, gamma: float = 0) -> bool: return self.score['ww'] < self.score['mm'] @@ -278,8 +268,9 @@ def __eq__(self, other): return False if self.non_i() and other.non_i(): return True - return self.is_valid(0) == other.is_valid(0) and \ - self.score['mm'] - self.score['ww'] == other.score['mm'] - other.score['ww'] + return (self.is_valid(0) == other.is_valid(0) and + (not self.is_valid(0) or + self.score['mm'] - self.score['ww'] == other.score['mm'] - other.score['ww'])) class ConstantAccessProfileQuery(Measure): @@ -297,12 +288,9 @@ def post_compute(self): self.score['is_valid'] = int(self.is_valid(1)) def calculate_score(self, data: dict, mut_type: str, m: Mutation = None, param_file: str = ''): - data = self.read_temp_file()[0] - data = data.strip().split('\n') - # noinspection PyTypeChecker - self.score[mut_type] = float(data[idx_to_array_index(m.query_n, m.qidxpos0) + 2].split('\t')[1]) + self.score[mut_type] = self.profile_value('qAcc', idx_to_array_index(m.query_n, m.qidxpos0)) + - # noinspection PyTypeChecker def is_valid(self, alpha: float, beta: float = None, gamma: float = 0) -> bool: return abs(self.score['ww'] - self.score['mm']) < alpha diff --git a/python/copomus/mutation.py b/python/copomus/mutation.py index b33d5372..a9e24b45 100644 --- a/python/copomus/mutation.py +++ b/python/copomus/mutation.py @@ -2,13 +2,14 @@ # Author: Fabio Gutmann from dataclasses import dataclass, field -from typing import List +from typing import Any +import re from copomus.indexing import idx_to_array_index @dataclass() class Mutation: - measures: List[any] + measures: list[Any] qw: str # wildtype sequences tw: str qidxpos0: int # indexing start of query @@ -24,6 +25,10 @@ class Mutation: def __post_init__(self): mut_pos_q = idx_to_array_index(self.query_n, self.qidxpos0) mut_pos_t = idx_to_array_index(self.target_n, self.tidxpos0) + if not 0 <= mut_pos_q < len(self.qw) or not 0 <= mut_pos_t < len(self.tw): + raise ValueError('mutation position is outside its sequence') + if not re.fullmatch('[ACGU]{2}', self.bp_mm): + raise ValueError('mutated bases must be A, C, G or U') self.bp_ww = f'{self.qw[mut_pos_q]}{self.tw[mut_pos_t]}' self.qm = self.qw[0:mut_pos_q] + self.bp_mm[0] + self.qw[mut_pos_q+1:] self.tm = self.tw[0:mut_pos_t] + self.bp_mm[1] + self.tw[mut_pos_t+1:] @@ -33,3 +38,23 @@ def __lt__(self, other): def __str__(self): return f'{self.bp_ww[0]}{self.query_n}{self.bp_mm[0]}&{self.bp_ww[1]}{self.target_n}{self.bp_mm[1]}' + + +def parse_mutation_encoding(encoding, query, target, qidxpos0, tidxpos0): + """Validate query-first substitutions without allowing negative array wrapping.""" + mutations = [] + for entry in encoding.split(','): + match = re.fullmatch(r'([ACGU])(-?\d+)([ACGU])&([ACGU])(-?\d+)([ACGU])', entry) + if match is None: + raise ValueError(f'invalid mutation encoding: {entry!r}') + qw, q, qm, tw, t, tm = match.groups() + q, t = int(q), int(t) + qi, ti = idx_to_array_index(q, qidxpos0), idx_to_array_index(t, tidxpos0) + if not 0 <= qi < len(query) or not 0 <= ti < len(target): + raise ValueError('mutation position is outside its sequence') + if query[qi] != qw or target[ti] != tw: + raise ValueError('mutation encoding does not match the wildtype sequences') + if qm == qw or tm == tw: + raise ValueError('both encoded bases must change') + mutations.append((q, t, qm + tm)) + return mutations diff --git a/python/copomus/mutation_generators.py b/python/copomus/mutation_generators.py index 0c97d526..356e51e2 100644 --- a/python/copomus/mutation_generators.py +++ b/python/copomus/mutation_generators.py @@ -4,7 +4,7 @@ from copy import deepcopy from copomus.indexing import idx_to_array_index from copomus.measures import get_measure -from copomus.mutation import Mutation +from copomus.mutation import Mutation, parse_mutation_encoding from typing import List, Tuple @@ -66,13 +66,7 @@ def gen_specific(query: str, target: str, qidxpos0: int, tidxpos0: int, mms: Lis for m in mms: mm.append(get_measure(m, alpha, beta)) - for mut in encoding.split(','): # split block mutations into single mutations - q, t = mut.split('&') - q_index = int(q[1:len(q)-1]) # mutation indices - t_index = int(t[1:len(t)-1]) - - bp_mm = f'{q[len(q)-1]}{t[len(t)-1]}' - + for q_index, t_index, bp_mm in parse_mutation_encoding(encoding, query, target, qidxpos0, tidxpos0): m = Mutation(deepcopy(mm), query, target, qidxpos0, tidxpos0, q_index, t_index, bp_mm) mutations.append(m) return mutations diff --git a/python/examples/regions.csv b/python/examples/regions.csv new file mode 100644 index 00000000..5da6d985 --- /dev/null +++ b/python/examples/regions.csv @@ -0,0 +1,6 @@ +id1;start1;end1;id2;start2;end2 +b0001;266;273;query;116;123 +b0002;204;231;query;85;111 +b0003;229;262;query;96;125 +b0004;265;300;query;10;38 +b0005;281;295;query;5;22 diff --git a/python/examples/regions.params b/python/examples/regions.params new file mode 100644 index 00000000..a56400af --- /dev/null +++ b/python/examples/regions.params @@ -0,0 +1,9 @@ +# Plot query regions in CSV order. +id = id2 +start = start2 +end = end2 +title = Query interaction regions +xvline = 50 +xmin = 1 +xmax = 130 +rows = 200 diff --git a/python/intarna_tools.py b/python/intarna_tools.py new file mode 100644 index 00000000..6f3b9cfd --- /dev/null +++ b/python/intarna_tools.py @@ -0,0 +1,121 @@ +"""Shared, dependency-free I/O for the IntaRNA Python companion tools.""" + +import argparse +from contextlib import contextmanager +from pathlib import Path +import os +import re +import shutil +import subprocess +import sys +import tempfile + + +def positive_integer(value): + """argparse type for counts that must be positive.""" + number = int(value) + if number < 1: + raise argparse.ArgumentTypeError('expected a positive integer') + return number + + +def nonnegative_integer(value): + """argparse type for thread counts (zero means automatic).""" + number = int(value) + if number < 0: + raise argparse.ArgumentTypeError('expected a nonnegative integer') + return number + + +def normalize_sequence(value, *, allow_n=False): + """Validate a nonempty literal RNA/DNA sequence and normalize to RNA.""" + alphabet = 'ACGUTN' if allow_n else 'ACGUT' + if not re.fullmatch(f'[{alphabet}{alphabet.lower()}]+', value): + raise ValueError(f'expected a nonempty sequence containing only {alphabet}') + return value.upper().replace('T', 'U') + + +def read_fasta(filename): + """Read exactly one FASTA record, including wrapped sequence lines.""" + with open(filename, encoding='utf-8-sig') as handle: + lines = [line.strip() for line in handle if line.strip()] + if not lines or not lines[0].startswith('>'): + raise ValueError(f'{filename}: expected a FASTA header') + if any(line.startswith('>') for line in lines[1:]): + raise ValueError(f'{filename}: expected exactly one FASTA record') + return normalize_sequence(''.join(lines[1:])) + + +def read_sequence(value): + """Accept a literal RNA/DNA sequence or a single-record FASTA file.""" + # Test sequences first: long literals can exceed the filesystem name limit. + if re.fullmatch('[ACGUTacgut]+', value): + return normalize_sequence(value) + return read_fasta(value) if Path(value).is_file() else normalize_sequence(value) + + +def find_program(name, script_path=__file__): + """Prefer executable siblings, then the source-tree IntaRNA, then PATH.""" + directory = Path(sys.executable if getattr(sys, "frozen", False) else script_path).resolve().parent + candidates = [directory / name] + if name == 'IntaRNA': + candidates.append(directory.parent / 'src' / 'bin' / name) + if script_path == __file__ and not getattr(sys, 'frozen', False): + # pip installs the shared module in site-packages, separate from bin. + candidates.append(Path(sys.argv[0]).absolute().parent / name) + for candidate in candidates: + found = shutil.which(str(candidate)) + if found: + return found + found = shutil.which(name) + if not found: + raise RuntimeError(f'{name} executable not found beside the script or on PATH') + return found + + +def run_command(command, input_text=None): + """Run an argument vector without a shell; report command failures once.""" + result = subprocess.run(command, input=input_text, text=True, + stdout=subprocess.PIPE, stderr=subprocess.PIPE) + if result.returncode: + detail = (result.stderr or result.stdout).strip() + raise RuntimeError(f'{Path(command[0]).name} exited with status {result.returncode}: {detail}') + if result.stderr: + sys.stderr.write(result.stderr) + return result.stdout + + +def csv_header(reader): + """Read a nonempty CSV header with distinct, nonempty names.""" + header = next(reader, []) + if not header or not all(header) or len(set(header)) != len(header): + raise ValueError('CSV must have a nonempty header with unique column names') + return header + + +def csv_rows(reader, header): + """Yield nonblank rows, rejecting truncated or extra fields.""" + for row in reader: + if not row: + continue + if len(row) != len(header): + raise ValueError(f'CSV line {reader.line_num} has the wrong number of fields') + yield row + + +@contextmanager +def atomic_text_output(destination): + """Publish UTF-8/LF text only on success, preserving an existing file's mode.""" + destination = Path(destination).resolve() + temporary = None + try: + with tempfile.NamedTemporaryFile(mode='w', newline='', encoding='utf-8', + dir=destination.parent, delete=False) as handle: + temporary = Path(handle.name) + yield handle + if destination.exists(): + shutil.copymode(destination, temporary) + os.replace(temporary, destination) + finally: + if temporary is not None: + temporary.unlink(missing_ok=True) diff --git a/python/intarnapvalue/Makefile.am b/python/intarnapvalue/Makefile.am index aef8d658..a2749a4c 100644 --- a/python/intarnapvalue/Makefile.am +++ b/python/intarnapvalue/Makefile.am @@ -5,3 +5,5 @@ noinst_PYTHON = \ dinucl_shuffle.py \ intarna_pvalue.py \ test.py + +EXTRA_DIST = README.md diff --git a/python/intarnapvalue/README.md b/python/intarnapvalue/README.md index 7c67eb11..1fd00b0f 100644 --- a/python/intarnapvalue/README.md +++ b/python/intarnapvalue/README.md @@ -1,86 +1,111 @@ # IntaRNApvalue -A tool for calculating p-values of IntaRNA energy scores. - -### How it works: -This tool shuffles one or more of the given sequences (depending on shuffle mode) while preserving mono- and dinucleotide frequences. -Thus the resulting sequences are very similar with similar properties to the original one. -It then calculates the IntaRNA interaction scores for these newly generated sequences and uses them to approximate the p-value of the original target/query combination. -This p-value is a score for the likelihood that a sequence with a better interaction score than the original ones can randomly occur. -It can thus be used as a measurement for how good an interaction between two sequences is. - -## Dependencies: -- IntaRNA -##### And if you want to run it with python or compile it yourself: -- Python 3 (>= 3.6) -- numpy -- scipy - -## Installation: -You can either run this tool directly with python, install it with setuptools or run it as a compiled binary. -##### Run with python directly: -You don't have to install it, if you just plan on running it from a certain directory. -Simply copy the "intarnapvalue" folder where you want to run it from. -You have to get the dependencies yourself in this case. - -To use the tool from anywhere however, you need to install it as a python module. -All dependencies (except IntaRNA) will be installed automatically. -This way it can also be installed in a virtual environment. -```console -python setup.py install -``` -##### Run as binary: -Go into the bin directory. You can find pre-compiled builds for linux and windows for x64 directly. -You don't need to install any dependencies for running these, on Windows however you will need Microsoft Visual C++ Redistributable. -If you want to compile it yourself, you need to get the dependencies and PyInstaller. -I had problems with numpy 1.17.0, so I had to use 1.16.4. -Simply run the build.py with the python binary that has the dependencies and PyInstaller installed. -You will find your binary in the build folder. - -## Usage: -Run it as a module without installing (you need to be in the parent dir of intarnapvalue): -```console -python3 -m intarnapvalue -``` +Estimate the significance of an observed IntaRNA energy using a background of +**dinucleotide-preserving shuffled sequences**. This differs from +`IntaRNA_CSV_p-value.py`, which fits the energies of a genome-wide target screen. + +## Installation and invocation -Run it as a module from anywhere with python (need to run setup.py first, see [Installation](#installation)): -```console -python3 -m intarnapvalue +Use Python **3.10+**, IntaRNA on PATH (or beside the companion scripts / in the +checkout's `src/bin`), NumPy >= 1.23 and SciPy >= 1.9 for fitted distributions. +Scores-only and empirical modes need only the standard library and IntaRNA. +From the repository root: + +```bash +python3 -m pip install ./python +python3 -m intarnapvalue --help +# Alternatively, without installing the package: +PYTHONPATH=python python3 -m intarnapvalue --help ``` -Or use as a compiled binary: -```console -./IntaRNApvalue + +`pip` installs NumPy and SciPy automatically. IntaRNA itself must be installed +separately. `--help` works without IntaRNA or scientific Python dependencies. +For an optional standalone executable, install PyInstaller and run +`python3 python/bin/build.py`; see `python/bin/README.md`. + +## Sampling and interpretation + +Each shuffle preserves length, endpoints, and mono/dinucleotide counts. Lowercase +and DNA input are normalized to uppercase RNA. FASTA inputs may wrap across lines +but must contain exactly one nonempty record; multiple records are rejected. +In `b` mode, each target and each query are independently shuffled, and IntaRNA +compares corresponding pairs in bounded batches, not all combinations. + +`--samples` is the number of **reported interactions** collected. Pairs with no +reported interaction under the selected constraints are discarded and replaced. +The resulting p-value is therefore **conditional on a reported interaction**; it +does not estimate the unconditional probability over all shuffled pairs. The +number of discarded pairs is reported on stderr. The default attempt limit is +100 times `--samples`; reaching it is an error, with no partial score output. +Use `--max-attempts` to change that limit. An IntaRNA command failure stops +immediately. A parameter file is applied equally to the observed and shuffled +predictions; output settings enforce one energy per pair. + +The fitted modes retain the historical models on **raw energies**: SciPy +`genextreme` (GEV), `gumbel_l` (left-tail Gumbel), or `norm` (Gaussian). The p-value +is the fitted CDF at the observed energy, replacing numerical PDF integration. +This preserves the former model, including its GEV orientation; the separate +CSV tool fits **negated** energies and has different statistical semantics. +At least three distinct finite energies are required for fitting. This is a +numerical minimum, not assurance of statistical adequacy. Inspect a sufficiently +large background for the intended application. Degenerate fits are errors. + +`-d empirical` (alias `none`) reports the fraction of sampled energies less than +or equal to the observed energy, including ties. It accepts constant samples and +may return zero; its resolution is `1 / samples`, with no pseudocount correction. +Zero interaction energies remain valid. If the original pair has no reported +interaction, p-value mode fails; scores mode can still sample the background. + +## Arguments + +| Option | Default | Meaning | +| --- | --- | --- | +| `-q`, `--query` | Required | Literal RNA/DNA or a single-record FASTA path | +| `-t`, `--target` | Required | Literal RNA/DNA or a single-record FASTA path | +| `-s`, `--samples` | Required | Positive number of favorable shuffled pairs | +| `-m`, `--shuffle-mode` | Required | `q` (query), `t` (target), or `b` (both) | +| `-d`, `--distribution` | `gev` | `gev`, `gumbel`, `gauss`, `empirical`, or `none` | +| `-o`, `--output` | `pvalue` | `pvalue` or `scores` (one energy per line) | +| `-p`, `--parameterFile` | None | IntaRNA parameters used for every prediction | +| `--threads` | `0` | Nonnegative IntaRNA thread count; 0 uses available CPUs | +| `--randSeed` | None | Seed string; repeat with the same inputs/options/software | +| `--max-attempts` | `100 * samples` | Maximum pairs tried; must be at least samples | + +```bash +python3 -m intarnapvalue \ + -q GCUGAAAAACAUAACCCAUAAAAUGCUAGCUGUACCAGGAACCA \ + -t GGUUUCUUCGCCUCUGCGUUCACCAAAGUGUUCACCC \ + -s 1000 -m b --threads 1 --randSeed example -d empirical +python3 -m intarnapvalue -q query.fasta -t target.fasta -s 1000 -m t -o scores > background.txt ``` -Or import it and use it from within other python code: +A seed uses an instance-local RNG and does not reseed Python's global generator. +Corrected both-strand sampling changes results from earlier releases. Batch +results are restored to input order so thread scheduling does not reorder scores. +Diagnostics use stderr; successful stdout contains only the requested numbers. + +The internal shuffle utility also remains directly runnable: +`python3 python/intarnapvalue/dinucl_shuffle.py ACGU --samples 10 --randSeed example` +from the repository root (default: 4000 shuffles). + +## Python API and tests + +Passing an argument list configures an instance without running predictions: + ```python from intarnapvalue.intarna_pvalue import IntaRNApvalue -IntaRNApvalue(['--flag1', 'arg1']) -``` -## Arguments: - -| Flag | Value | Default | Description | -| ------------------ |:---------------------- | :------ | -------------------- | -| -h, --help | | | Gives detailed command help. | -| -q, --query | Sequence or FASTA file | | Query as a raw sequence or a file in FASTA format. Takes the first sequence from the file. | -| -t, --target | Sequence or FASTA file | | Target as a raw sequence or a file in FASTA format. Takes the first sequence from the file. | -| -c, --cardinality | Integer | | How many sequence pairs are randomly permuted and considered for p-value calculation. | -| -m, --shuffle-mode | {q, t, b} | | Which sequence will be shuffled: Query, Target or both. | -| -d, --distribution | {gev, gumbel, gauss} | gev | (optional) The distribution used for p-value calculation: Generalized Extreme Value Distribution, Gumbel Distribution or Gauss. | -| -o, --output | {pvalue, scores} | pvalue | (optional) If set to p-value, will only output p-value. If set to scores, will output every score from randomly generated sequences, but no p-value. | -| --threads | 0 - {max threads} | 0 | (optional) How many threads IntaRNA uses for score calculation. If set to 0 it will use all available threads. | -| --randSeed | Integer | None | (optional) The seed used for generating random sequences. | -| -p, --parameterFile | file name | None | (optional) parameter file to be used in IntaRNA calls to further guide predictions. | - -## Example -An example use-case could look like this with sequences in raw format: -```console -$ python3 -m intarnapvalue --query GCUGAAAAACAUAACCCAUAAAAUGCUAGCUGUACCAGGAACCA --target GGUUUCUUCGCCUCUGCGUUCACCAAAGUGUUCACCC --scores 10000 --shuffle-mode b --threads 0 -0.2421277140114205 +analysis = IntaRNApvalue(['-q', 'GGGGGGGG', '-t', 'CCCCCCCC', + '-s', '10', '-m', 't', '--randSeed', 'example']) +observed = analysis.get_single_score(analysis.query, analysis.target) # float or None +scores, discarded = analysis.get_scores() +if observed is not None: + pvalue = analysis.calculate_pvalue_empirical(observed, scores) ``` -Or you could specify any of the sequences as files: -```console -$ python3 -m intarnapvalue --query query.fasta --target target.fasta --scores 10000 --shuffle-mode b --threads 0 -0.2421277140114205 -``` \ No newline at end of file + +The Python API raises exceptions for failed predictions and invalid data. CLI +argument errors use exit code 2; runtime errors use exit code 1. +From the repository root, run `python3 -m unittest discover -s python/tests -v`. +Tests cover both/single-strand sampling, sample counts, bounded failures, +parameter forwarding, RNG isolation, FASTA validation, zero energies, shuffle +composition, and agreement between CDFs and the former numerical integration. diff --git a/python/intarnapvalue/__main__.py b/python/intarnapvalue/__main__.py index a4eeb080..609ccacd 100644 --- a/python/intarnapvalue/__main__.py +++ b/python/intarnapvalue/__main__.py @@ -1,8 +1,8 @@ #!/usr/bin/env python3 - # Copyright 2019 # Author: Fabio Gutmann -from intarnapvalue import intarna_pvalue +from intarnapvalue.intarna_pvalue import main -i = intarna_pvalue.IntaRNApvalue() +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/python/intarnapvalue/dinucl_shuffle.py b/python/intarnapvalue/dinucl_shuffle.py index ff892989..ca419d9a 100644 --- a/python/intarnapvalue/dinucl_shuffle.py +++ b/python/intarnapvalue/dinucl_shuffle.py @@ -9,8 +9,15 @@ # right fashion. +import argparse +from pathlib import Path import sys import random +from collections import deque +# Preserve direct invocation of this historical standalone utility. +if __name__ == '__main__' and not __package__: + sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from intarna_tools import normalize_sequence, positive_integer def compute_count_and_lists(s): @@ -50,13 +57,14 @@ def compute_count_and_lists(s): assert (nucl_total == len(s)) assert (dinucl_total == len(s)-1) return nucl_cnt, dinucl_cnt, nucl_dict - - -def choose_edge(x, dinucl_cnt): + + +def choose_edge(x, dinucl_cnt, rng=None): + rng = random if rng is None else rng num_in_list = 0 for y in ['A', 'C', 'G', 'U']: num_in_list += dinucl_cnt[x][y] - z = random.random() + z = rng.random() denom = dinucl_cnt[x]['A'] + dinucl_cnt[x]['C'] + dinucl_cnt[x]['G'] + dinucl_cnt[x]['U'] numerator = dinucl_cnt[x]['A'] if z < float(numerator)/float(denom): @@ -93,9 +101,9 @@ def connected_to_last(edge_list, nucl_list, last_ch): if x != last_ch and d[x] == 0: return 0 return 1 - -def eulerian(s): + +def eulerian(s, rng=None): nucl_cnt, dinucl_cnt, nucl_dict = compute_count_and_lists(s) # compute nucleotides appearing in s nucl_list = [] @@ -113,16 +121,17 @@ def eulerian(s): edge_list = [] for x in nucl_list: if x != last_ch: - edge_list.append([x, choose_edge(x, dinucl_cnt)]) + edge_list.append([x, choose_edge(x, dinucl_cnt, rng)]) ok = connected_to_last(edge_list, nucl_list, last_ch) return ok, edge_list, nucl_list, last_ch -def shuffle_edge_list(l): +def shuffle_edge_list(l, rng=None): + rng = random if rng is None else rng n = len(l) barrier = n for i in range(n-1): - z = int(random.random() * barrier) + z = int(rng.random() * barrier) tmp = l[z] l[z] = l[barrier - 1] l[barrier - 1] = tmp @@ -130,11 +139,15 @@ def shuffle_edge_list(l): return l -def dinucl_shuffle(s): +def dinucl_shuffle(s, rng=None): + """Preserve endpoints and mono/dinucleotide counts, using an optional RNG.""" + s = normalize_sequence(s) + if len(s) == 1: + return s ok = 0 edge_list, nucl_list = [], [] # initialize arrays while not ok: - ok, edge_list, nucl_list, last_ch = eulerian(s) + ok, edge_list, nucl_list, last_ch = eulerian(s, rng) nucl_cnt, dinucl_cnt, nucl_dict = compute_count_and_lists(s) # remove last edges from each vertex list, shuffle, then add back @@ -142,27 +155,39 @@ def dinucl_shuffle(s): for x, y in edge_list: nucl_dict[x].remove(y) for x in nucl_list: - shuffle_edge_list(nucl_dict[x]) + shuffle_edge_list(nucl_dict[x], rng) for x, y in edge_list: nucl_dict[x].append(y) # construct the eulerian path + nucl_dict = {base: deque(edges) for base, edges in nucl_dict.items()} l = [s[0]] prev_ch = s[0] for i in range(len(s)-2): - ch = nucl_dict[prev_ch][0] + ch = nucl_dict[prev_ch].popleft() l.append(ch) - del nucl_dict[prev_ch][0] prev_ch = ch l.append(s[-1]) t = ''.join(l) return t - + + +def main(argv=None): + """Print composition-preserving shuffles of one RNA/DNA sequence.""" + parser = argparse.ArgumentParser(description=main.__doc__) + parser.add_argument('sequence') + parser.add_argument('-s', '--samples', type=positive_integer, default=4000) + parser.add_argument('--randSeed', help='random seed string') + args = parser.parse_args(argv) + try: + sequence = normalize_sequence(args.sequence) + except ValueError as error: + parser.error(str(error)) + rng = random.Random(args.randSeed) + for _ in range(args.samples): + print(dinucl_shuffle(sequence, rng)) + return 0 + if __name__ == '__main__': - if len(sys.argv) != 2: - print("Usage: ./dinucl_shuffle.py ") - sys.exit(1) - seq = sys.argv[1] - for _ in range(4000): - print(dinucl_shuffle(seq)) + raise SystemExit(main()) diff --git a/python/intarnapvalue/intarna_pvalue.py b/python/intarnapvalue/intarna_pvalue.py index c09ae63c..a5120531 100644 --- a/python/intarnapvalue/intarna_pvalue.py +++ b/python/intarnapvalue/intarna_pvalue.py @@ -3,307 +3,236 @@ # Copyright 2019 # Author: Fabio Gutmann -from subprocess import PIPE, Popen, run -import random -import os +"""Estimate interaction-energy significance against dinucleotide shuffles.""" + import argparse +import csv +import io +import math +from pathlib import Path +import random import sys -from typing import List, Tuple -import numpy as np -from scipy.integrate import quad as integ -from scipy.stats import norm as gauss -from scipy.stats import genextreme as gev -from scipy.stats import gumbel_l as gum -from intarnapvalue import dinucl_shuffle +from tempfile import TemporaryDirectory + +from intarna_tools import (csv_header, csv_rows, find_program, nonnegative_integer, + positive_integer, read_fasta, read_sequence, run_command) +from intarnapvalue.dinucl_shuffle import dinucl_shuffle class IntaRNApvalue: + """Configure from argv; explicit argument lists do not run predictions. + + ``get_scores()`` returns favorable shuffled energies and the number of + rejected pairs. The fitted distribution is conditional on a reported + interaction under the supplied IntaRNA constraints. + """ + def __init__(self, test_args=None): - self.bin = self.find_binary() - self.query = '' - self.target = '' - self.samples = 0 - self.shuffle_query = False - self.shuffle_target = False - self.threads = '' - self.dist = '' - self.output = '' - self.parameterFile = '' self.process_cmd_args(test_args) - if not test_args: + self.bin = self.find_binary() + if test_args is None: self.main() - def main(self) -> None: - """The main function""" - if self.output == 'scores': # output scores and exit the process - scores, non_interactions = self.get_scores() - print('\n'.join(iter([str(x) for x in scores]))) - sys.exit(0) - - original_score = self.get_single_score(self.query, self.target) - if not original_score: # exit if given seq has no interaction and not scores output only mode - print('The query/target combination you specified has no favorable interaction') - sys.exit(1) - - scores, non_interactions = self.get_scores() - pvalue = 0.0 - if self.dist == 'gauss': - pvalue = self.calculate_pvalue_gauss(original_score, scores) - elif self.dist == 'none': - pvalue = self.calculate_pvalue_empirical(original_score, scores) - elif self.dist == 'gumbel': - pvalue = self.calculate_pvalue_gumbel(original_score, scores) - elif self.dist == 'gev': - pvalue = self.calculate_pvalue_gev(original_score, scores) - print(pvalue) - @staticmethod - def find_binary() -> str: - """Tries to find the IntaRNA executable and returns path to binary or exits with error code 1 if not found""" - if not run('IntaRNA --version', shell=True, stdout=PIPE, stderr=PIPE).returncode: - #print(run('IntaRNA --version', shell=True, stdout=PIPE).stdout) - return 'IntaRNA' - - # if binary not found in path, search in parent of this script recursively - bin_name = 'IntaRNA.exe' if os.name == 'nt' else 'IntaRNA' - for dir_path, dir_name, file_names in os.walk(os.path.abspath(os.path.join(os.curdir, '..'))): - if bin_name in file_names: - #print(run(os.path.join(dir_path, bin_name)+' --version', shell=True, stdout=PIPE).stdout) - return os.path.join(dir_path, bin_name) - - print('Error: Cannot find IntaRNA binary executable, please add it to your PATH') - sys.exit(1) - - def process_cmd_args(self, test_args=None) -> None: - """Processes all commandline args and sets them as instance variables. - - >>> i = IntaRNApvalue(['-q', 'AGGAUG', '-t', 'UUUAUCGUU', '-s', '10', '-m', 'b', '-d', 'gauss', - ... '--threads', '3']) - >>> i.query - 'AGGAUG' - >>> i.target - 'UUUAUCGUU' - >>> i.samples - 10 - >>> i.shuffle_query - True - >>> i.shuffle_target - True - >>> i.threads - '3' - >>> i.dist - 'gauss' - >>> i.output - 'pvalue' - >>> i = IntaRNApvalue(['-q', 'Z', '-t', 'UUUAUCGUU', '-s', '10', '-m', 'b', '--threads', '3']) - Traceback (most recent call last): - ... - SystemExit: 1 - >>> open('test.fasta', 'w').write('>someseq\\nGACU') - 13 - >>> i = IntaRNApvalue(['-q', 'test.fasta', '-t', 'UUUAUCGUU', '-s', '10', '-m', 'b', '--threads', '3']) - >>> i.query - 'GACU' - >>> os.remove('test.fasta') - """ - parser = argparse.ArgumentParser(description='Calculates p-values of IntaRNA energy scores.') - parser.add_argument('-q', '--query', dest='query', type=str, help='Query sequence', required=True) - parser.add_argument('-t', '--target', dest='target', type=str, help='Target sequence', required=True) - parser.add_argument('-s', '--samples', dest='samples', type=int, required=True, - help='How many sequence pairs are randomly permuted and considered for p-value calculation.') - parser.add_argument('-m', '--shuffle-mode', dest='shuffleMode', required=True, choices=['q', 't', 'b'], - help='Which sequences are going to be shuffled: both, query only or target only.') - parser.add_argument('-p', '--parameterFile', dest='parameterFile', default="", - help='Optional parameter file for IntaRNA provide further IntaRNA parameters and prediction constraints.') - parser.add_argument('-d', '--distribution', dest='dist', choices=['gev', 'gumbel', 'gauss'], default='gev', - help='[gev] Which distribution is fitted and used to calculate the p-value.') - parser.add_argument('-o', '--output', dest='output', choices=['pvalue', 'scores'], default='pvalue', - help='[pvalue] If set to "scores", outputs all IntaRNA scores from random sequences to STDOUT. ' - 'This is useful for piping the scores. Otherwise outputs just the p-value.') - parser.add_argument('--threads', type=str, default='0', help='[0] Sets the amount of threads used for IntaRNA.') - parser.add_argument('--randSeed', type=str, default=None, - help='Random seed to make sequence shuffling and thus output deterministic.') - + def find_binary(): + """Use the same executable lookup as the other companion tools.""" + return find_program('IntaRNA') + + def process_cmd_args(self, test_args=None): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('-q', '--query', required=True, help='RNA/DNA sequence or single-record FASTA') + parser.add_argument('-t', '--target', required=True, help='RNA/DNA sequence or single-record FASTA') + parser.add_argument('-s', '--samples', type=positive_integer, required=True, + help='number of favorable shuffled pairs to collect') + parser.add_argument('-m', '--shuffle-mode', dest='shuffle_mode', required=True, + choices=['q', 't', 'b'], help='shuffle query, target, or both independently') + parser.add_argument('-p', '--parameterFile', default='', help='IntaRNA parameter file for ALL predictions') + parser.add_argument('-d', '--distribution', dest='dist', + choices=['gev', 'gumbel', 'gauss', 'empirical', 'none'], default='gev', + help='fitted distribution; empirical (alias none) counts scores <= observed') + parser.add_argument('-o', '--output', choices=['pvalue', 'scores'], default='pvalue', + help='print one p-value or one shuffled energy per line') + parser.add_argument('--threads', type=nonnegative_integer, default=0, + help='IntaRNA threads; zero selects all available CPUs') + parser.add_argument('--randSeed', default=None, help='random seed string for reproducible shuffling') + parser.add_argument('--max-attempts', type=positive_integer, + help='maximum shuffled pairs tried (default: 100 times --samples)') args = parser.parse_args(test_args) - - allowed = ['G', 'A', 'C', 'U', 'T'] - if os.path.exists(args.query): - args.query = self.read_fasta_file(args.query) - - if os.path.exists(args.target): - args.target = self.read_fasta_file(args.target) - - if args.parameterFile and not os.path.exists(args.parameterFile): - print("Cannot find provided parameter file '"+args.parameterFile+"'") - sys.exit(1) - - if False in [n in allowed for n in args.query.upper()] or False in [n in allowed for n in args.target.upper()]: - print('A sequence you specified contains illegal characters, allowed: G, A, C, U (T)') - sys.exit(1) - - self.shuffle_query = True if args.shuffleMode in ['b', 'q'] else False - self.shuffle_target = True if args.shuffleMode in ['b', 't'] else False - for key, value in args.__dict__.items(): - self.__setattr__(key, value) - - random.seed(a=args.randSeed) - - @staticmethod - def read_fasta_file(filename: str) -> str: - """Reads a FASTA file and returns the sequence. Exits the program if something goes wrong. - - >>> with open('test.fasta', 'w') as f: - ... f.write('>somename\\nGACUGGAGUGC') - 21 - >>> IntaRNApvalue.read_fasta_file('test.fasta') - 'GACUGGAGUGC' - >>> IntaRNApvalue.read_fasta_file('test2.fasta') - Traceback (most recent call last): - ... - SystemExit: 1 - >>> os.remove('test.fasta') - """ try: - with open(filename) as f: - data = f.read().split('\n') - if not data[0].startswith('>'): - print('Error: {} is not in FASTA format!'.format(filename)) - sys.exit(1) - f.close() - return data[1].upper().replace('T', 'U') - except FileNotFoundError as e: # In theory this can never happen, lets prevent it anyways - print(e) - sys.exit(1) - except IndexError: - print('Error: {} is not in FASTA format!'.format(filename)) - sys.exit(1) + args.query = read_sequence(args.query) + args.target = read_sequence(args.target) + if args.parameterFile and not Path(args.parameterFile).is_file(): + raise ValueError(f'Cannot find provided parameter file: {args.parameterFile}') + args.max_attempts = args.max_attempts or 100 * args.samples + if args.max_attempts < args.samples: + raise ValueError('--max-attempts must be at least --samples') + except (OSError, ValueError) as error: + parser.error(str(error)) + for key, value in vars(args).items(): + setattr(self, key, value) + self.shuffle_query = args.shuffle_mode in ('q', 'b') + self.shuffle_target = args.shuffle_mode in ('t', 'b') + self.rng = random.Random(args.randSeed) + + read_fasta_file = staticmethod(read_fasta) @staticmethod - def shuffle_sequence(seq: str, n: int) -> List[str]: - """Shuffles a sequence n times and returns a list of sequences, duplicate entries are possible - - >>> random.seed('IntaRNA') - >>> IntaRNApvalue.shuffle_sequence('AGGAUGGGGGA', 5) - ['AUGGAGGGGGA', 'AUGGGGAGGGA', 'AGGGGAUGGGA', 'AUGGGGGAGGA', 'AUGGGAGGGGA'] - """ - return [dinucl_shuffle.dinucl_shuffle(seq.upper()) for _ in range(n)] + def shuffle_sequence(seq, n, rng=None): + """Generate n independent, composition-preserving RNA shuffles.""" + if n < 0: + raise ValueError('shuffle count cannot be negative') + return [dinucl_shuffle(seq, rng=rng) for _ in range(n)] @staticmethod - def to_fasta(sequences: List[str]) -> str: - """Combines a list of sequences into a string in FASTA format - - >>> IntaRNApvalue.to_fasta(['AUGGAGGGGGA', 'AUGGGGAGGGA', 'AGGGGAUGGGA', 'AUGGGGGAGGA', 'AUGGGAGGGGA']) - '>0\\nAUGGAGGGGGA\\n>1\\nAUGGGGAGGGA\\n>2\\nAGGGGAUGGGA\\n>3\\nAUGGGGGAGGA\\n>4\\nAUGGGAGGGGA\\n' - """ - fasta_str = '' - n = 0 - for seq in sequences: - fasta_str += '>{}\n{}\n'.format(n, seq) - n += 1 - return fasta_str - - def get_scores(self) -> Tuple[List[float], int]: - """Calculates n IntaRNA energy scores from random sequences with given parameters as class variables""" - scores = [] - missing = self.samples - non_interactions = 0 - paramFile = "" - if self.parameterFile : - paramFile = "--parameterFile="+self.parameterFile - - while missing > 0: - if self.shuffle_query and self.shuffle_target: # shuffle both - query = self.shuffle_sequence(self.query, 1)[0] # get a random query - target = 'STDIN' - shuffles = self.to_fasta(self.shuffle_sequence(self.query, missing)) - elif self.shuffle_query and not self.shuffle_target: # only shuffle query - query = 'STDIN' - target = self.target # target stays the same - shuffles = self.to_fasta(self.shuffle_sequence(self.query, missing)) - else: # only shuffle target - query = self.query # query stays the same - target = 'STDIN' - shuffles = self.to_fasta(self.shuffle_sequence(self.target, missing)) - - p = Popen([self.bin, '-q', query, '-t', target, '--outMode=C', '--outCsvCols=E', '--threads', self.threads, paramFile], - stdout=PIPE, stdin=PIPE, universal_newlines=True) - stdout, stderr = p.communicate(input=shuffles) # send shuffles as STDIN - if p.returncode: # If IntaRNA exits with a returncode != 0, skip this iteration - print("IntaRNA failed (ret={})".format(p.returncode)) - continue - stdout = stdout.split('\n') # split on newline - del stdout[0], stdout[-1] # remove first element aka 'E' and trailing newline element - scores.extend(stdout) # add elements to scores - missing = self.samples - len(scores) - non_interactions += missing # count non-interactions - - # return list with all elements as float and amount of non-interactions - return [float(x) for x in scores], non_interactions - - def get_single_score(self, query, target) -> float: - """Gets an IntaRNA score to a single query/target combination""" - o = run('{} -q {} -t {} --outMode=C --outCsvCols=E --threads {}'.format(self.bin, query, target, self.threads), - stdout=PIPE, stdin=PIPE, shell=True, universal_newlines=True).stdout - if o.startswith('E') and o != 'E\n': # Check that we got a result - return float(o.split('\n')[1]) - else: - return 0 # no interaction + def to_fasta(sequences): + """Number the records so batch results can be reordered deterministically.""" + return ''.join(f'>{index}\n{seq}\n' for index, seq in enumerate(sequences)) + + def _predict(self, query, target, input_text=None, *, pairwise=False): + """Apply identical constraints to observed and shuffled pairs; fail once.""" + command = [self.bin, '-q', query, '-t', target, '--outMode=C', + '--outCsvCols=id1,id2,E', '--outSep=;', '--outNumber=1', + '--outPerRegion=0', f'--outPairwise={int(pairwise)}', + f'--threads={self.threads}'] + if self.parameterFile: + command.append(f'--parameterFile={self.parameterFile}') + reader = csv.reader(io.StringIO(run_command(command, input_text)), delimiter=';', strict=True) + header = csv_header(reader) + if header != ['id1', 'id2', 'E']: + raise ValueError('IntaRNA returned an unexpected CSV header') + rows = [] + for target_id, query_id, value in csv_rows(reader, header): + score = float(value) + if not math.isfinite(score): + raise ValueError('IntaRNA returned a nonfinite energy') + rows.append((target_id, query_id, score)) + return rows + + def get_single_score(self, query, target): + """Return the observed energy, including zero, or None for no interaction.""" + rows = self._predict(query, target) + if len(rows) > 1: + raise ValueError('IntaRNA returned multiple energies for one pair') + return rows[0][2] if rows else None + + def get_scores(self): + """Collect favorable scores in bounded batches; reject invalid/duplicate IDs.""" + scores, attempts = [], 0 + # Two files allow IntaRNA's pairwise mode to pair independent query/target + # shuffles by index, avoiding a Cartesian product and repeated query draws. + with TemporaryDirectory(prefix='IntaRNApvalue-') as directory: + query_file = Path(directory) / 'query.fa' + target_file = Path(directory) / 'target.fa' + while len(scores) < self.samples and attempts < self.max_attempts: + count = min(1000, self.samples - len(scores), self.max_attempts - attempts) + if self.shuffle_query and self.shuffle_target: + query_file.write_text(self.to_fasta(self.shuffle_sequence(self.query, count, self.rng)), + encoding='utf-8') + target_file.write_text(self.to_fasta(self.shuffle_sequence(self.target, count, self.rng)), + encoding='utf-8') + rows = self._predict(str(query_file), str(target_file), pairwise=True) + elif self.shuffle_query: + fasta = self.to_fasta(self.shuffle_sequence(self.query, count, self.rng)) + rows = self._predict('STDIN', self.target, fasta) + else: + fasta = self.to_fasta(self.shuffle_sequence(self.target, count, self.rng)) + rows = self._predict(self.query, 'STDIN', fasta) + batch = {} + for target_id, query_id, score in rows: + identifier = query_id if self.shuffle_query else target_id + index = int(identifier) + if not 0 <= index < count or index in batch: + raise ValueError('IntaRNA returned an invalid or duplicate shuffled-pair ID') + if self.shuffle_query and self.shuffle_target and target_id != query_id: + raise ValueError('IntaRNA returned an unpaired shuffle result') + batch[index] = score + scores.extend(batch[index] for index in sorted(batch)) + attempts += count + if len(scores) != self.samples: + raise RuntimeError(f'Only {len(scores)} of {self.samples} favorable shuffled interactions ' + f'found in {attempts} attempts; adjust constraints or --max-attempts') + return scores, attempts - len(scores) @staticmethod - def calculate_pvalue_empirical(original_score, scores: list = None) -> float: - """Calculates a p-value to a target/query combination empirical with a given amount of shuffle iterations - - >>> i = IntaRNApvalue(['-q', 'AGGAUG', '-t', 'UUUAUCGUU', '--scores', '10', '-m', 'b', '--threads', '3']) - >>> i.calculate_pvalue_empirical(-10.0, [-1.235, -1.435645, -6.234234, -12.999, -15.23, -6.98, -6.23, -2.78]) - 0.25 - """ - return [score <= original_score for score in scores].count(True) / len(scores) + def _sample(original_score, scores, *, fit=False): + """Reject nonfinite, empty or degenerate fitting data before SciPy runs.""" + if not math.isfinite(original_score): + raise ValueError('observed energy must be finite') + values = list(scores) + if not values or not all(math.isfinite(x) for x in values): + raise ValueError('scores must contain finite energies') + if fit and (len(values) < 3 or len(set(values)) < 3): + raise ValueError('distribution fitting needs at least three distinct energies; use empirical otherwise') + return values @staticmethod - def calculate_pvalue_gauss(original_score, scores: list) -> float: - """Calculates a p-value to a target/query combination by int. with a given amount of shuffle iterations by - fitting a gaussian distribution and integrating from -inf to the original score - - >>> i = IntaRNApvalue(['-q', 'AGGAUG', '-t', 'UUUAUCGUU', '--scores', '10', '-m', 'b', '--threads', '3']) - >>> i.calculate_pvalue_gauss(-10.0, [-1.235, -1.435645, -6.234234, -12.999, -15.23, -6.98, -6.23, -2.78]) - 0.2429106747265256 - """ - loc, scale = gauss.fit(scores) - - def f(x): - return gauss.pdf(x, loc=loc, scale=scale) - return integ(f, -np.inf, original_score)[0] + def calculate_pvalue_empirical(original_score, scores): + """Fraction of favorable shuffled energies at least as low as observed.""" + values = IntaRNApvalue._sample(original_score, scores) + return sum(score <= original_score for score in values) / len(values) @staticmethod - def calculate_pvalue_gumbel(original_score: float, scores: list) -> float: - """Calculates a p-value to a target/query combination by int. with a given amount of shuffle iterations by - fitting a gumbel distribution and integrating from -inf to the original score - - >>> i = IntaRNApvalue(['-q', 'AGGAUG', '-t', 'UUUAUCGUU', '--scores', '10', '-m', 'b', '--threads', '3']) - >>> i.calculate_pvalue_gumbel(-10.0, [-1.235, -1.435645, -6.234234, -12.999, -15.23, -6.98, -6.23, -2.78]) - 0.19721934073203196 - """ - loc, scale = gum.fit(scores) - - def f(x): - return gum.pdf(x, loc=loc, scale=scale) - return integ(f, -np.inf, original_score)[0] + def _fitted_pvalue(original_score, scores, name): + from scipy import stats + values = IntaRNApvalue._sample(original_score, scores, fit=True) + distribution = getattr(stats, name) + # Keep the historical fits on raw energies (including the GEV shape). + # The genome-wide CSV tool fits negated energies: a different model. + params = distribution.fit(values) + if not all(math.isfinite(x) for x in params) or params[-1] <= 0: + raise ValueError('distribution fit returned invalid parameters') + if not all(math.isfinite(float(x)) for x in distribution.logpdf(values, *params)): + raise ValueError('distribution fit does not cover all sampled energies') + result = float(distribution.cdf(original_score, *params)) + if not math.isfinite(result) or not 0 <= result <= 1: + raise ValueError('distribution fit returned an invalid p-value') + return result @staticmethod - def calculate_pvalue_gev(original_score: float, scores: list) -> float: - """Calculates a p-value to a target/query combination by int. with a given amount of shuffle iterations by - fitting a generalized extreme value distribution and integrating from -inf to the original score + def calculate_pvalue_gauss(original_score, scores): + return IntaRNApvalue._fitted_pvalue(original_score, scores, 'norm') - >>> i = IntaRNApvalue(['-q', 'AGGAUG', '-t', 'UUUAUCGUU', '--scores', '10', '-m', 'b', '--threads', '3']) - >>> i.calculate_pvalue_gev(-10.0, [-1.235, -1.435645, -6.234234, -12.999, -15.23, -6.98, -6.23, -2.78]) - 0.17611816922560236 - """ - shape, loc, scale = gev.fit(scores) + @staticmethod + def calculate_pvalue_gumbel(original_score, scores): + return IntaRNApvalue._fitted_pvalue(original_score, scores, 'gumbel_l') - def f(x): - return gev.pdf(x, shape, loc=loc, scale=scale) - return integ(f, -np.inf, original_score)[0] + @staticmethod + def calculate_pvalue_gev(original_score, scores): + return IntaRNApvalue._fitted_pvalue(original_score, scores, 'genextreme') + + def main(self): + """Publish output only after every required prediction and fit succeeds.""" + original_score = None + if self.output != 'scores': + original_score = self.get_single_score(self.query, self.target) + if original_score is None: + raise ValueError('The query/target combination has no favorable interaction') + scores, non_interactions = self.get_scores() + if non_interactions: + print(f'Skipped {non_interactions} shuffled pairs without a reported interaction; ' + 'p-values are conditional on reported interactions.', file=sys.stderr) + if self.output == 'scores': + print('\n'.join(map(str, scores))) + else: + calculators = {'gauss': self.calculate_pvalue_gauss, 'gumbel': self.calculate_pvalue_gumbel, + 'gev': self.calculate_pvalue_gev, 'empirical': self.calculate_pvalue_empirical, + 'none': self.calculate_pvalue_empirical} + print(calculators[self.dist](original_score, scores)) + + +def main(argv=None): + """CLI entry point; --help works without IntaRNA, NumPy or SciPy installed.""" + try: + tool = IntaRNApvalue(sys.argv[1:] if argv is None else argv) + tool.main() + except ImportError as error: + print(f'Missing Python dependency: {error}. Install numpy and scipy.', file=sys.stderr) + return 1 + except (OSError, ValueError, RuntimeError, csv.Error) as error: + print(f'Error: {error}', file=sys.stderr) + return 1 + return 0 if __name__ == '__main__': - i = IntaRNApvalue() + raise SystemExit(main()) diff --git a/python/intarnapvalue/test.py b/python/intarnapvalue/test.py index 9bbea0fe..9e33e70e 100644 --- a/python/intarnapvalue/test.py +++ b/python/intarnapvalue/test.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 - # Copyright 2019 # Author: Fabio Gutmann -from intarnapvalue import intarna_pvalue +"""Compatibility entry point for the source-tree p-value regression tests.""" -if __name__ == '__main__': - import doctest - doctest.testmod(intarna_pvalue) +from pathlib import Path +import unittest +if __name__ == '__main__': + tests = Path(__file__).resolve().parents[1] / 'tests' + suite = unittest.defaultTestLoader.discover(str(tests), pattern='test_intarnapvalue.py') + raise SystemExit(not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful()) diff --git a/python/plotRegions.example.png b/python/plotRegions.example.png new file mode 100644 index 00000000..cdd0961b Binary files /dev/null and b/python/plotRegions.example.png differ diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 00000000..09977b5b --- /dev/null +++ b/python/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools>=61"] +build-backend = "setuptools.build_meta" diff --git a/python/requirements-output-tools.txt b/python/requirements-output-tools.txt new file mode 100644 index 00000000..7bc0ac84 --- /dev/null +++ b/python/requirements-output-tools.txt @@ -0,0 +1,5 @@ +# Optional scientific dependencies for IntaRNA Python companion tools. +# Requires Python >= 3.10. The C++ predictor does not need these packages. +numpy>=1.23 +scipy>=1.9 +matplotlib>=3.6 diff --git a/python/setup.py b/python/setup.py index 96ac0a86..695cff8a 100644 --- a/python/setup.py +++ b/python/setup.py @@ -1,24 +1,23 @@ #!/usr/bin/env python3 - +"""Install IntaRNA's companion modules and command-line tools.""" from setuptools import setup setup( name='intarnapvalue', version='0.1', - description='Calculates p-values of IntaRNA energy scores', + description='IntaRNA energy significance, mutation selection and output tools', author='Fabio Gutmann', - python_requires='>=3.6.0', + python_requires='>=3.10', url='https://github.com/BackofenLab/IntaRNA/', - packages=['intarnapvalue'], -# install_requires=['scipy', 'numpy'], - include_package_data=True, + packages=['intarnapvalue', 'copomus'], + py_modules=['intarna_tools'], + scripts=['CopomuS.py', 'IntaRNA-fuse.py', 'IntaRNA_CSV_p-value.py', 'IntaRNA_plotRegions.py'], + entry_points={'console_scripts': ['IntaRNApvalue=intarnapvalue.intarna_pvalue:main']}, + install_requires=['numpy>=1.23', 'scipy>=1.9'], + extras_require={'plot': ['matplotlib>=3.6']}, classifiers=[ - 'Programming Language :: Python', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.6', - 'Programming Language :: Python :: 3.7', - 'Programming Language :: Python :: Implementation :: CPython', - 'Programming Language :: Python :: Implementation :: PyPy' - ] + 'Programming Language :: Python :: 3 :: Only', + ], ) diff --git a/python/tests/test_common_tools.py b/python/tests/test_common_tools.py new file mode 100644 index 00000000..e1df35a0 --- /dev/null +++ b/python/tests/test_common_tools.py @@ -0,0 +1,100 @@ +"""Contracts shared by standalone and installed Python companion commands.""" + +import importlib.util +from pathlib import Path +import stat +import subprocess +import sys +import tempfile +import unittest +from unittest.mock import patch + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT)) +from intarna_tools import atomic_text_output, find_program, run_command + + +class CommonToolsTests(unittest.TestCase): + def test_atomic_publication_and_rollback(self): + with tempfile.TemporaryDirectory() as directory: + output = Path(directory) / 'result.csv' + output.write_text('old') + output.chmod(0o640) + with self.assertRaises(RuntimeError): + with atomic_text_output(output) as handle: + handle.write('partial') + raise RuntimeError('failed computation') + self.assertEqual(output.read_text(), 'old') + with atomic_text_output(output) as handle: + handle.write('α\n') + self.assertEqual(output.read_bytes(), 'α\n'.encode()) + self.assertEqual(stat.S_IMODE(output.stat().st_mode), 0o640) + self.assertEqual(list(Path(directory).iterdir()), [output]) + + def test_subprocess_arguments_stdin_and_failure(self): + payload = '$(touch SHOULD_NOT_EXIST); quoted spaces' + result = run_command([sys.executable, '-c', 'import sys; print(sys.argv[1]); print(sys.stdin.read())', payload], 'ACGU') + self.assertEqual(result, payload + '\nACGU\n') + with self.assertRaisesRegex(RuntimeError, 'status 7: specific failure'): + run_command([sys.executable, '-c', 'import sys; print("specific failure", file=sys.stderr); sys.exit(7)']) + + def test_lookup_siblings_source_tree_path_and_frozen_executable(self): + with tempfile.TemporaryDirectory() as directory: + # macOS temporary directories can already contain symlink components. + root = Path(directory).resolve() + script = root / 'python' / 'script.py' + source = root / 'src' / 'bin' / 'IntaRNA' + sibling = script.parent / 'IntaRNA' + source.parent.mkdir(parents=True) + sibling.parent.mkdir() + linked_directory = root / 'linked-python' + linked_directory.symlink_to(script.parent, target_is_directory=True) + scripts = (script, linked_directory / script.name) + for path in (source, sibling): + path.write_text('#!/bin/sh\nexit 0\n') + path.chmod(0o755) + with patch.dict('os.environ', {'PATH': ''}): + for script_path in scripts: + with self.subTest(script=script_path, lookup='sibling'): + self.assertEqual(find_program('IntaRNA', script_path), str(sibling)) + sibling.unlink() + for script_path in scripts: + with self.subTest(script=script_path, lookup='source'): + self.assertEqual(find_program('IntaRNA', script_path), str(source)) + source.unlink() + for script_path in scripts: + with self.subTest(script=script_path, lookup='missing'), self.assertRaises(RuntimeError): + find_program('IntaRNA', script_path) + sibling.write_text('#!/bin/sh\nexit 0\n') + sibling.chmod(0o755) + for script_path in scripts: + with self.subTest(script=script_path, lookup='frozen'), \ + patch.object(sys, 'frozen', True, create=True), \ + patch.object(sys, 'executable', str(script_path.parent / 'IntaRNApvalue')): + self.assertEqual(find_program('IntaRNA'), str(sibling)) + + def test_optional_builder_uses_paths_and_propagates_failures(self): + spec = importlib.util.spec_from_file_location('companion_builder', ROOT / 'bin' / 'build.py') + builder = importlib.util.module_from_spec(spec) + spec.loader.exec_module(builder) + with patch.object(builder.subprocess, 'run') as run: + self.assertEqual(builder.main(), 0) + command = run.call_args.args[0] + self.assertEqual(command[:3], [sys.executable, '-m', 'PyInstaller']) + self.assertEqual(command[command.index('--paths') + 1], str(ROOT)) + self.assertTrue(run.call_args.kwargs['check']) + temporary = Path(command[command.index('--workpath') + 1]) + self.assertFalse(temporary.exists()) + with patch.object(builder.subprocess, 'run', side_effect=subprocess.CalledProcessError(7, 'PyInstaller')), \ + self.assertRaises(subprocess.CalledProcessError): + builder.main() + + def test_all_standalone_help_without_optional_dependencies(self): + for name in ('IntaRNA-fuse.py', 'CopomuS.py', 'IntaRNA_CSV_p-value.py', 'IntaRNA_plotRegions.py'): + result = subprocess.run([sys.executable, '-S', str(ROOT / name), '--help'], capture_output=True, text=True) + self.assertEqual(result.returncode, 0, result.stderr) + self.assertIn('usage:', result.stdout) + + +if __name__ == '__main__': + unittest.main() diff --git a/python/tests/test_copomus.py b/python/tests/test_copomus.py new file mode 100644 index 00000000..11ffcfc3 --- /dev/null +++ b/python/tests/test_copomus.py @@ -0,0 +1,156 @@ +"""CopomuS regressions for coordinates, mutation selection, profiles and ranking.""" + +import csv +import io +from pathlib import Path +import sys +import tempfile +import unittest +from unittest.mock import patch + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from CopomuS import CopomuS, main +from copomus.IntaRNA import IntaRNA +from copomus.candidate_filters import get_filter, filter_lp_mfe +from copomus.candidate_selectors import get_mutation_enc +from copomus.indexing import array_index_to_idx, idx_to_array_index, add_sub_idx +from copomus.measures import MinDeltaEnergy, MFECover, Measure +from copomus.mutation import Mutation, parse_mutation_encoding +from copomus.mutation_generators import gen_specific + + +class CopomuSTests(unittest.TestCase): + def test_signed_coordinates_round_trip(self): + for offset, expected in ((-3, [-3, -2, -1, 1, 2]), (0, [0, 1, 2, 3, 4]), + (2, [2, 3, 4, 5, 6])): + for i, value in enumerate(expected): + self.assertEqual(array_index_to_idx(i, offset), value) + self.assertEqual(idx_to_array_index(value, offset), i) + self.assertEqual(add_sub_idx(-1, 1), 1) + self.assertEqual(add_sub_idx(1, -1), -1) + with self.assertRaises(ValueError): + idx_to_array_index(0, -1) + + def test_zero_energies_ids_and_quoted_csv(self): + rows = IntaRNA.output_to_dict('id1;id2;start1;E;ED1;other\n"a;b";001;0;0.0;NA;\n') + self.assertEqual(rows, [{'id1': 'a;b', 'id2': '001', 'start1': 0, 'E': 0.0, + 'ED1': None, 'other': None}]) + for empty in ('id1\n', 'E\n', 'E;ED1\n\n'): + self.assertIs(IntaRNA.output_to_dict(empty), False) + for bad in ('E;E\n1;2\n', 'E;ED1\n1\n', 'E\nnan\n', 'E\ninf\n'): + with self.assertRaises(ValueError): + IntaRNA.output_to_dict(bad) + + def test_version_and_no_empty_subprocess_arguments(self): + self.assertEqual(IntaRNA.version_tuple('4.0.0-dev'), (4, 0, 0)) + with patch('copomus.IntaRNA.find_program', return_value='/tmp/with spaces/IntaRNA'), \ + patch('copomus.IntaRNA.run_command', return_value='E\n0\n') as run: + self.assertEqual(IntaRNA().run('GC', 'GC', 1, 1, 'E', 0), {'E': 0}) + self.assertNotIn('', run.call_args.args[0]) + IntaRNA().run('GC', 'GC', 1, 1, 'E', 0, param_file='with spaces.txt') + self.assertIn('--parameterFile=with spaces.txt', run.call_args.args[0]) + + def test_explicit_mutations_and_bounds(self): + args = ('GC', 'GC', 1, 1) + self.assertEqual(get_mutation_enc(*args, '', 2, enc='G1C&C2G'), [(1, 2)]) + mutations = gen_specific(*args, ['E'], 1, 1, 'G1C&C2G') + self.assertEqual((mutations[0].qm, mutations[0].tm), ('CC', 'GG')) + for bad in ('G0C&C2G', 'C-1G&C2G', 'G3C&C2G', 'G1N&C2G', 'G1G&C2G', + 'A1C&C2G', 'G1C&C2G,', 'bad'): + with self.subTest(bad=bad), self.assertRaises(ValueError): + parse_mutation_encoding(bad, *args) + with self.assertRaises(ValueError): + Mutation([], 'GC', 'GC', 1, 1, -1, 2, 'CG') + mut = gen_specific('AGC', 'GCU', -2, -2, [], 1, 1, 'A-2U&U1A')[0] + self.assertEqual((mut.qm, mut.tm), ('UGC', 'GCA')) + + def test_filters_across_zero_and_gc_alias(self): + self.assertIs(get_filter('GC'), get_filter('CG')) + self.assertEqual(get_filter('CG')([(1, 2)], 'GC', 'GC', 1, 1), []) + pairs = [(-1, 1), (1, -1)] + self.assertEqual(filter_lp_mfe(pairs, 'GC', 'GC', -1, -1), pairs) + zero_pairs = [(0, 1), (1, 0)] + self.assertEqual(filter_lp_mfe(zero_pairs, 'GC', 'GC', 0, 0), zero_pairs) + + def test_zero_differences_total_order_and_missing_scores(self): + measures = [] + for values in ((-3, 0, 0, -3), (-3, -3, -3, -3), (-3, -4, -4, -3)): + measure = MinDeltaEnergy() + measure.score.update(zip(('ww', 'wm', 'mw', 'mm'), values)) + measure.post_compute() + measures.append(measure) + self.assertEqual([m.score['min'] for m in measures], [3, 0, -1]) + self.assertEqual(sorted(measures[::-1]), measures) + missing = MinDeltaEnergy() + missing.post_compute() + self.assertTrue(missing.non_i()) + cover = MFECover() + cover.calculate_score(False, 'wm') + self.assertIs(cover.score['wm'], False) + measure = Measure() + measure.score.update(ww=-5, mm=-2, wm=-1.5, mw=-1.5) + self.assertTrue(measure.is_valid(1, 0)) + + def test_profiles_shared_and_wildtype_reused(self): + app = CopomuS() + app.measures = ['EDqi', 'cEDqi', 'Eqi', 'Eti'] + muts = gen_specific('GC', 'GC', 1, 1, app.measures, 1, 1, 'G1C&C2G,C2G&G1C') + paths_seen = [] + def prediction(*args): + for option in args[-1]: + pipe, filename = option.removeprefix('--out=').split(':', 1) + path = Path(filename) + paths_seen.append(path) + path.write_text('# header\n# columns\n1\t0\n2\t1\n' if pipe == 'qAcc' + else 'idx;nt;E\n1;G;-2\n2;C;-3\n') + return {'id1': 'target'} + with patch.object(IntaRNA, 'run', side_effect=prediction) as run: + app.run_tests(muts) + self.assertEqual(run.call_count, 7) # one wildtype + three per mutation + self.assertTrue(all(not path.exists() for path in paths_seen)) + for mut in muts: + self.assertEqual(mut.measures[0].score, mut.measures[1].score | {'is_valid': 0}) + self.assertIsNotNone(mut.measures[0].score['ww']) + + def test_profile_na_and_malformed_rows(self): + measure = Measure() + measure.profiles = {'qAcc': '# header\n# columns\n1\tNA\n', + 'qMinE': 'idx;nt;E\n1;G;0\n'} + self.assertIsNone(measure.profile_value('qAcc', 0)) + self.assertEqual(measure.profile_value('qMinE', 0), 0) + with self.assertRaises(ValueError): + measure.profile_value('qAcc', 1) + measure.profiles['qAcc'] = '# header\n# columns\n1\tinf\n' + with self.assertRaises(ValueError): + measure.profile_value('qAcc', 0) + + def test_cli_validation_and_failed_output_preservation(self): + with tempfile.TemporaryDirectory() as directory: + output = Path(directory) / 'results.tsv' + output.write_text('previous results') + for extra in (['-q', '', '-t', 'GC'], ['-q', 'GC', '-t', 'GC', '--threads=-1'], + ['-q', 'GC', '-t', 'GC', '-d', ';;']): + with patch('sys.stderr', new_callable=io.StringIO), self.assertRaises(SystemExit): + CopomuS().parse_cmd_args([*extra, '-o', str(output)]) + self.assertEqual(output.read_text(), 'previous results') + with patch.object(IntaRNA, 'check_binary'), \ + patch.object(IntaRNA, 'run', side_effect=RuntimeError('prediction failed')), \ + patch('sys.stderr', new_callable=io.StringIO): + self.assertEqual(main(['-q', 'GC', '-t', 'GC', '--mutationEncoding', 'G1C&C2G', + '-o', str(output)]), 1) + self.assertEqual(output.read_text(), 'previous results') + + def test_cli_explicit_mutation_zero_energy_csv(self): + stdout = io.StringIO() + with patch.object(IntaRNA, 'check_binary'), patch.object(IntaRNA, 'run', return_value={'E': 0}), \ + patch('sys.stdout', stdout): + self.assertEqual(main(['-q', 'GC', '-t', 'GC', '--mutationEncoding', 'G1C&C2G', + '-m', 'minDeltaE']), 0) + row = list(csv.DictReader(io.StringIO(stdout.getvalue()), delimiter='\t'))[0] + self.assertEqual(row['mutation'], 'G1C&C2G') + self.assertEqual(row['minDeltaE_min'], '0') + self.assertNotIn('\r', stdout.getvalue()) + + +if __name__ == '__main__': + unittest.main() diff --git a/python/tests/test_intarnapvalue.py b/python/tests/test_intarnapvalue.py new file mode 100644 index 00000000..55fa7fbc --- /dev/null +++ b/python/tests/test_intarnapvalue.py @@ -0,0 +1,195 @@ +"""Sampling and statistical regressions, independent of an installed IntaRNA.""" + +from collections import Counter +import io +import json +import math +import os +from pathlib import Path +import random +import subprocess +import sys +import tempfile +import unittest +from unittest.mock import patch + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT)) +from intarnapvalue.intarna_pvalue import IntaRNApvalue, main +from intarnapvalue.dinucl_shuffle import dinucl_shuffle +from intarna_tools import read_sequence + +QUERY = 'AGGAUGGGGGA' +TARGET = 'CCGCUUCUCGUUCCU' +ARGS = ['-q', QUERY, '-t', TARGET, '-s', '5', '-m', 'b', '--randSeed', 'reproducible'] + + +def fasta_sequences(text): + return text.splitlines()[1::2] + + +class PvalueTests(unittest.TestCase): + def tool(self, *args): + with patch.object(IntaRNApvalue, 'find_binary', return_value='/tmp/IntaRNA with spaces'): + return IntaRNApvalue([*ARGS, *args]) + + def test_normalized_wrapped_single_fasta_and_invalid_inputs(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / 'seq.fa' + path.write_text('\ufeff>sequence\nacgt\ntgca\n') + self.assertEqual(read_sequence(str(path)), 'ACGUUGCA') + self.assertEqual(read_sequence('a' * 10000), 'A' * 10000) + for invalid in ('', '>empty\n', 'ACGU\n', '>one\nAC\n>two\nGU\n', '>one\nACN\n'): + path.write_text(invalid) + with self.assertRaises(ValueError): + read_sequence(str(path)) + for args in (['-q', ''], ['-q', 'AX'], ['-s', '0'], ['--threads=-1'], ['--max-attempts', '4']): + with patch('sys.stderr', new_callable=io.StringIO), self.assertRaises(SystemExit): + self.tool(*args) + + def test_shuffle_composition_endpoints_and_singletons(self): + rng = random.Random(13) + for seq in ('A', 'acgt', 'G' * 1000, QUERY, TARGET): + normalized = seq.upper().replace('T', 'U') + for _ in range(30): + shuffled = dinucl_shuffle(seq, rng) + self.assertEqual(len(shuffled), len(normalized)) + self.assertEqual((shuffled[0], shuffled[-1]), (normalized[0], normalized[-1])) + self.assertEqual(Counter(shuffled), Counter(normalized)) + self.assertEqual(Counter(zip(shuffled, shuffled[1:])), + Counter(zip(normalized, normalized[1:]))) + for seq in ('', 'N', 'AX'): + with self.assertRaises(ValueError): + dinucl_shuffle(seq) + + def test_local_rng_deterministic_and_no_global_reseed(self): + state = random.getstate() + a, b = self.tool(), self.tool() + self.assertEqual(random.getstate(), state) + self.assertEqual(a.shuffle_sequence(QUERY, 20, a.rng), b.shuffle_sequence(QUERY, 20, b.rng)) + self.assertEqual(random.getstate(), state) + + def test_both_strands_independent_pairwise_and_cleanup(self): + tool = self.tool() + seen = [] + def predict(query, target, input_text=None, *, pairwise=False): + self.assertTrue(pairwise) + self.assertIsNone(input_text) + q, t = Path(query), Path(target) + seen.extend((q, t)) + queries, targets = fasta_sequences(q.read_text()), fasta_sequences(t.read_text()) + self.assertEqual(len(queries), 5) + self.assertEqual(len(targets), 5) + self.assertGreater(len(set(queries)), 1) + for sequences, original in ((queries, QUERY), (targets, TARGET)): + for seq in sequences: + self.assertEqual(Counter(seq), Counter(original)) + self.assertEqual(Counter(zip(seq, seq[1:])), Counter(zip(original, original[1:]))) + return [(str(i), str(i), -i) for i in reversed(range(5))] + with patch.object(tool, '_predict', side_effect=predict): + self.assertEqual(tool.get_scores(), ([0, -1, -2, -3, -4], 0)) + self.assertTrue(all(not p.exists() for p in seen)) + + def test_each_single_strand_mode_and_retry_counts(self): + for mode in ('q', 't'): + tool = self.tool('-m', mode) + attempts = [] + def predict(query, target, input_text=None): + self.assertEqual((query, target), ('STDIN', TARGET) if mode == 'q' else (QUERY, 'STDIN')) + seqs = fasta_sequences(input_text) + attempts.append(len(seqs)) + # Reject one pair on the first pass, then fill the remainder. + start = 1 if len(attempts) == 1 else 0 + return [('fixed' if mode == 'q' else str(i), str(i) if mode == 'q' else 'fixed', -i) + for i in range(start, len(seqs))] + with patch.object(tool, '_predict', side_effect=predict): + scores, skipped = tool.get_scores() + self.assertEqual((len(scores), skipped), (5, 1)) + self.assertEqual(attempts, [5, 1]) + + def test_bounded_no_interactions_and_fail_fast(self): + tool = self.tool('--max-attempts', '7') + with patch.object(tool, '_predict', return_value=[]) as run, self.assertRaisesRegex(RuntimeError, '7 attempts'): + tool.get_scores() + self.assertEqual(run.call_count, 2) + with patch.object(tool, '_predict', side_effect=RuntimeError('command failed')) as run, \ + self.assertRaisesRegex(RuntimeError, 'command failed'): + tool.get_scores() + self.assertEqual(run.call_count, 1) + for rows in ([('0', '0', -1), ('0', '0', -2)], [('0', '1', -1)], [('5', '5', -1)]): + with patch.object(tool, '_predict', return_value=rows), self.assertRaises(ValueError): + tool.get_scores() + + def test_parameter_file_all_predictions_and_zero_energy(self): + tool = self.tool() + tool.parameterFile = '/tmp/parameter file with spaces' + with patch('intarnapvalue.intarna_pvalue.run_command', return_value='id1;id2;E\nt;q;0\n') as run: + self.assertEqual(tool.get_single_score(QUERY, TARGET), 0) + tool._predict('STDIN', TARGET, '>0\nAAA\n') + for call in run.call_args_list: + command = call.args[0] + self.assertIn('--parameterFile=/tmp/parameter file with spaces', command) + self.assertIn('--outNumber=1', command) + self.assertNotIn('', command) + for bad in ('E\n-1\n', 'id1;id2;E\na;b;nan\n', 'id1;id2;E\na;b\n'): + with patch('intarnapvalue.intarna_pvalue.run_command', return_value=bad), self.assertRaises(ValueError): + tool.get_single_score(QUERY, TARGET) + with patch('intarnapvalue.intarna_pvalue.run_command', return_value='id1;id2;E\n'): + self.assertIsNone(tool.get_single_score(QUERY, TARGET)) + + def test_cdf_matches_independent_legacy_integration(self): + from scipy import integrate, stats + scores = [-1.235, -1.435645, -6.234234, -12.999, -15.23, -6.98, -6.23, -2.78] + for name, calculate in (('norm', IntaRNApvalue.calculate_pvalue_gauss), + ('gumbel_l', IntaRNApvalue.calculate_pvalue_gumbel), + ('genextreme', IntaRNApvalue.calculate_pvalue_gev)): + distribution = getattr(stats, name) + params = distribution.fit(scores) + expected = integrate.quad(lambda x: distribution.pdf(x, *params), -math.inf, -10)[0] + self.assertAlmostEqual(calculate(-10, scores), expected, places=7) + self.assertLessEqual(calculate(-20, scores), calculate(-10, scores)) + self.assertEqual(IntaRNApvalue.calculate_pvalue_empirical(-10, scores), 0.25) + self.assertEqual(IntaRNApvalue.calculate_pvalue_empirical(0, [0]), 1) + for scores in ([], [float('nan')], [float('inf')], [-1, -1, -1], [-1, -2]): + for method in (IntaRNApvalue.calculate_pvalue_gev, IntaRNApvalue.calculate_pvalue_gauss, + IntaRNApvalue.calculate_pvalue_gumbel): + with self.assertRaises(ValueError): + method(-1, scores) + + def test_standalone_shuffle_cli(self): + script = ROOT / 'intarnapvalue' / 'dinucl_shuffle.py' + for seq, expected in (('a', 'A\nA\n'), ('t', 'U\nU\n')): + result = subprocess.run([sys.executable, '-S', str(script), seq, '--samples', '2'], + cwd=tempfile.gettempdir(), capture_output=True, text=True) + self.assertEqual(result.returncode, 0, result.stderr) + self.assertEqual(result.stdout, expected) + result = subprocess.run([sys.executable, '-S', str(script), 'AX'], capture_output=True, text=True) + self.assertEqual(result.returncode, 2) + self.assertNotIn('Traceback', result.stderr) + + def test_cli_help_without_binary_or_dependencies(self): + result = subprocess.run([sys.executable, '-S', '-m', 'intarnapvalue', '--help'], + env={**os.environ, 'PYTHONPATH': str(ROOT), 'PATH': ''}, + capture_output=True, text=True) + self.assertEqual(result.returncode, 0, result.stderr) + self.assertIn('--samples', result.stdout) + + def test_cli_empirical_scores_and_errors(self): + with patch.object(IntaRNApvalue, 'find_binary', return_value='IntaRNA'), \ + patch.object(IntaRNApvalue, 'get_single_score', return_value=0), \ + patch.object(IntaRNApvalue, 'get_scores', return_value=([-1, 0, -2, -3, -4], 2)), \ + patch('sys.stdout', new_callable=io.StringIO) as stdout, \ + patch('sys.stderr', new_callable=io.StringIO) as stderr: + self.assertEqual(main([*ARGS, '-d', 'empirical']), 0) + self.assertEqual(stdout.getvalue(), '1.0\n') + self.assertIn('conditional', stderr.getvalue()) + with patch.object(IntaRNApvalue, 'find_binary', side_effect=RuntimeError('binary missing')), \ + patch('sys.stdout', new_callable=io.StringIO) as stdout, \ + patch('sys.stderr', new_callable=io.StringIO) as stderr: + self.assertEqual(main(ARGS), 1) + self.assertEqual(stdout.getvalue(), '') + self.assertIn('binary missing', stderr.getvalue()) + + +if __name__ == '__main__': + unittest.main() diff --git a/python/tests/test_native_tools.py b/python/tests/test_native_tools.py new file mode 100644 index 00000000..1c34d9b0 --- /dev/null +++ b/python/tests/test_native_tools.py @@ -0,0 +1,73 @@ +"""Optional integration tests against a native IntaRNA binary on PATH.""" + +import csv +import io +import os +from pathlib import Path +import shutil +import subprocess +import sys +import tempfile +import unittest + +ROOT = Path(__file__).resolve().parents[1] +QUERY = 'GCUGAAAAACAUAACCCAUAAAAUGCUAGCUGUACCAGGAACCA' +TARGET = 'GGUUUCUUCGCCUCUGCGUUCACCAAAGUGUUCACCC' + + +@unittest.skipUnless(shutil.which('IntaRNA'), 'native IntaRNA is not on PATH') +class NativeToolsTests(unittest.TestCase): + def run_tool(self, command, *, success=True): + result = subprocess.run([sys.executable, *command], cwd=ROOT.parent, + env={**os.environ, 'PYTHONPATH': str(ROOT)}, + capture_output=True, text=True, timeout=90) + self.assertEqual(result.returncode == 0, success, result.stderr) + self.assertNotIn('Traceback', result.stderr) + return result + + def test_sampling_modes_reproducible_across_threads(self): + for mode in ('q', 't', 'b'): + args = ['-m', 'intarnapvalue', '-q', QUERY, '-t', TARGET, '-s', '5', + '-m', mode, '--randSeed', 'native', '-o', 'scores'] + a = self.run_tool([*args, '--threads', '1']) + b = self.run_tool([*args, '--threads', '2']) + self.assertEqual(a.stdout, b.stdout) + self.assertEqual(len([float(x) for x in a.stdout.splitlines()]), 5) + + def test_native_no_interaction_attempt_limit(self): + result = self.run_tool(['-m', 'intarnapvalue', '-q', 'AAAAAAAA', '-t', 'AAAAAAAA', + '-s', '1', '-m', 'b', '-o', 'scores', '--threads', '1', + '--max-attempts', '2'], success=False) + self.assertEqual(result.stdout, '') + self.assertIn('2 attempts', result.stderr) + + def test_native_parameter_file_and_empirical(self): + with tempfile.TemporaryDirectory() as directory: + params = Path(directory) / 'params with spaces.txt' + params.write_text('seedBP=3\nqAcc=N\ntAcc=N\n') + result = self.run_tool(['-m', 'intarnapvalue', '-q', QUERY, '-t', TARGET, + '-s', '5', '-m', 'b', '--threads', '1', '-d', 'empirical', + '--randSeed', 'test', '-p', str(params)]) + self.assertTrue(0 <= float(result.stdout) <= 1) + + def test_native_copomus_profiles_signed_indices_and_explicit_mutation(self): + result = self.run_tool([str(ROOT / 'CopomuS.py'), '-q', QUERY, '-t', TARGET, + '--threads', '1', '-m', 'EDqi', '-m', 'cEDqi', '-m', 'Eqi', '-m', 'Eti', + '--qIdxPos0', '-20', '--tIdxPos0', '-15', '-f', 'lpMfe']) + rows = list(csv.DictReader(io.StringIO(result.stdout), delimiter='\t')) + self.assertTrue(rows) + for row in rows: + self.assertEqual(row['EDqi_ww'], row['cEDqi_ww']) + self.assertNotEqual(row['qIndex'], '0') + self.assertNotEqual(row['tIndex'], '0') + encoding = rows[0]['mutation'] + result = self.run_tool([str(ROOT / 'CopomuS.py'), '-q', QUERY, '-t', TARGET, + '--threads', '1', '--qIdxPos0', '-20', '--tIdxPos0', '-15', + '--mutationEncoding', encoding]) + explicit = list(csv.DictReader(io.StringIO(result.stdout), delimiter='\t')) + self.assertEqual(len(explicit), 1) + self.assertEqual(explicit[0]['mutation'], encoding) + + +if __name__ == '__main__': + unittest.main() diff --git a/python/tests/test_output_tools.py b/python/tests/test_output_tools.py new file mode 100644 index 00000000..d4275962 --- /dev/null +++ b/python/tests/test_output_tools.py @@ -0,0 +1,258 @@ +"""Regression tests for the Python replacements of the R output tools.""" + +import csv +import importlib.util +import math +import os +from pathlib import Path +import subprocess +import sys +import tempfile +import unittest +from unittest.mock import patch + +import numpy as np + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT)) + + +def load_script(name, filename): + spec = importlib.util.spec_from_file_location(name, ROOT / filename) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +pvalue = load_script("csv_pvalue", "IntaRNA_CSV_p-value.py") +plot = load_script("plot_regions", "IntaRNA_plotRegions.py") + + +class OutputToolsTests(unittest.TestCase): + def setUp(self): + self.directory = tempfile.TemporaryDirectory() + self.addCleanup(self.directory.cleanup) + self.path = Path(self.directory.name) + + def write(self, name, text): + path = self.path / name + path.write_text(text, encoding="utf-8") + return path + + def run_cli(self, script, *arguments, success=True): + result = subprocess.run([sys.executable, str(ROOT / script), *map(str, arguments)], + capture_output=True, text=True, + env={**os.environ, "MPLCONFIGDIR": str(self.path / "mpl"), + "DISPLAY": ""}) + self.assertEqual(result.returncode == 0, success, result.stdout + result.stderr) + self.assertNotIn("Traceback", result.stderr) + return result + + def energy_table(self, column="E"): + # Deterministic GEV quantiles: independent inverse-CDF formula, xi=0.2. + probabilities = (np.arange(1000) + 0.5) / 1000 + energies = -(10 + 2 * ((-np.log(probabilities)) ** -0.2 - 1) / 0.2) + source = self.path / "energies.csv" + with source.open("w", newline="", encoding="utf-8") as handle: + writer = csv.writer(handle, delimiter=";") + writer.writerow(["id", column, "annotation"]) + for index, energy in enumerate(energies): + writer.writerow(["{:04}".format(index), energy, 'RNA; "α"']) + return source, energies + + def test_gev_tail_sign_gumbel_and_endpoints(self): + energies = np.array([-3.0, -8.0, -12.0, -40.0]) + for xi in (-0.2, 0.2): + # Original R/evir formula, using explicit xi rather than SciPy's c. + y = 1 + xi * (-energies[:3] - 10) / 2 + expected = -np.expm1(-y ** (-1 / xi)) + np.testing.assert_allclose(pvalue.gev_pvalues(energies[:3], (-xi, 10, 2)), expected) + expected = -np.expm1(-np.exp((energies + 10) / 2)) + np.testing.assert_allclose(pvalue.gev_pvalues(energies, (0, 10, 2)), expected) + np.testing.assert_equal(pvalue.gev_pvalues([-100, 100], (0.2, 10, 2)), [0, 1]) + # Tiny upper tails must not be lost to 1-cdf cancellation. + self.assertGreater(pvalue.gev_pvalues([-60], (0, 0, 1))[0], 0) + + def test_fit_recovers_known_distribution(self): + _, energies = self.energy_table() + params = pvalue.fit_gev(energies) + np.testing.assert_allclose(params, [-0.2, 10, 2], atol=0.015, rtol=0) + expected = 1 - (np.arange(1000) + 0.5) / 1000 + np.testing.assert_allclose(pvalue.gev_pvalues(energies, params), expected, atol=0.002) + + def test_bh_unsorted_ties_and_bounds(self): + np.testing.assert_allclose(pvalue.adjust_bh([0.04, 0.01, 0.03, 0.01, 1, 0]), + [0.048, 0.02, 0.045, 0.02, 1, 0]) + + def test_legacy_r_reference(self): + # R 4.4.3, original R/IntaRNA_CSV_p-value.R at 2c4d0ec, 1000 inverse-CDF + # quantiles per shape, location=10, scale=2. Compare rows 0/99/499/899/999. + # Optimizers differ; BH can amplify the seven-decimal tail differences. + references = { + -0.2: ([0.9994943, 0.9006078, 0.5007587, 0.1001753, 0.0004651], + [0.9994943, 0.9994715, 0.9994715, 0.9918347, 0.4651]), + 0.0: ([0.9995020, 0.9006536, 0.5005593, 0.1003176, 0.0004918], + [0.9995020, 0.9995020, 0.9991204, 0.9932436, 0.4918]), + 0.2: ([0.9995087, 0.9006690, 0.5003434, 0.1003388, 0.0004986], + [0.9995087, 0.9995087, 0.9986894, 0.9934535, 0.4986]), + } + quantiles = (np.arange(1000) + 0.5) / 1000 + for xi, (expected_p, expected_fdr) in references.items(): + with self.subTest(xi=xi): + energies = (-(10 + 2 * ((-np.log(quantiles)) ** -xi - 1) / xi) + if xi else -(10 - 2 * np.log(-np.log(quantiles)))) + actual = np.round(pvalue.gev_pvalues(energies, pvalue.fit_gev(energies)), 7) + fdr = np.round(pvalue.adjust_bh(actual), 7) + indices = [0, 99, 499, 899, 999] + np.testing.assert_allclose(actual[indices], expected_p, atol=1e-4, rtol=0) + np.testing.assert_allclose(fdr[indices], expected_fdr, atol=0.002, rtol=0) + + def test_csv_default_in_place_preserves_fields(self): + source, energies = self.energy_table() + source.chmod(0o640) + self.run_cli("IntaRNA_CSV_p-value.py", source) + with source.open(newline="", encoding="utf-8") as handle: + reader = csv.DictReader(handle, delimiter=";") + rows = list(reader) + self.assertEqual(reader.fieldnames, ["id", "E", "annotation", "p-value", "fdr"]) + self.assertEqual(len(rows), 1000) + self.assertEqual(rows[0]["id"], "0000") + self.assertEqual(rows[0]["annotation"], 'RNA; "α"') + self.assertEqual(float(rows[0]["E"]), energies[0]) + self.assertLess(float(rows[-1]["p-value"]), float(rows[0]["p-value"])) + self.assertTrue(all(0 <= float(row["p-value"]) <= float(row["fdr"]) <= 1 for row in rows)) + self.assertEqual(source.stat().st_mode & 0o777, 0o640) + self.assertNotIn(b"\r", source.read_bytes()) + + def test_csv_separate_output_custom_energy_and_rounding(self): + source, _ = self.energy_table("normalized energy") + before = source.read_bytes() + output = self.path / "extended.csv" + self.run_cli("IntaRNA_CSV_p-value.py", source, output, "normalized energy") + self.assertEqual(source.read_bytes(), before) + with output.open(newline="") as handle: + rows = list(csv.DictReader(handle, delimiter=";")) + expected = pvalue.adjust_bh([float(row["p-value"]) for row in rows]) + self.assertEqual([row["fdr"] for row in rows], [format(value, ".7f") for value in expected]) + + def test_invalid_csv_does_not_replace_input_or_output(self): + invalid = ["", "id;score\na;-5\n", "E;E\n-5;-6\n", "E\n", + "E\n-5\n", "E\n-5\n-5\n-5\n", "E\n-5\n-6\nNA\n", + "E\n-5\n-6\nnan\n", "E\n-5\n-6\ninf\n", "E\n-5;extra\n", + "E;p-value\n-5;0.1\n"] + for index, content in enumerate(invalid): + with self.subTest(content=content): + source = self.write("invalid{}.csv".format(index), content) + output = self.write("output.csv", "keep me\n") + for destination in (source, output): + before = destination.read_bytes() + with self.assertRaises(ValueError): + pvalue.extend_csv(source, destination) + self.assertEqual(destination.read_bytes(), before) + + def test_fit_failure_does_not_overwrite(self): + source, _ = self.energy_table() + before = source.read_bytes() + with patch.object(pvalue, "fit_gev", side_effect=ValueError("fit failed")): + with self.assertRaisesRegex(ValueError, "fit failed"): + pvalue.extend_csv(source, source) + self.assertEqual(source.read_bytes(), before) + + def test_publish_failure_cleans_temporary_output(self): + source, _ = self.energy_table() + before = source.read_bytes() + with patch("intarna_tools.os.replace", side_effect=OSError("cannot replace")): + with self.assertRaises(OSError): + pvalue.extend_csv(source, source) + self.assertEqual(source.read_bytes(), before) + self.assertEqual(list(self.path.iterdir()), [source]) + + def test_cli_errors_and_help(self): + for script in ("IntaRNA_CSV_p-value.py", "IntaRNA_plotRegions.py"): + self.assertIn("usage:", self.run_cli(script, "--help").stdout) + self.run_cli(script, success=False) + bad = self.write("bad.csv", "E\nNaN\n") + self.assertIn("nonfinite", self.run_cli("IntaRNA_CSV_p-value.py", bad, success=False).stderr) + self.run_cli("IntaRNA_plotRegions.py", bad, "1", self.path / "out.png", success=False) + self.run_cli("IntaRNA_plotRegions.py", bad, "1", self.path / "out.xyz", success=False) + + def test_region_selectors_and_default_row_limit(self): + source = ROOT / "examples/regions.csv" + self.assertEqual(plot.read_regions(source, plot.read_parameters("1"))[0], ("b0001", 266, 273)) + self.assertEqual(plot.read_regions(source, plot.read_parameters("2"))[0], ("query", 116, 123)) + source = self.write("many.csv", "id1;start1;end1\n" + "001;2;5\n" * 205) + regions = plot.read_regions(source, plot.read_parameters("1")) + self.assertEqual(len(regions), 200) + self.assertEqual(regions[0][0], "001") + + def test_parameters_custom_columns_comments_defaults_and_limits(self): + path = self.write("plot.params", "# comment\nid = name\nstart = left\nend = right\n" + "title = custom title\nxmin = -5\nxmax = 8\nxvline = 0\nrows = 2\n") + parameters = plot.read_parameters(str(path)) + source = self.write("custom.csv", "name;left;right\nx;4;-2\nx;0;0\ny;10;15\n") + regions = plot.read_regions(source, parameters) + self.assertEqual(regions, [("x", 4, -2), ("x", 0, 0)]) + figure = plot.create_figure(regions, parameters) + import matplotlib.pyplot as plt + self.addCleanup(plt.close, figure) + coverage, intervals = figure.axes + self.assertEqual(coverage.get_title(), "custom title") + self.assertEqual(intervals.get_xlim(), (-5, 8)) + self.assertEqual([tick.get_text() for tick in intervals.get_yticklabels()], ["x", "x"]) + self.assertEqual(list(intervals.lines[0].get_xdata()), [4, -2]) + self.assertEqual(list(intervals.lines[-1].get_xdata()), [0, 0]) + minimal = self.write("minimal.params", "id=name\nstart=left\nend=right\n") + self.assertEqual(plot.read_parameters(str(minimal))["rows"], 200) + + def test_invalid_parameters_and_regions(self): + for content in ("id=x\n", "id=x\nstart=s\nend=e\nrows=0\n", + "id=x\nstart=s\nend=e\nrows=2.5\n", "id=x\nid=y\n", + "id=x\nstart=s\nend=e\nxmin=nan\n", "id=x\nunknown=3\n"): + with self.subTest(content=content): + path = self.write("invalid.params", content) + with self.assertRaises(ValueError): + plot.read_parameters(str(path)) + for content in ("", "id1;start1;end1\n", "id1;start1\na;3\n", + "id1;start1;end1\na;3;NA\n", "id1;start1;end1\na;3.5;4\n", + "id1;start1;end1\na;3;4;extra\n"): + with self.subTest(content=content): + with self.assertRaises(ValueError): + plot.read_regions(self.write("invalid.csv", content), plot.read_parameters("1")) + + def test_density_inclusive_reversed_clipped_and_count_scaled(self): + regions = [("a", -1, 1), ("b", 1, 0)] + x, actual = plot.coverage_density(regions, -1, 1) + # Five bases including both endpoints, with repeated covered positions. + bases = np.array([-1, 0, 1, 1, 0]) + bandwidth = 0.9 * min(np.std(bases, ddof=1), 1 / 1.34) * len(bases) ** -0.2 + expected = sum(np.exp(-0.5 * ((x - base) / bandwidth) ** 2) + / (math.sqrt(2 * math.pi) * bandwidth) for base in bases) + np.testing.assert_allclose(actual, expected) + # X limits remove out-of-range samples before density estimation. + np.testing.assert_allclose(plot.coverage_density(regions, 0, 1)[1], + plot.coverage_density([("a", 0, 1), ("b", 1, 0)], 0, 1)[1]) + self.assertTrue(np.all(plot.coverage_density(regions, 10, 20)[1] == 0)) + + def test_single_base_region(self): + parameters = plot.read_parameters("1") + figure = plot.create_figure([("one", 0, 0)], parameters) + import matplotlib.pyplot as plt + self.addCleanup(plt.close, figure) + self.assertEqual(figure.axes[0].get_xlim(), (-0.5, 0.5)) + self.assertTrue(np.all(np.isfinite(figure.axes[0].lines[0].get_ydata()))) + + def test_every_output_format_without_display(self): + signatures = {".png": b"\x89PNG", ".pdf": b"%PDF", ".svg": b"