IBIS PWM & AAA Submissions - Method Write-up

No models, just interpretable rules

Team name

callitmagic

Primary Disciplines

A2G-PWM, G2A-PWM, A2G-AAA, G2A-AAA

Summary

This submission focuses on using interpretable rules and motif-based predictions for transcription factor binding site identification. I used tools from the MEME Suite (including MEME - Multiple EM for Motif Elicitation - for motif discovery, AME - Analysis of Motif Enrichment - for PWM selection based on enrichment analysis and FIMO - Find Individual Motif Occurrences - for scanning sequences for selected motifs) to generate position weight matrices (PWMs) and identify motif occurrences in genomic and artificial data. Some other tools were also tested (see in the Methods section), but not included in the final submission.

My approach utilizes predefined parameter grids for motif generation, cleaning the results based on e-value thresholds, and using top-performing PWMs to make predictions for the IBIS AAA track. I emphasize interpretability, ensuring that each motif identified is linked back, without relying on complex machine-learning models. My pipeline also includes steps to handle genomic and artificial datasets differently, ensuring robustness across different data types.

Implementation & Software

• PWM generation and motif analysis: MEME Suite (MEME, FIMO, AME), bash
• Data handling: R (with packages such as Biostrings, dplyr, purrr, and XML)
• Code repository: https://github.com/useitforgood/ibis_final (requires access, complete script)

For the final PWMs: No additional models were used. A greedy search was conducted among numerous motifs generated by MEME on differently pre-processed and aggregated data. The decision rule for top-4 motifs selection was based on AME results using control sequences from differently pre-processed and aggregated data (the exact decision rule was selected based on LB results).

For the final AAAs: No additional models were used in the final submission. Scorings were made using FIMO for the extended list of best-selected motifs according to the decision rule, followed by aggregation of these to make predictions (the aggregation rule was selected based on LB results).

This scheme has given better results than ml models built on k-mers of different lengths, especially for underrepresented TFs, and than CNN for genome data.

Methods and steps

Full script Rmd:

https://github.com/useitforgood/ibis_final/blob/master/selected_for_final_submission.Rmd

Full script knitted html:

https://github.com/useitforgood/ibis_final/blob/master/selected_for_final_submission.html

Steps:
Data Preparation:

• Downloaded genomic data (hg38) and processed raw IBIS data files for each discipline (CHS, GHTS, PBM, SMS, HTS).
  • Resulting raw data processing:
    https://github.com/useitforgood/ibis_final/blob/master/functions/convert_data.R
    * for genome data:
    get WINDOW = 60 nucleotides before and after the absolute peak summit position
    * for artificial data:
    leave as-is
  • Tested, but not included:
    * for genome data:
    get the whole sequence from START to END position (required more time to calculate motifs with no significant accuracy gain and accuracy loss for some TFs)
    * for artificial data:
    several additional steps were tested, including CD-HIT for HT-Selex data clustering, fastp for quality trimming with adapter cleaning, and adapter concatenation (all these did not result in a noticeable improvement in the lb predictions and were not worth the additional complexity and runtime they introduced).
• Split datasets into genome-based and artificial data.

PWMs/motifs Generation:

• Defined parameter grids for motif discovery using the MEME tool, tuning for variables such as motif length, the usage of GHTS/CHS/other data, trimming, reverse complement searching, etc. This tuning felt necessary because the quality and binding characteristics of the TFs could vary, and different sequences could respond differently to specific methods. Therefore, it was not clear in advance how best to trim or filter sequences, which made systematic exploration essential.
  
  Here’s an explanation of each parameter used in the grid for genome data for motif generation:

1. length:
   Specifies the number of nucleotides to extract from the center of each sequence for motif discovery. If a value is provided, the sequence is trimmed to this length, centered around the midpoint of the original sequence.

• N: Extracts N nucleotides before and after the center.
• NA: Uses the full sequence (WINDOW*2 b.p.) without trimming.

1. GHTS:
   Indicates whether Genomic High-Throughput SELEX (GHTS) data should be included in the motif discovery:
   • TRUE: GHTS data is used.
   • FALSE: GHTS data is not used.
2. CHS:
   Indicates whether ChIP-Seq (CHS) data should have been included in the motif discovery:
   • TRUE: CHS data is used.
   • FALSE: CHS data is not used.
3. meme_rev:
   Consider both the given strand and the reverse complement strand when searching for motifs in a complementable alphabet (https://meme-suite.org/meme/doc/meme.html)
   • TRUE: Reverse complement is included.
   • FALSE: Reverse complement is not included.
4. pileup:
   Filter sequences based on the pileup (coverage) of a region. The percentile (25%, 50%, 75%, or 90%) is used to retain only the sequences that have a pileup greater than the specified threshold:

• 25, 50, 75, 90: Retains sequences with a pileup value above the given percentile.
• NA: No pileup filtering is applied.

1. maxw:
   Specifies the maximum width of the motifs to be discovered:
   • 15: Motif width is limited to 15 bases.
   • 30: Motif width can extend up to 30 bases.
2. seed:
   Sets the seed for random number generation to ensure reproducibility in the analysis:
   • 53: The fixed seed value used for all experiments.

Here’s an explanation of each parameter used in the grid for artificial data for motif generation:

1. SMS:
   Indicates whether SMS data should be included in the motif discovery:
   • TRUE: SMS data is used.
   • FALSE: SMS data is not used.
     NB! Actually, SMS was not used for final submission, it was always FALSE.
2. PBM:
   Indicates whether Protein Binding Microarray (PBM) data should be included in the motif discovery:
   • TRUE: PBM data is used.
   • FALSE: PBM data is not used.
3. HTS:
   Indicates whether High-Throughput Sequencing (HTS) data should be included:
   • TRUE: HTS data is used.
   • FALSE: HTS data is not used.
4. PBM_filter:
   Specifies the filtering method for PBM data, determining which sequences to retain based on signal intensity:
   • "max": Retains sequences with the maximum signal intensity for each probe.
   • "2sd": Retains sequences whose signal intensity is more than two standard deviations above the mean.
   • NA: No PBM filtering is applied.
5. PBM_subset:
   Determines if a specific subset of PBM data should be used:
   • "QNZS": Uses the QNZS subset of PBM data.
   • "SD": Uses the SD subset of PBM data.
   • NA: No PBM subset is selected, and all PBM data is used.
6. HTS_subset:
   Specifies whether a particular subset of HTS data should be used, filtering the data based on rounds of selection:
   • "C4": Uses the C4 subset of HTS data.
   • "C3": Uses both C3 and C4 subsets of HTS data.
   • NA: No HTS subset is selected, and all HTS data is used.
7. meme_rev:
   Consider both the given strand and the reverse complement strand when searching for motifs in a complementable alphabet (https://meme-suite.org/meme/doc/meme.html)
   • TRUE: Reverse complement is included.
   • FALSE: Reverse complement is not included.
8. seed:
   Sets the seed for random number generation to ensure reproducibility in the analysis:
   • 53: The fixed seed value used for all experiments.

• Processed the genome and artificial data separately to generate 6 motifs from training data per grid row.
• Cleaned motifs based on e-value thresholds, keeping only motifs with e-value <= 0.05.
  NB! The last two steps could have been combined, but memesuite ignores -nmotifs when e-value is set.
  
  **Tested, but not included:
  **objfun, shuf parameters, as significant gain on lb was only for 1 TF.
  
  The pre-calculated results of MEME processing are stored in the GitHub repository under the folders:
  meme_comparison/data/tmp_memedir_*
  If you want to recalculate and check the results for a specific grid row, simply save the files, delete the respective folder or subfolder, and compare the new file and the previous one.

Motif Enrichment Analysis (AME):

• Generated control sequences from CHS, GHTS, PBM, SMS, and HTS datasets to run AME comparisons for each TF.
  • For PBM data, signal intensities were calculated, and sequences with signal intensity above a threshold (mean + 2 standard deviations) were selected.
  • HTS sequences were filtered based on rounds of selection, such as C3 or C4 rounds.
    NB! Also HTS resulting data was halfed randomly with dplyr::sample_frac to fit into limits. See the “If you want to recalculate the results” section.
  • Positive and negative sequence sets were created for each transcription factor (TF).
  • For artificial data, positive and negative sequences were trimmed, taking 20 nucleotides from each side of the midpoint of the sequence.
    For artificial data, CHS/GHTS sequences were excluded from negative control sequences.
    **Tested, but not included:
    **Using all methods data simultaneously, running on artificial data without sequence trimming, and testing with different HTS subsets.
    (required significantly more processing time and offered minimal or even negative performance gains)
• Used control and experimental sequences to measure motif enrichment, identifying statistically significant motifs for transcription factors.
  • For each TF, the control (negative) and experimental (positive) sequences were used to measure motif enrichment.
  • The pre-generated filtered PWMs from the MEME tool were input for AME to test the enrichment of these motifs and get the respective metrcs.
  • After running AME, the results were processed and saved as .RDS files, containing motifs, their metrics, and associated scores.
    The score was calculated as the mean between the True Positive Rate (TPR) and the True Negative Rate (TNR) per motif.

The pre-calculated results of AME processing are stored in the GitHub repository under the folder:
meme_comparison/data/ame_*.
If you want to recalculate the results for a specific transcription factor, delete all mentions of that TF from the corresponding folder.
The same motifs will highly likely emerge in the top (tested), but they will have updated scores based on your new reference sequences.

• Selected “best” PWMs.
  • For each data type, the pre-calculated AME results were loaded from the corresponding directories (meme_comparison/data/ame_{type}/ for genome data and meme_comparison/data/ame_sep_cut_{type}/ for artificial data).
  • Results from different paths (data_and data_cut) were combined.
  • The top 12 motifs were selected for each transcription factor based on the descending order of score, ensuring distinct motifs, for FIMO Predictions step.
    *only 4, 8, 12 were tested, 12 just have given the best results on lb. No specific function was used for N selection.
  • The top 4 motifs PWMs per transcription factor were exported to text files for submission.

FIMO Predictions (AAA track):

• Leveraged FIMO to scan test sequences using the best PWMs derived from AME. Predicted binding sites were selected based on motif scores and statistical significance.
  • The top 12 motifs derived from the AME were used to scan test sequences.
    *only 4, 8, 12 were tested, 12 just have given the best results on lb. No specific function was used for N selection.
  • For each transcription factor and motif, FIMO was executed to scan the test sequences.
  • Binding site predictions were filtered based on p-values (p-value < 0.001), ensuring only significant motif occurrences were kept.
  • The results (motif occurrences and their scores) were saved for each transcription factor.
  • For each experiment, the FIMO prediction results (stored in .RDS files) were loaded and combined for all transcription factors.
  • The predictions were aggregated by sequence and transcription factor, computing the maximum score for each sequence.
  • The scores were scaled.
    
    Tested, but not included:
    Various aggregation functions were tested instead of using the maximum value, including mean and median calculations. Additionally, several metrics for the score and different ways of normalizing/scaling the score were explored.
    Despite these variations, all approaches have given roughly similar results in terms of performance. Ultimately, the best-performing method was chosen based on results from the leaderboard.

Results & Discussion

Despite experimenting with several machine learning approaches available within my resource constraints (including an xgboost/rf/h2o ensembles on k-mers of different lengths and a CNN on genomic data), this method proved to be the most effective in terms of lb scores and resources. I suspect that, in a less resource-limited environment, e.g. CNN would likely outperform this approach for the AAA track for the artificial data. However, I find the metrics I achieved to be quite acceptable.