<p>I want to map TFBS (transcription factor binding sites) on to the sequences from dataset containing 8448 (210 base) sequences, but it is  not giving all the hits I need</p>
<pre><code>import pickle
from pickle import load, dump
from collections import defaultdict
from Bio.Seq import Seq

#first element is BREu = SSRCGCC, length = 7 bases, size of matrix = 4*7, max score = 3.33
# first row = A, second row = T, third row = G and forth row = C
BREu = [
    [0.0, 0.0, 0.5, 0, 0, 0, 0],
    [0.0, 0.0, 0.0, 0, 0, 0, 0],
    [0.5, 0.5, 0.5, 0, 1, 0, 0],
    [0.5, 0.5, 0.0, 1, 0, 1, 1]]

#second element is ETS = MGGAART, length = 7 bases, size of matrix = 4*7, max score = 6
ETS = [
    [0.5, 0, 0, 1, 1, 0.5, 0],
    [0.0, 0, 0, 0, 0, 0.0, 1],
    [0.0, 1, 1, 0, 0, 0.5, 0],
    [0.5, 0, 0, 0, 0, 0.0, 0]]

#Third element is INR = YYANWYY, length = 7 bases, size of matrix = 4*7, max score = 3.75
INR = [
    [0.0, 0.0, 1, 0.25, 0.5, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.5, 0.5, 0.5],
    [0.0, 0.0, 0, 0.25, 0.0, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.0, 0.5, 0.5]]

#forth element is NFY = CCAAT, length = 5 bases, size of matrix = 4*5, max score = 5
NFY = [
    [0, 0, 1, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0],
    [1, 1, 0, 0, 0]]

#seventh element is SP1_NC = GGGNGGG, length = 7 bases, size of matrix = 4*7, max score = 7
#Gives some possibilities of having any bases at position 4 while still favoring G
#threshold score &gt;= 6.5 : keep all degenerated motif with degenerated = 6.5 and canonical = 7

SP1 = [
    [0, 0, 0, .5, 0, 0, 0],
    [0, 0, 0, .5, 0, 0, 0],
    [1, 1, 1, .5, 1, 1, 1],
    [0, 0, 0, 1., 0, 0, 0]]

#sixth element is TATA = TATAWAAG, length = 8 bases, size of matrix = 4*8, max score = 7.5
#threshold score &gt; 20 : keep all degenerated motif with degenerated = 20 and canonical = 23
TATA = [
    [0, 4, 0, 4, 4, 1, 1, 0],
    [4, 0, 4, 0, 4, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0]]

# define matrix and thresholds as tuple
motifs = dict(
    BREu = (BREu, 5.5),
    ETS = (ETS, 6),
    INR = (INR, 3.75),
    NFY = (NFY, 5),
    SP1 = (SP1, 6.25),
    TATA = (TATA, 20))

def get_motif_score(seq, motif):

    motif_size = len(motif[0])
    bases_encodings = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
    coded_seq = [
        bases_encodings[b] if b in 'ATGC' else None
        for b in seq]

    return [
        sum(
            0 if coded_b is None else motif[coded_b][i]
            for i, coded_b in enumerate(coded_seq[start:start+motif_size])
        ) for start in range(len(seq) - motif_size)]


def get_motifs(seq):
    &quot;&quot;&quot;find motif occurences in given sequence
    
    return list of (motif sequence, start, end, score)

    &quot;&quot;&quot;
    seq_size = len(seq)
    motif_hits = {}

    for motif_name, (motif_matrix, threshold) in motifs.items():
        motif_size = len(motif_matrix[0])
        rcseq = str(Seq(seq).reverse_complement())
        scores = get_motif_score(seq, motif_matrix)
        rc_scores = get_motif_score(rcseq, motif_matrix)
        
        motif_hits[motif_name] = [
            (seq[i:i+motif_size], i, i+motif_size, score,'+')
            for i, score in enumerate(scores)
            if score &gt;= threshold
        ] + [
            (rcseq[i:i+motif_size], seq_size-(i+motif_size), (seq_size-i), score,'-')
            for i, score in enumerate(rc_scores)
            if score &gt;= threshold
        ]

    return motif_hits
</code></pre>

CC BY-SA 4.0
2022-09-08T17:36:08.913Z
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Community
CC BY-SA 4.0
2022-09-09T13:17:35.363Z
What is the use and how we should interpret the matrix in your code ???
pippo1980
CC BY-SA 4.0
2022-09-09T13:25:32.217Z
Hi pippo1930, The matrices are self designed to map the motifs on basis of arbitary score, with the function of getting score, I am calculating the score on the window size equal to the length of the motif (length of matrix row), then if the score is above the threshold as indicated, it means that the motif exist. note that the sequence are non - consensus, which means that some basis are not very particular (that is why the score to map that particular base is less). I am not sure if that answer your question and thank you for taking interest to solve this problem.
Bioinfo_nerd
CC BY-SA 4.0
2022-09-09T13:28:29.803Z
Its a really hard alghiritm to understand, here usually you should provide a test input with the expected results vs actual ones to help people understand the issue
pippo1980
CC BY-SA 4.0
2022-09-09T13:34:36.460Z
How of curiosity could you show me how a BREu motif can get a > 5,5 score
pippo1980
CC BY-SA 4.0
2022-09-09T14:04:29.990Z
print(get_motifs('TATAWAAG')) gives {'BREu': [], 'ETS': [], 'INR': [], 'NFY': [], 'SP1': [], 'TATA': []}
pippo1980
CC BY-SA 4.0
2022-09-09T14:52:40.327Z
lower the threshold ???
pippo1980
CC BY-SA 4.0
2022-09-09T14:56:55.110Z
or change the matrix values
pippo1980
<python><phylogenetics>
I want to map TFBS on to the sequences from dataset containing 8448 (210 base) sequences. its not giving all the hits I need
# I want to map TFBS on to the sequences from dataset containing 8448 (210 base) sequences. its not giving all the hits I need

(Asked by: username_0 on 2022-09-08)

I want to map TFBS (transcription factor binding sites) on to the sequences from dataset containing 8448 (210 base) sequences, but it is  not giving all the hits I need
import pickle
from pickle import load, dump
from collections import defaultdict
from Bio.Seq import Seq

#first element is BREu = SSRCGCC, length = 7 bases, size of matrix = 4*7, max score = 3.33
# first row = A, second row = T, third row = G and forth row = C
BREu = [
    [0.0, 0.0, 0.5, 0, 0, 0, 0],
    [0.0, 0.0, 0.0, 0, 0, 0, 0],
    [0.5, 0.5, 0.5, 0, 1, 0, 0],
    [0.5, 0.5, 0.0, 1, 0, 1, 1]]

#second element is ETS = MGGAART, length = 7 bases, size of matrix = 4*7, max score = 6
ETS = [
    [0.5, 0, 0, 1, 1, 0.5, 0],
    [0.0, 0, 0, 0, 0, 0.0, 1],
    [0.0, 1, 1, 0, 0, 0.5, 0],
    [0.5, 0, 0, 0, 0, 0.0, 0]]

#Third element is INR = YYANWYY, length = 7 bases, size of matrix = 4*7, max score = 3.75
INR = [
    [0.0, 0.0, 1, 0.25, 0.5, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.5, 0.5, 0.5],
    [0.0, 0.0, 0, 0.25, 0.0, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.0, 0.5, 0.5]]

#forth element is NFY = CCAAT, length = 5 bases, size of matrix = 4*5, max score = 5
NFY = [
    [0, 0, 1, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0],
    [1, 1, 0, 0, 0]]

#seventh element is SP1_NC = GGGNGGG, length = 7 bases, size of matrix = 4*7, max score = 7
#Gives some possibilities of having any bases at position 4 while still favoring G
#threshold score >= 6.5 : keep all degenerated motif with degenerated = 6.5 and canonical = 7

SP1 = [
    [0, 0, 0, .5, 0, 0, 0],
    [0, 0, 0, .5, 0, 0, 0],
    [1, 1, 1, .5, 1, 1, 1],
    [0, 0, 0, 1., 0, 0, 0]]

#sixth element is TATA = TATAWAAG, length = 8 bases, size of matrix = 4*8, max score = 7.5
#threshold score > 20 : keep all degenerated motif with degenerated = 20 and canonical = 23
TATA = [
    [0, 4, 0, 4, 4, 1, 1, 0],
    [4, 0, 4, 0, 4, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0]]

# define matrix and thresholds as tuple
motifs = dict(
    BREu = (BREu, 5.5),
    ETS = (ETS, 6),
    INR = (INR, 3.75),
    NFY = (NFY, 5),
    SP1 = (SP1, 6.25),
    TATA = (TATA, 20))

def get_motif_score(seq, motif):

    motif_size = len(motif[0])
    bases_encodings = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
    coded_seq = [
        bases_encodings[b] if b in 'ATGC' else None
        for b in seq]

    return [
        sum(
            0 if coded_b is None else motif[coded_b][i]
            for i, coded_b in enumerate(coded_seq[start:start+motif_size])
        ) for start in range(len(seq) - motif_size)]


def get_motifs(seq):
    """find motif occurences in given sequence
    
    return list of (motif sequence, start, end, score)

    """
    seq_size = len(seq)
    motif_hits = {}

    for motif_name, (motif_matrix, threshold) in motifs.items():
        motif_size = len(motif_matrix[0])
        rcseq = str(Seq(seq).reverse_complement())
        scores = get_motif_score(seq, motif_matrix)
        rc_scores = get_motif_score(rcseq, motif_matrix)
        
        motif_hits[motif_name] = [
            (seq[i:i+motif_size], i, i+motif_size, score,'+')
            for i, score in enumerate(scores)
            if score >= threshold
        ] + [
            (rcseq[i:i+motif_size], seq_size-(i+motif_size), (seq_size-i), score,'-')
            for i, score in enumerate(rc_scores)
            if score >= threshold
        ]

    return motif_hits

**Tags:** python, phylogenetics
> **Comments:**
> **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
> **username_2**: What is the use and how we should interpret the matrix in your code ???  
> **username_0**: Hi pippo1930, The matrices are self designed to map the motifs on basis of arbitary score, with the function of getting score, I am calculating the score on the window size equal to the length of the motif (length of matrix row), then if the score is above the threshold as indicated, it means that the motif exist. note that the sequence are non - consensus, which means that some basis are not very particular (that is why the score to map that particular base is less). I am not sure if that answer your question and thank you for taking interest to solve this problem.  
> **username_2**: Its a really hard alghiritm to understand, here usually you should provide a test input with the expected results vs actual ones to help people understand the issue  
> **username_2**: How of curiosity could you show me how a BREu motif can get a > 5,5 score  
> **username_2**: print(get_motifs('TATAWAAG')) gives {'BREu': [], 'ETS': [], 'INR': [], 'NFY': [], 'SP1': [], 'TATA': []}  
> **username_2**: lower the threshold ???  
> **username_2**: or change the matrix values  

---

## 0 answers

