Source code for tfmp

#!/usr/bin/env python
"""
pytfmpval contains Python wrappers for the TFM-Pvalue program. It allows for
accurate calculation of p-values associated with scores for a given motif PWM.
Naturally, it can also be used to calculate the score corresponding to a given
p-value for a transcription factor motif PWM.

The original program can be downloaded at:
http://bioinfo.lifl.fr/tfm-pvalue/tfm-pvalue.php

Additionally, see this paper for how the p-values and thresholds are calculated:
Efficient and accurate P-value computation for Position Weight Matrices
H. Touzet and J.S. Varre
Algorithms for Molecular Biology 2007, 2:15

Copyright (C) 2017  Jared Andrews

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

from __future__ import absolute_import, division, print_function, unicode_literals
import pytfmpval.pytfmpval as tfm
import psutil
from math import ceil


[docs]def create_matrix(matrix_file, bg=[0.25, 0.25, 0.25, 0.25], mat_type="counts", log_type="nat"): """ From a JASPAR formatted motif matrix count file, create a Matrix object. This function also converts it to a log-odds (position weight) matrix if necessary. Args: matrix_file (str): White-space delimited string of row-concatenated motif matrix. bg (list of floats): Background nucleotide frequencies for [A, C, G, T]. mat_type (str): Type of motif matrix provided. Options are: "counts", "pfm", "pwm". "counts" is for raw count matrices for each base at each position. "pfm" is for position frequency matrices (frequencies already calculated. "pwm" is for position weight matrices (also referred to as position-specific scoring matrices.) log_type (str): Base to use for log. Default is to use the natural log. "log2" is the other option. This will affect the scores and p-values. Returns: m (pytfmpval Matrix): Matrix in pwm format. """ a, c, g, t = bg[0], bg[1], bg[2], bg[3] m = tfm.Matrix(a, c, g, t) m.readJasparMatrix(matrix_file) if mat_type.upper() == "COUNTS": if log_type.upper() == "NAT": m.toLogOddRatio() elif log_type.upper() == "LOG2": m.toLog2OddRatio() else: print("Improper log type argument, using natural log.") m.toLogOddRatio() return m
[docs]def read_matrix(matrix, bg=[0.25, 0.25, 0.25, 0.25], mat_type="counts", log_type="nat"): """ From a string of space-delimited counts create a Matrix object. Break the string into 4 rows corresponding to A, C, G, and T. This function also converts it to a log-odds (position weight) matrix if necessary. Args: matrix_file (str): White-space delimited string of row-concatenated motif matrix. bg (list of floats): Background nucleotide frequencies for [A, C, G, T]. mat_type (str): Type of motif matrix provided. Options are: "counts", "pfm", "pwm". "counts" is for raw count matrices for each base at each position. "pfm" is for position frequency matrices (frequencies already calculated). "pwm" is for position weight matrices (also referred to as position-specific scoring matrices.) log_type (str): Base to use for log. Default is to use the natural log. "log2" is the other option. This will affect the scores and p-values. Returns: m (pytfmpval Matrix): Matrix in pwm format. """ try: a, c, g, t = bg[0], bg[1], bg[2], bg[3] if (len(matrix.split()) % 4) != 0: raise ValueError("Uneven rows in motif matrix. Ensure rows of equal length in input.") m = tfm.Matrix(a, c, g, t) m.readMatrix(matrix) if mat_type.upper() == "COUNTS": if log_type.upper() == "NAT": m.toLogOddRatio() elif log_type.upper() == "LOG2": m.toLog2OddRatio() else: print("Improper log type argument, using natural log.") m.toLogOddRatio() return m except ValueError as error: print(repr(error))
[docs]def score2pval(matrix, req_score, mem_thresh=2.0): """ Determine the p-value for a given score for a specific motif PWM. Args: matrix (pytfmpval Matrix): Matrix in pwm format. req_score (float): Requested score for which to determine the p-value. mem_thresh (float): Memory in GBs to remain free to system. Once passed, the closest p-val approximation will be returned instead of the exact p-val. Should only occur rarely with very long and degenerate motifs. Used to help ensure the system won't run out of memory due to these outliers. This is only calculated after each pass, each of which is more time and memory intensive than the last, so changing this value isn't recommended unless accuracy out to the 8th decimal place is really necessary. Returns: pv (float): The calculated p-value corresponding to the score. """ mem_thresh = mem_thresh * 1000 * 1024 * 1024 granularity = 0.1 max_granularity = 1e-10 decrgr = 10 # Factor to increase granularity by after each iteration. pv = tfm.doublep() ppv = tfm.doublep() while granularity > max_granularity: matrix.computesIntegerMatrix(granularity) max_s = int(req_score * matrix.granularity + matrix.offset + matrix.errorMax + 1) min_s = int(req_score * matrix.granularity + matrix.offset - matrix.errorMax - 1) score = int(req_score * matrix.granularity + matrix.offset) matrix.lookForPvalue(score, min_s, max_s, ppv, pv) mem = psutil.virtual_memory() if mem.available <= mem_thresh: print("Memory usage threshold passed, returning closest approximation.") return pv.value() if ppv.value() == pv.value(): return pv.value() granularity = granularity / decrgr print("Max granularity exceeded. Returning closest approximation.") return pv.value()
[docs]def pval2score(matrix, pval, mem_thresh=2.0): """ Determine the score for a given p-value for a specific motif PWM. Args: matrix (pytfmpval Matrix): Matrix in pwm format. pval (float): p-value for which to determine the score. mem_thresh (float): Memory in GBs to remain free to system. Once passed, the closest p-val approximation will be returned instead of the exact p-val. Should only occur rarely with very long and degenerate motifs. Used to help ensure the system won't run out of memory due to these outliers. This is only calculated after each pass, each of which is more time and memory intensive than the last, so changing this value isn't recommended unless accuracy out to the 8th decimal place is really necessary. Returns: score (float): The calculated score corresponding to the p-value. """ mem_thresh = mem_thresh * 1000 * 1024 * 1024 init_granularity = 0.1 max_granularity = 1e-10 decrgr = 10 # Factor to increase granularity by after each iteration. pv = tfm.doublep() # Initialize as a c++ double. ppv = tfm.doublep() matrix.computesIntegerMatrix(init_granularity) max_s = int(matrix.maxScore + ceil(matrix.errorMax + 0.5)) min_s = int(matrix.minScore) granularity = init_granularity while granularity > max_granularity: matrix.computesIntegerMatrix(granularity) score = matrix.lookForScore(min_s, max_s, pval, pv, ppv) min_s = int((score - ceil(matrix.errorMax + 0.5)) * decrgr) max_s = int((score + ceil(matrix.errorMax + 0.5)) * decrgr) if ppv.value() == pv.value(): break mem = psutil.virtual_memory() if mem.available <= mem_thresh: print("Memory usage threshold passed, returning closest score approximation.") break granularity = granularity / decrgr if granularity <= max_granularity: print("Max granularity exceeded. Returning closest score approximation.") final_score = (score - matrix.offset) / matrix.granularity return final_score