Examples, scripts, usage

Two ways of using this package

The dipwmsearch Python package provides

a programming interface that enables one to use its functions inside Python scripts,
a script for searching occurrences of a di-PWM in a sequence.

This gives you two ways of using dipwmsearch. The use inside your own scripts (corresponding to the first way) is documented in the sections below.

The direct use as a script (second way) is exemplified just here. On the command line, just launch python followed by the module name without extension, and all requested parameters: at least a di-PWM file (here ANDR_HUMAN.H11DI.0.A.dpwm) and a FASTA file containing a sequence (here chr15.fa). The last parameter (the real value 0.935), the ratio threshold, is optional.

python -m dipwmsearch ANDR_HUMAN.H11DI.0.A.dpwm chr15.fa 0.935

Output

45518148     GGTACATGCTGTTCTCTTCTCT          28.113017946130746      +
81388496     AGAACTGTCTGTTCCTCTTTGT          28.160033168868388      +
85414568     GGTACAGAGTGTTCCTTTATGT          29.51111375013996       +
95429642     AGCACAGACTGTTCTTTTCTGA          28.03013129722457       +
53515706     AGTACTTTATGTTCTTTTCTGT          28.550965810527835      -
64603163     AGAACATTCTGTACTCTTTTGT          27.488908042195362      -

In all output lines, there are four fields separated by tabulation:

starting position of the occurrence
word or occurrence found at that position
score of the word
strand (+ for Watson or - for Crick strand)

The output can be redirected to a CSV file where the separator is a ‘’ symbol. This file can easily be imported into a spreadsheet, into many database system, or even be read or analysed using Pandas library, using awk, etc.

A starting example

We detail the minimal python lines to use the package, load a di-PWM from a file, then search for and print out all di-PWM occurrences.

Detail of each step

Import the package:
```
import dipwmsearch as ds
```
To parse a diPWM file and create a diPWM object in memory:
```
diP = ds.create_diPwm(diPwm_path_file)
```

In the variable text put the sequence in which you want to search for occurrences

text = "AAATGACATATGTAAATACCTACAACATAAAATTACCCCAGTAAATTTTCTATTTATAAAATTACTACAA"

Use a search function in a loop (it performs the enumeration of valid words and the search with Aho-Corasick algorithm). You must provide the diPWM (variable diP), a sequence to search in as a Python string (variable text), and a threshold ratio, which is a value between 0 and 1. Here we choose 0.95 (in any case a value above 0.8).
```
for start_position, word, score in ds.search_aho_ratio(diP, text, 0.95):
    print(f'{start_position}\t{word}\t{score}')
```
The function is used in a for loop because it returns one occurrence at a time (as an iterator). For each occurrence, it outputs the starting position (variable start_position), the matching word (variable word), and its score (variable score).

Choosing a threshold ratio

To search for occurrences of a di-PWM, one should provide a score threshold. The score is dependent on the content of the di-PWM score matrix. An expert user of a given matrix is able to determine an appropriate score threshold. It means that one cannot use the same score threshold with two different motifs to get “equivalent” results for both. Users prefer to specify a score threshold indirect using a parameter that automatically set score threshold adapated to the motif. There are two alternative ways to do so, using

either a ratio threshold
or a P-value threshold.

Within your Python script, you can use either of them.

Assume we are searching for a motif of length .

The ratio threshold is a real value between 0 and 1 that represents a percentage of the best score. In practice, experts that produce the motifs advise to use a threshold ratio above 0.8 for classical PWM and di-PWM motifs.

An internal function diPwm.set_threshold_from_ratio(ratio) computes the score threshold given a ratio threshold.

Assume you observe an occurrence of a motif with a score . Its P-value is the probability that any word of nucleotides has a score larger than or equal to . The other way round, for a desired P-value of one computes approximately the score such that 99% of the words of length achieves a score smaller than , and 1% achieve a score higher than or equal to .

An internal function diPwm.set_threshold_from_pvalue( pvalue ) computes the score threshold given a P-value threshold.

Attention: choosing a ratio threshold, a P-value threshold, or a score threshold that is too low, can trigger the enumeration of too many valid words. The program may run out of memory. This is especially true if you use less optimized search functions like ds.search_aho_ratio() instead of the block optimized search strategy implemented in the function ds.search_block_optimized_ratio(). In general, the minimum recommended threshold ratio should be above 0.8, although for more selective di-PWM one can use much lower thresholds (for instance a ratio threshold of 0.5 is possible with ATF4_HUMAN.H11DI.0.A.dpwm).

Searching for a di-PWM in a sequence stored in FASTA file

This script takes three parameters on the command line:

the di-PWM file
a FASTA formatted file containing the sequence in which you want to find occurrences
the ratio threshold as an integer value usually 70 and 100.

It uses the SeqIO and Seq modules from BioPython to load the sequence from the FASTA file and store it into a SeqRecord object. To use this script the BioPython must be installed beforehand. Take a look at the documentation on `SeqIO <https://biopython.org/wiki/SeqIO>`_ and `Seq <https://biopython.org/wiki/Seq>`_ if you want more information.

Script

#!/usr/bin/env python3
#_*_ coding:utf-8 _*_

"""
Search for and report all occurrences of a di-PWM motif in a sequence
using package dipwmsearch
and its optimized function ds.search_block_optimized_ratio
Search only on one strand.

Requires as input three parameters on the command line:
1. a di-PWM file
2. a FASTA file containing a sequence
3. a ratio (ie. a value between 0.8 and 1)
"""

# example search with block enumeration strategy and Aho-Corasick automaton
# from simple fasta file
import sys, os
import dipwmsearch as ds

from Bio import SeqIO
from Bio.Seq import Seq

# check number of arguments
if (len(sys.argv) <= 3):
    print( sys.argv[0], " not enough arguments." )
    print( __doc__ )
    sys.exit(1)

# create diPWM object from path
diP = ds.create_diPwm(sys.argv[1])

# create object SeqIO from fasta file
file = open(sys.argv[2])
seqRecord = SeqIO.read(file, "fasta")

# convert sequence text in uppercase
mySeq = str(seqRecord.seq.upper())

# print for each solution : starting position in the sequence, word, score
for i, word, score in ds.search_block_optimized_ratio(diP, mySeq, float(sys.argv[3])):
    print(f'{i}\t{word}\t{score}')

sys.exit(0)

Usage

The command line needed to execute this script in a terminal is below (provided all input files are in the current directory), followed by its output.

python ex_block_opt_fasta_single.py ANDR_HUMAN.H11DI.0.A.dpwm chr15.fa 0.94

Output:

45518148    GGTACATGCTGTTCTCTTCTCT  28.113017946130746
81388496    AGAACTGTCTGTTCCTCTTTGT  28.160033168868388
85414568    GGTACAGAGTGTTCCTTTATGT  29.51111375013996
95429642    AGCACAGACTGTTCTTTTCTGA  28.03013129722457

Searching for a di-PWM using a P-value instead of a ratio threshold

In the previous script, instead of calling the function ds.search_block_optimized_ratio and providing a ratio threshold as third parameter, you call the function ds.search_block_optimized_pvalue and provides a P-value threshold as third parameter. Otherwise, all other lines of the script remain the same.

Searching for a di-PWM on both strand of a sequence

This script is very similar to the previous one, but it does it on both strands (the Watson and Crick strands) and report the overall running time for the search. For each occurrence reported in the output, its line ends with a symbol + if the occurrence was found on the Watson strand, or with a symbol - if the occurrence was found on the other strand.

Script

#!/usr/bin/env python3
#_*_ coding:utf-8 _*_

# docstring
"""
Search for and report all occurrences of a di-PWM motif in a sequence
using package dipwmsearch
and its optimized function ds.search_block_optimized_ratio

Print the motif name
Search occurrences on both strands (report strand information with symbols + or -).
Report the running time.

Requires as input three parameters on the command line:
1. a di-PWM file
2. a FASTA file containing a sequence
3. a ratio (ie. a value between 0.8 and 1)
"""

import os, sys, time
from Bio import SeqIO
from Bio.Seq import Seq
from dipwmsearch.diPwm import diPWM, create_diPwm
from dipwmsearch.Block import search_block_optimized

################################## function search_block_optimized_ratio_time
def search_block_optimized_ratio_time(pathDiPwm, pathFasta, ratio):
    start_time = time.time()

    # Creation objet diPwm
    diP = create_diPwm(pathDiPwm)
    nameDiPwm = os.path.basename(pathDiPwm)
    print(f'Name diPWM : {nameDiPwm}')

    # Threshold défini en fonction du ratio
    threshold = diP.set_threshold_from_ratio(ratio)

    # Ouverture du fichier Fasta et conversion en objet seqIO
    file = open(pathFasta)
    seqRecord = SeqIO.read(file, "fasta")
    seq_name = os.path.basename(pathFasta)[:-3]

    # Convert sequence text in uppercase
    text = str(seqRecord.seq.upper())

    # Create seq Reverse Complement et convert in text
    seqRC = seqRecord.seq.reverse_complement()
    textRC = str(seqRC.upper())

    # Recherche dnas le brin +
    for position, word, score in search_block_optimized(diP, text, threshold):
            print(f'{position} \t {word} \t {score} \t +')
    # Recherche dans le brin -
    for position, word, score in search_block_optimized(diP, textRC, threshold):
            print(f'{position} \t {word} \t {score} \t -')

    # Calcul du temps :
    time_search = time.time() - start_time
    print(f'overall running time: {time_search:.2f} s.')


################################## MAIN
# check number of arguments
if (len(sys.argv) <= 3):
    print( sys.argv[0], " not enough arguments." )
    print( __doc__ )
    sys.exit(1)

pathDiPwm = (sys.argv[1])

pathFasta = sys.argv[2]

ratio = float(sys.argv[3])

search_block_optimized_ratio_time(pathDiPwm, pathFasta, ratio)

Usage

The command line needed to execute this script in a terminal is below (provided all input files are in the current directory), followed by its output.

python ex_search_block_opt_both_strands.py ANDR_HUMAN.H11DI.0.A.dpwm chr15.fa 0.935

Output

Name diPWM : ANDR\_HUMAN.H11DI.0.A.dpwm
45518148     GGTACATGCTGTTCTCTTCTCT          28.113017946130746      +
81388496     AGAACTGTCTGTTCCTCTTTGT          28.160033168868388      +
85414568     GGTACAGAGTGTTCCTTTATGT          29.51111375013996       +
95429642     AGCACAGACTGTTCTTTTCTGA          28.03013129722457       +
53515706     AGTACTTTATGTTCTTTTCTGT          28.550965810527835      -
64603163     AGAACATTCTGTACTCTTTTGT          27.488908042195362      -
overall running time: 3.92 s.

Computing P-values of occurrences

The search functions of the package outputs for each occurrence the starting position in the sequence, the valid word, its score (and the strand). It does not output the P-value. Once your search results is saved in a file (say the “results” file), you have the possibility, later on, to add the P-values to your output using additional scripts located in the util subdirectory. For this you will

compute the P-values of valid words above your ratio threshold and stored them in a file (say a P-value file)
then merge your search output file with the results file using the found word as a key.

An example shows the process. Assume you have organized your files as follows in your working directory:

in the pwm subdirectory, use the di-PWM file: ATF4_HUMAN.H11DI.0.A.dpwm
in the seq subdirectory, use the chr15.fa
in the util subdirectory, are located the python, Awk scripts you need.

You have already search for occurrences of this di-PWM in the sequence with the following command

python -m dipwmsearch pwm/ATF4_HUMAN.H11DI.0.A.dpwm seq/chr15.fa 0.9 > chr15_ATF4_90_res.tsv

and your output file in current directory is thus: chr15_ATF4_90_res.tsv

Then you execute the two following command lines:

util/enumerate_valid_words.py pwm/ATF4_HUMAN.H11DI.0.A.dpwm 0.90 | LC_ALL=C sort -g -r -k2,2 | ./util/comp_pvalue.awk  > ATF4_HUMAN_VWS_pvalue_90.tsv

util/merge_pvalue_occ.awk ATF4_HUMAN_VWS_pvalue_90.tsv chr15_ATF4_90_res.tsv > chr15_ATF4_90_res_pval.tsv

Explanations The first line

enumerates valid words for the same ratio threshold used in the search (here 0.90)
sorts the valid words by decreasing score (using the bash sort command)
executes an AWK script (comp_pvalue.awk) to add a column with the P-value for each valid word
then redirects the output in a TSV formatted text file.

The second line executes the script merge_pvalue_occ.awk which

loads the table of valid words and their corresponding P-value in memory
performs a merge between your search results (in file chr15_ATF4_90_res.tsv) with the table, to add a new column containing the P-value of each found occurrence
redirects the output in the file chr15_ATF4_90_res_pval.tsv

This script also prints on the standard error a summary of its work including the number of words for which it found no P-value. If this number is zero, then everything went fine.

The file chr15_ATF4_90_res_pval.tsv differs from the original search output (file chr15_ATF4_90_res.tsv) only by its last column which gives the P-value of the occurrence on that line. Here are the four first lines of this output file; the columns in order contain: position, found word, score, strand (+/-), and P-value.

19874075    AAAATGATGCAAGT  22.512495133029873      +       0.000002484769
19874599    AACATGATGAAACA  20.23612924522815       +       0.000006336719
19902355    GAAATGATGAAATG  20.14393194381094       +       0.000006578863
19902363    GAAATGATGAAATG  20.14393194381094       +       0.000006578863

Requirements:

gawk must be installed
all scripts in util subdirectory must be executable

Remarks: The P-value computation assumes all nucleotides are equiprobable. It does not yet consider the model in which the nucleotides are distributed as in the target genome.

Simplicity

Everything should be made as simple as possible,
but not any simpler ---Albert Einstein