RBP motif extraction and visualization

RBP motif extraction and visualization#

BRIDGE aims to extract motifs from RBP datasets and generate visualization graphics. BRIDGE uses a series of scripts and tools to achieve this, including FASTA file processing, sequence-structure Position Weight Matrix (PWM) generation and final motif logo visualization.

Overview of step-by-step instructions#

├── discover_motif   # Discover motifs from attention maps
├── convert_format   # Convert .txt files to MEME format files
├── tomtom           # Compare motifs to mCrossBase database motifs using Tomtom
├── struct_motif     # Automatically extract motif-related sequences
├── seq_struct2pwm   # Generate combined sequence-structure PWM matrices

Prerequisites#

Ensure all required R packages and Python libraries are installed.
Modify file paths and file names in the scripts to match your specific dataset.

Dependencies#

R Packages

dendsort
dendextend
circlize
ComplexHeatmap
grid
ggplot2
ggseqlogo
ggsci
gridExtra
seqLogo
Cairo

Python Libraries

argparse
pathlib
numpy
Bio (Biopython)

Important

The steps below support large-scale dataset analysis. Here we use two datasets as examples; for the remaining datasets, please download them from Figshare first and place them under the corresponding directories.

Discover motifs from attention maps#

Download & prepare paper-scale inputs (Figshare)

To reproduce the paper-scale motif extraction without re-running the model, you need two Figshare archives:

Datasets archive

Contains the raw sequence sets for each RBP×cell dataset (e.g., pos.fa / neg.fa).
These sequences define what biological instances we are mining motifs from.

Attention archive

Contains precomputed model outputs for each dataset:
- atten.npy: per-sequence, per-position attention scores
- pred_results.npy: per-sequence probability of positive class
These files are the direct inputs to downstream motif discovery (so you can run motif extraction via np.load(...)).

Where to place them

Download both archives from Figshare and extract them into:

motif_construction/datasets/

After extraction, the expected layout is:

motif_construction/datasets/
  `<DATASET_ID>/`
    `dev.tsv`
    `atten.npy`
    `pred_results.npy`

<DATASET_ID> is the dataset key (e.g., AARS_K562).
Downstream scripts/notebooks will iterate over datasets/*/ and load:
- sequences from dev.tsv
- attention outputs from atten.npy and pred_results.npy

Next, follow the instructions below to run motif discovery on the precomputed attention maps.

%%bash
set -euo pipefail

Motif_Construction_ROOT=$(cd "$(dirname "$0")/../../../../motif_construction"; pwd)
export KMER=1
echo $(pwd)
export LD_LIBRARY_PATH="$CONDA_PREFIX/lib:$LD_LIBRARY_PATH"
cd $Motif_Construction_ROOT

# Define the base directory
export BASE_DIR="datasets/"

# Loop through each folder in the BASE_DIR
for filename in $BASE_DIR*/; do
    if [ -d "$filename" ]; then
        echo "Processing folder: $filename"

        # Set the paths for the current folder
        export DATA_PATH=$filename
        export ATTENTION_PATH=$filename
        export OUTPUT_PATH=$filename/result/
        # Create the motif result directory if it doesn't exist
        mkdir -p $OUTPUT_PATH

        # Run the motif finding script for the current folder
        python motif/Discovery_motifs.py \
            --input_rna_dir $DATA_PATH \
            --attention_dir $ATTENTION_PATH \
            --motif_length 6 \
            --min_region_len 6 \
            --fdr_cutoff 0.05 \
            --min_motif_count 1 \
            --align_all_ties \
            --output_motif_dir $OUTPUT_PATH \
            --verbose
    fi
done

/fs1/private/user/wangyubo/code/BRIDGE/docs/tutorials/notebooks/motif_analysis
Processing folder: datasets/AARS_K562/

=== Motif Analysis ===
 * Positives: 1000
 * Negatives: 2000
 * Raw motifs found: 1201
 * Applying statistical filtering …
[+] Motif AGGGGG: N=3000; K=1000; n=114; x=66; p=3.60e-08
[+] Motif AGGAAA: N=3000; K=1000; n=135; x=72; p=7.73e-07
[+] Motif TGAAGA: N=3000; K=1000; n=243; x=116; p=8.53e-07
[+] Motif AAAAGA: N=3000; K=1000; n=117; x=66; p=1.44e-07
 * Motif sequences saved to → datasets/AARS_K562//result/

>>> Final motifs (top): ['AGGGGG', 'TGGGGA', 'CAAAAC', 'AGGAAA', 'TGAAGA', 'AAAAGA', 'GAGGATGAA', 'CCCAAGGACT']
Processing folder: datasets/AATF_K562/

=== Motif Analysis ===
 * Positives: 1000
 * Negatives: 2000
 * Raw motifs found: 1202
 * Applying statistical filtering …
[+] Motif GGGCGGG: N=3000; K=1000; n=56; x=35; p=6.16e-06
[+] Motif TGCGGC: N=3000; K=1000; n=130; x=67; p=9.12e-06
[+] Motif GGGGGC: N=3000; K=1000; n=184; x=99; p=2.92e-09
[+] Motif CGGGTT: N=3000; K=1000; n=28; x=21; p=7.10e-06
[+] Motif GGCCGTG: N=3000; K=1000; n=57; x=37; p=8.52e-07
[+] Motif GTGCGG: N=3000; K=1000; n=98; x=57; p=2.58e-07
[+] Motif CGCTGC: N=3000; K=1000; n=179; x=87; p=9.28e-06
[+] Motif TGGTGGG: N=3000; K=1000; n=68; x=41; p=3.71e-06
[+] Motif CCGGGG: N=3000; K=1000; n=156; x=84; p=4.89e-08
[+] Motif GCGGCC: N=3000; K=1000; n=189; x=103; p=5.10e-10
[+] Motif GGCCGCG: N=3000; K=1000; n=61; x=44; p=4.45e-10
[+] Motif GGGCGG: N=3000; K=1000; n=162; x=98; p=3.04e-13
[+] Motif GGGGGT: N=3000; K=1000; n=90; x=52; p=1.18e-06
[+] Motif TGGCGG: N=3000; K=1000; n=136; x=70; p=6.06e-06
[+] Motif GCGGGGA: N=3000; K=1000; n=41; x=30; p=1.76e-07
[+] Motif GGCGGA: N=3000; K=1000; n=110; x=64; p=4.51e-08
[+] Motif GGTCGC: N=3000; K=1000; n=60; x=37; p=5.27e-06
[+] Motif GGGCGC: N=3000; K=1000; n=131; x=76; p=2.72e-09
[+] Motif GGTGAG: N=3000; K=1000; n=100; x=59; p=7.76e-08
[+] Motif GGGCTG: N=3000; K=1000; n=254; x=122; p=2.78e-07
[+] Motif CGGGGC: N=3000; K=1000; n=161; x=87; p=2.21e-08
[+] Motif GGTGGG: N=3000; K=1000; n=206; x=112; p=9.87e-11
[+] Motif CTGGGG: N=3000; K=1000; n=247; x=119; p=3.23e-07
[+] Motif CGGGGA: N=3000; K=1000; n=122; x=68; p=1.74e-07
[+] Motif GCGGGG: N=3000; K=1000; n=153; x=92; p=2.87e-12
[+] Motif GTGGGG: N=3000; K=1000; n=186; x=92; p=1.93e-06
[+] Motif CGGGCC: N=3000; K=1000; n=167; x=89; p=3.65e-08
[+] Motif CGCCGC: N=3000; K=1000; n=154; x=84; p=2.21e-08
[+] Motif GGCCGG: N=3000; K=1000; n=192; x=110; p=1.61e-12
[+] Motif AGGCCGG: N=3000; K=1000; n=40; x=28; p=2.01e-06
[+] Motif GGCGGCGG: N=3000; K=1000; n=66; x=47; p=2.27e-10
 * Motif sequences saved to → datasets/AATF_K562//result/

>>> Final motifs (top): ['TGGTGG', 'CGGGAA', 'TGCGGC', 'GGGGGC', 'GCGGAC', 'TGGGGT', 'CGGGTT', 'CGCTGC', 'TGGGCC', 'CCGGGG']

Export motifs to MEME format#

import os
import numpy as np

alphabet = ['A', 'C', 'G', 'U']
letter_to_idx = {letter: idx for idx, letter in enumerate(alphabet)}

def txt_to_meme(input_path, output_path):
    '''
    Convert a.txt file to a MEME format file.
    :param input_path: Path to the input.txt file.
    :param output_path: Path to the output.meme file.
    '''
    with open(input_path, 'r') as file:
        sequences = [line.strip().replace('T', 'U') for line in file.readlines()]

    # Initialise a matrix to store the probability for each position.
    num_sequences = len(sequences)
    sequence_length = len(sequences[0])
    probability_matrix = np.zeros((sequence_length, len(alphabet)))

    # Compute the probability matrix.
    for seq in sequences:
        for i, base in enumerate(seq):
            probability_matrix[i][letter_to_idx[base]] += 1

    # Convert the frequency to probability (probability of each position is the number of times the base appears / total number of sequences).
    probability_matrix /= num_sequences

    # Write the probability matrix to a MEME file.
    with open(output_path, 'w') as f:
        # MEME file header
        f.write("MEME version 5\n")
        f.write("ALPHABET= ACGU\n")
        f.write("Background letter frequencies: A 0.25 C 0.25 G 0.25 U 0.25\n")
        motif_name = os.path.basename(input_path).replace('.txt', '')  # Use the file name as the MOTIF name.
        f.write(f"MOTIF {motif_name}\n")
        f.write(f"letter-probability matrix: alength= 4 w= {sequence_length}\n")
        
        # Write the probability matrix.
        for i in range(sequence_length):
            f.write(" ".join([f"{prob:.14f}" for prob in probability_matrix[i]]) + "\n")

    print(f"Meme file written to {output_path}")

def process_directory(base_directory):
    # Process all directories under the base directory.
    for root, dirs, files in os.walk(base_directory):
        if 'result' in root:  # Only process directories containing 'result'.
            for file in files:
                if file.endswith('.txt'):  # Only process .txt files.
                    input_file = os.path.join(root, file)
                    output_file = os.path.join(root, file.replace('.txt', '.meme'))  # Output as .meme file.
                    txt_to_meme(input_file, output_file)

# Main program, set the base directory to be processed.
base_directory = "../../../../motif_construction/datasets"
process_directory(base_directory)

Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_AGGAAA_1.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_GAGGATGAA_1.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_CCCAAGGACT_1.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_AGGGGG_1.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_TGAAGA_2.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_TGGGGA_4.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_CAAAAC_1.meme
Meme file written to ../../../../motif_construction/datasets/AARS_K562/result/motif_AAAAGA_1.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CGGCCC_3.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GCGTGG_9.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TGGGCC_9.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CGGGAA_6.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TTCGCG_2.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TAGGGC_4.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGCCGTG_7.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGGGGC_19.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CGGGTT_3.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TGGTGG_16.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGGTCC_2.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GTCTGG_1.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGGCGG_16.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TGGCGC_5.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GCGGAC_3.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CCTGCGC_3.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGTCGC_2.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGTTGGG_9.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CCGGGG_11.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_GGTGAG_5.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TGCGGC_6.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_TGGGGT_10.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CGCTGC_2.meme
Meme file written to ../../../../motif_construction/datasets/AATF_K562/result/motif_CAGGTAA_1.meme

Compare motifs against mCrossBase and ATtRACT using Tomtom#

This section runs Tomtom (from the MEME suite) to match BRIDGE-discovered motifs against two external motif databases:

ATtRACT (RBP motif collection)
mCrossBase (cell/RBP-specific motif collection)

The goal is to produce per-motif similarity reports (Tomtom outputs) that can later be used for:

validating that discovered motifs resemble known motifs,
finding database “hits” for each discovered motif,
enabling downstream intersection / selection logic for paper figures.

What this script does

For each database × RBP × cell line:

Look for BRIDGE motif results in: motif_construction/datasets/<RBP>_<CELL>/result/*.meme
For every discovered motif file *.meme, run: tomtom <DB_MOTIF_FILE> <DISCOVERED_MOTIF_FILE>
Save Tomtom outputs to a dedicated folder: motif_construction/datasets/<RBP>_<CELL>_<tmp_suffix>/<meme_file>_tomtom_results/

Prerequisites

Install Tomtom (MEME suite)

You need the tomtom executable available in your PATH:

conda install -c conda-forge -c bioconda meme

Verify: tomtom -h

Inputs#

A) BRIDGE discovered motif files (query motifs) For each dataset: motif_construction/datasets/<RBP>_<CELL>/result/*.meme

Each *.meme file is a set of motifs discovered by BRIDGE and written in MEME motif format.

B) Database motif files (target motifs) Depending on the database:

ATtRACT: motif_construction/db/ATtRACT/<RBP>.meme
mCrossBase: motif_construction/db/mCrossBase/<CELL>.<RBP>.meme

Outputs#

For each dataset and each query MEME file, Tomtom writes a result directory:

ATtRACT example: motif_construction/datasets/HNRNPU_K562_tmp_results/<query_file>.meme_tomtom_results/

mCrossBase example: motif_construction/datasets/HNRNPU_K562_tmp_results_mcross/<query_file>.meme_tomtom_results/

Tomtom typically generates files like:

tomtom.tsv (main match table)
tomtom.xml
tomtom.html
PNG motif alignment images

import os
import subprocess
from pathlib import Path

# List of cell lines
cell_lines = ['K562', 'HepG2', 'HEK293', 'H9', 'HEK293T', 'HeLa']

# Base directory containing example results
BRIDGE_HOME = Path.cwd().resolve().parents[3]
BRIDGE_HOME = str(BRIDGE_HOME) + "/"
base_directory = os.path.join(BRIDGE_HOME,'motif_construction/datasets')
# Directory containing input FASTA files (used to derive RBP names)

# Database configurations: ATtRACT and mCrossBase. You need to download these sources from our Figshare page.
databases = {
    'ATtRACT': {
        'path': BRIDGE_HOME +'motif_construction/db/ATtRACT',
        'tmp_suffix': '_tmp_results',
        'min_overlap': '5',
        'db_file_pattern': '{rbp}.meme'
    },
    'mCrossBase': {
        'path': BRIDGE_HOME +'motif_construction/db/mCrossBase',
        'tmp_suffix': '_tmp_results_mcross',
        'min_overlap': '6',
        'db_file_pattern': '{cell}.{rbp}.meme'
    }
}

# Get all RBP names
rbp_files = [d for d in os.listdir(base_directory) if os.path.isdir(os.path.join(base_directory, d)) and "_" in d]
rbps = set(f.split('_')[0] for f in rbp_files)

# Iterate over databases, RBPs, and cell lines
for db_name, db_conf in databases.items():
    meme_directory = db_conf['path']
    print(f"\n=== Running Tomtom for database: {db_name} ===")

    for rbp in rbps:
        for cell in cell_lines:
            # Construct the path to the result folder
            result_directory = os.path.join(base_directory, f"{rbp}_{cell}", "result")

            # Skip if the result folder does not exist
            if not os.path.exists(result_directory):
                # print(f"Skipping {rbp}_{cell}: result folder does not exist.")
                continue

            # Temporary directory for Tomtom intermediate results
            tmp_dir = os.path.join(base_directory, f"{rbp}_{cell}{db_conf['tmp_suffix']}")
            os.makedirs(tmp_dir, exist_ok=True)

            # Iterate through all MEME files in the result folder
            for file in os.listdir(result_directory):
                if not file.endswith('.meme'):
                    continue

                meme_file = os.path.join(result_directory, file)

                # Output directory for Tomtom results
                out_directory = os.path.join(tmp_dir, f"{file}_tomtom_results")
                os.makedirs(out_directory, exist_ok=True)

                # Determine the corresponding motif database file
                db_meme = db_conf['db_file_pattern'].format(rbp=rbp, cell=cell)
                db_meme_path = os.path.join(meme_directory, db_meme)

                # Skip if the database motif file is missing
                if not os.path.exists(db_meme_path):
                    print(f"Database file not found: {db_meme_path}, skipping...")
                    continue

                # Construct the Tomtom command
                tomtom_cmd = [
                    'tomtom', '-no-ssc', '-oc', out_directory,
                    '-verbosity', '1',
                    '-min-overlap', db_conf['min_overlap'],
                    '-dist', 'pearson', '-evalue', '-thresh', '0.05',
                    '-time', '300', '-png',
                    db_meme_path,
                    meme_file
                ]

                # Execute Tomtom
                print(f"Running Tomtom for {file} in {cell} ({db_name})...")
                subprocess.run(tomtom_cmd)

=== Running Tomtom for database: ATtRACT ===
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/ATtRACT/AARS.meme, skipping...

=== Running Tomtom for database: mCrossBase ===
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Database file not found: /fs1/private/user/wangyubo/code/BRIDGE/motif_construction/db/mCrossBase/K562.AATF.meme, skipping...
Running Tomtom for motif_TGAAGA_2.meme in K562 (mCrossBase)...
Running Tomtom for motif_AGGGGG_1.meme in K562 (mCrossBase)...
Running Tomtom for motif_CAAAAC_1.meme in K562 (mCrossBase)...
Running Tomtom for motif_AGGAAA_1.meme in K562 (mCrossBase)...
Running Tomtom for motif_CCCAAGGACT_1.meme in K562 (mCrossBase)...
Running Tomtom for motif_GAGGATGAA_1.meme in K562 (mCrossBase)...
Running Tomtom for motif_TGGGGA_4.meme in K562 (mCrossBase)...
Running Tomtom for motif_AAAAGA_1.meme in K562 (mCrossBase)...

Extract related structure motifs#

It iterates through all subdirectories under the given base path, each representing one RBP-cell combination. For every .txt motif file found in the result/ folder, it performs the following steps:

Parse motif names and numeric indices from filenames (e.g., motif_TGCGGC_3.txt).
Apply filtering criteria:
- The motif sequence (the middle part) must be exactly 6 letters long.
- The final numeric suffix must be greater than 2 (optional, can be enabled).
Match each motif against the corresponding annotation files (utils/motif_prior/output/{rbp_cell_folder}/output/fg_alphabet_annotations.tab), which are saved in the figshare repository and contain the aligned structures for each motif.
Extract the aligned six-letter sequence from the same position as the motif match.
Write the extracted structures into separate FASTA files, ensuring that each motif has its own dedicated output folder to avoid overwriting.

Output files are saved in: ../../../../motif_construction/datasets/Struct_motifs/{RBP_CELL}/{MOTIF}_{INDEX}/

import re
import os

# Base directory containing subdirectories for each RBP_CELL
BRIDGE_ROOT="../../../../"
base_dir = BRIDGE_ROOT+'motif_construction/datasets/'
print(base_dir)
# Iterate through all subdirectories (each representing an RBP_CELL)
for rbp_cell_folder in os.listdir(base_dir):
    rbp_cell_path = os.path.join(base_dir, rbp_cell_folder)
    
    # Proceed only if it is a directory
    if os.path.isdir(rbp_cell_path):
        result_dir = os.path.join(rbp_cell_path, 'result')
        
        # Ensure the result directory exists
        if os.path.exists(result_dir):
            # Iterate through all .txt files inside the "result" directory
            for motif_file_name in os.listdir(result_dir):
                if motif_file_name.endswith('.txt'):
                    motif_file_path = os.path.join(result_dir, motif_file_name)
                    
                    # Extract motif sequence name (e.g., motif_TGGGAG_3.txt → TGGGAG)
                    motif_name_match = re.match(r'motif_([A-Z]+)_\d+\.txt', motif_file_name)
                    if motif_name_match:
                        motif_name = motif_name_match.group(1)
                        # Replace T with U to convert DNA to RNA format
                        motif_name = motif_name.replace('T', 'U')
                        
                        # Split the motif filename to extract numeric and sequence parts
                        motif_parts = motif_file_name.split('_')
                        last_num = int(motif_parts[-1].replace('.txt', ''))  # Last numeric part of the filename
                        
                        # This section filters motif files based on naming conventions
                        # Conditions:
                        #   (1) The last number must be > 2
                        #   (2) The motif sequence must be exactly 6 letters long
                        if len(motif_parts) >= 3:
                            try:
                                # Uncomment below if you want to strictly enforce last_num > 2
                                # if last_num <= 2:
                                #     continue
                                
                                if len(motif_name) != 6:
                                    continue  # Skip motifs that are not 6 letters long
                            except ValueError:
                                continue

                        # Construct path to the annotation file for this RBP_CELL
                        annotations_file = BRIDGE_ROOT+f'utils/motif_prior/output/{rbp_cell_folder}/output/fg_alphabet_annotations.tab'

                        # Proceed only if the annotation file exists
                        if os.path.exists(annotations_file):
                            with open(annotations_file, 'r') as annotations_file_obj:
                                annotations_lines = annotations_file_obj.readlines()

                            # Iterate through each line in the annotation file
                            for line in annotations_lines:
                                fields = line.strip().split('\t')
                                
                                # Skip malformed lines
                                if len(fields) < 2:
                                    continue

                                sequence = fields[1]      # Full RNA sequence
                                struct7  = fields[-1]     # 7-state structural alphabet aligned to `sequence` (E, L, H, R, B, M, T)

                                # Search for the motif sequence within the second field
                                match = re.search(motif_name, sequence)
                                if match:
                                    # Determine the position of the motif within the sequence
                                    position = match.start()
                                    
                                    # Aligned 6-char structure motif
                                    six_letters = struct7[position:position+6]
                                    
                                    # Create a FASTA header line
                                    fasta_header = f">{rbp_cell_folder}_motif_{motif_name}_pos_{position+1}"
                                    fasta_sequence = six_letters
                                    
                                    # Define output directory (unique per motif to avoid overwriting)
                                    output_dir = BRIDGE_ROOT+f"motif_construction/datasets/Struct_motifs/{rbp_cell_folder}/{motif_name}_{last_num}"
                                    os.makedirs(output_dir, exist_ok=True)
                                    
                                    # Define output FASTA file path
                                    output_fa_path = f"{output_dir}/{rbp_cell_folder}_{motif_name}_{last_num}.fa"
                                    
                                    # Append the FASTA entry to the output file
                                    with open(output_fa_path, 'a') as fa_file:
                                        fa_file.write(f"{fasta_header}\n{fasta_sequence}\n")
                            print(f"[Processed] {rbp_cell_folder}")
                        else:
                            print(f"Warning: Annotation file not found for {rbp_cell_folder}: {annotations_file}")
        else:
            print(f"Warning: Result directory not found for {rbp_cell_folder}")

../../../../motif_construction/datasets/
Warning: Result directory not found for AATF_K562_tmp_results
Warning: Result directory not found for Struct_motifs
Warning: Result directory not found for AARS_K562_tmp_results
[Processed] AARS_K562
[Processed] AARS_K562
[Processed] AARS_K562
[Processed] AARS_K562
[Processed] AARS_K562
[Processed] AARS_K562
Warning: Result directory not found for AATF_K562_tmp_results_mcross
Warning: Result directory not found for AARS_K562_tmp_results_mcross
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab
Warning: Annotation file not found for AATF_K562: ../../../../utils/motif_prior/output/AATF_K562/output/fg_alphabet_annotations.tab

Generate sequence-structure PWMs#

For each RBP-cell combination directory under the base path, it:

Searches for all FASTA files that contain RNA secondary structure motifs.
Matches each structural motif with its corresponding sequence motif file (motif_*.txt).
Converts the paired data into an 11*N PWM matrix, where:
- The first 4 rows represent nucleotide frequencies (A, C, G, U).
- The last 7 rows represent structural element frequencies (L, R, E, H, B, M, T).
- N is the motif length (typically 6).
Normalizes both parts and writes the resulting PWM matrix to an output file specific to each RBP-cell.

Output directory:

`motif_construction/datasets/Seq_struct_motif/{RBP_CELL}.txt`

import os
import numpy as np
from Bio import SeqIO


#  PWM Construction Functions
def seq_to_pwm(sequences, motif_len, structure_data):
    """
    Convert a list of RNA sequences and structural annotations into an 11*N PWM matrix.

    Parameters:
        sequences (list of str): List of RNA sequences (ACGU).
        motif_len (int): Length of the motif.
        structure_data (list of str): List of structural annotation strings (L, R, E, H, B, M, T).

    Returns:
        np.ndarray: PWM matrix with 11 rows (4 for ACGU + 7 for structure) and motif_len columns.
    """
    pwm = np.zeros((11, motif_len))  # 11 rows: 4 (ACGU) + 7 (structure)
    bases = ['A', 'C', 'G', 'T']     # Note: 'T' is used for compatibility with motif files
    structures = ['L', 'R', 'E', 'H', 'B', 'M', 'T']

    # Process sequence (ACGU) part
    for seq in sequences:
        for i, base in enumerate(seq):
            if base in bases:
                pwm[bases.index(base), i] += 1

    # Normalize nucleotide part
    pwm[:4, :] /= len(sequences)

    # Process structure part
    for struct in structure_data:
        for i, s in enumerate(struct):
            if s in structures:
                pwm[4 + structures.index(s), i] += 1

    # Normalize structure part
    pwm[4:, :] /= len(structure_data)

    return pwm


def save_pwm_to_file(pwm, output_file):
    """
    Save a PWM matrix to a text file.
    """
    np.savetxt(output_file, pwm, fmt='%.6f', delimiter="\t")


def process_fasta_and_motif(fasta_file, motif_file, output_file, motif_name):
    """
    Process one pair of FASTA and motif files to generate a PWM matrix.

    Parameters:
        fasta_file (str): Path to the FASTA file containing structural motifs.
        motif_file (str): Path to the motif file containing sequence motifs.
        output_file (str): File to save the resulting PWM matrix.
        motif_name (str): Name of the motif being processed.
    """
    # Load structure strings from the FASTA file
    structures = [str(record.seq) for record in SeqIO.parse(fasta_file, "fasta")]

    # Load sequences from the motif file
    with open(motif_file, 'r') as f:
        sequences = [line.strip() for line in f.readlines()]

    # Assume all motifs have the same length
    motif_len = len(sequences[0])

    # Generate PWM
    pwm = seq_to_pwm(sequences, motif_len, structures)

    # Write PWM matrix to the output file
    with open(output_file, 'a') as out_f:
        out_f.write(f"{motif_name}\n")
        np.savetxt(out_f, pwm, fmt='%.6f', delimiter="\t")
        out_f.write('\n')


def find_and_process_files(base_dir, result_dir, output_file):
    """
    Search for all FASTA files in the given directory and process each one
    along with its corresponding motif file to produce PWM matrices.

    Parameters:
        base_dir (str): Directory containing structural FASTA files.
        result_dir (str): Directory containing sequence motif files.
        output_file (str): Output file for the combined PWM results.
    """
    for root, _, files in os.walk(base_dir):
        fasta_files = [f for f in files if f.endswith(".fa")]

        for fasta_file in fasta_files:
            # Extract motif name from the FASTA filename, e.g., AARS_K562_CAAAAC_1.fa
            motif_name = os.path.basename(fasta_file).split('.')[0]
            motif = motif_name.split('_')[2] + '_' + motif_name.split('_')[3]
            motif = motif.replace('U', 'T')  # Ensure compatibility with motif file names

            motif_file = os.path.join(result_dir, f"motif_{motif}.txt")

            # Only process if the motif file exists
            if os.path.exists(motif_file):
                fasta_file_path = os.path.join(root, fasta_file)
                process_fasta_and_motif(fasta_file_path, motif_file, output_file, motif_name)


def process_rbp_cells(base_parent_dir, result_parent_dir, output_parent_dir):
    """
    Automatically detect all RBP-cell directories under the base parent directory
    and generate PWM matrices for each one.

    Parameters:
        base_parent_dir (str): Base directory containing subdirectories for each RBP-cell.
        result_parent_dir (str): Parent directory containing motif result folders.
        output_parent_dir (str): Directory where output files will be stored.
    """
    # Automatically get all RBP-cell folders
    rbp_cells = [
        d for d in os.listdir(result_parent_dir)
        if os.path.isdir(os.path.join(result_parent_dir, d))
    ]

    for rbp_cell in rbp_cells:
        base_dir = os.path.join(base_parent_dir, rbp_cell)
        result_dir = os.path.join(result_parent_dir, rbp_cell, "result")
        os.makedirs(output_parent_dir, exist_ok=True)
        output_file = os.path.join(output_parent_dir, f"{rbp_cell}.txt")

        print(f"Processing {rbp_cell}...")

        # Remove existing output file to start fresh
        if os.path.exists(output_file):
            os.remove(output_file)

        find_and_process_files(base_dir, result_dir, output_file)

        print(f"PWM matrix generation completed for {rbp_cell}.")

#  Main Script Execution
BRIDGE_ROOT="../../../../"
base_parent_dir = BRIDGE_ROOT+"motif_construction/datasets/Struct_motifs"
result_parent_dir = BRIDGE_ROOT+"motif_construction/datasets"
output_parent_dir = BRIDGE_ROOT+"motif_construction/datasets/Seq_struct_motif"

# Start processing all detected RBP-cell directories
process_rbp_cells(base_parent_dir, result_parent_dir, output_parent_dir)

print("PWM matrix generation completed for all RBP-cell datasets.")

Processing AATF_K562_tmp_results...
PWM matrix generation completed for AATF_K562_tmp_results.
Processing Struct_motifs...
PWM matrix generation completed for Struct_motifs.
Processing AARS_K562_tmp_results...
PWM matrix generation completed for AARS_K562_tmp_results.
Processing AARS_K562...
PWM matrix generation completed for AARS_K562.
Processing AATF_K562_tmp_results_mcross...
PWM matrix generation completed for AATF_K562_tmp_results_mcross.
Processing AARS_K562_tmp_results_mcross...
PWM matrix generation completed for AARS_K562_tmp_results_mcross.
Processing AATF_K562...
PWM matrix generation completed for AATF_K562.
PWM matrix generation completed for all RBP-cell datasets.