Software variables

Assuming that you are in the JASPAR_2022_motif_discovery_and_curation_pipeline folder:

• bin : bin
• python : python3.6
• RSAT : /lsc/rsat

You can see where RSAT program are using the following command: echo $RSAT

Organism specific variables

• data_folder : This is the folder where the input file provided as narrowPeak must be placed. See Expected folder structure for more information about the folder and file names.
• genome_fasta: The path to the .fa file with the organism’s genome sequence.
• genome_size: the path to the .chrom.sizes file.
• TF_Experiment_map: The path to the experiment table. The pipeline expects at least three columns in the following order: 1) Experiment ID, 2) Condition/Cell-type, and 3) TF.
• out_dir: The path to results folder where all the analysis will be stored, one folder per experiment. Don’t forget to indicate the genome in the path. Example: ReMap2020/Hs_hg38

sacCer3 (UCSC)

Download genome from UCSC through command-line:

## Download the twoBitToFa tool for format conversion
## It is not necessary to download twoBitToFa every time you download a genome
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa 
chmod +x twoBitToFa

## Download chromosome sizes file
wget http://hgdownload.cse.ucsc.edu/goldenPath/sacCer3/bigZips/sacCer3.chrom.sizes

## Download the 2bit version of sacCer3 genome
wget http://hgdownload.cse.ucsc.edu/goldenPath/sacCer3/bigZips/sacCer3.2bit

## Get the fasta file
./twoBitToFa sacCer3.2bit sacCer3.fa

Once the .fa and .chrom.sizes were downloaded, set the genome_fasta and genome_size variables in the corresponding config.yaml file

Expected folder structure

The Snakefile will launch a motif analysis + centrality enrichment test for each discovered motif of each dataset (experiment). To launch the pipeline we only required a narrowPeak file that must be located in the ‘data_folder’ specified in the config file (see configuration file section).

Example:

cd data_folder

.
├── ENCSR000BST_GATA3_MCF-7
│   └── ENCSR000BST_GATA3_MCF-7_peaks.narrowPeak
└── ENCSR000ATT_FOXA1_MCF-7
  └── ENCSR000ATT_FOXA1_MCF-7_peaks.narrowPeak

Every EXPERIMENT_FOLDER in the data folder must contain a narrowPeak file with the same name and the suffix _peaks, see tree above. Example:

EXPERIMENT_FOLDER = ENCSR000BST_GATA3_MCF-7
narrowPeak        = ENCSR000BST_GATA3_MCF-7_peaks.narrowPeak

This is the data structure required by the pipeline.

Launching the snakemake workflow

There are two main task/steps in this pipeline:

• Step 1: motif discovery. The program RSAT peak-motifs (8) is called, this program runs four motif discovery algorithms (over-represented k-mers, over-represented spaced-k-mers, positionally biased k-mers, k-mers enriched in central windows) to discover significant motifs on each set of ChIP-seq peaks. This can be launched by running the Snakefile, you can parallelize the process using the argument –cores.

The Snakefile contain a series of rules, for the motif discovery steps, it searches for datasets in narrowPeak format.

• Step 2: selecting centrally enriched motifs. Note that given the ChIP-seq dataset quality and the number of peaks may vary among the datasets, many datasets may not produce any significant motif. A priori, the number of discovered motifs is unknown, so we use checkpoints within the Snakefile.

Selecting config file

This pipeline is generic for any narrowPeak file. In order to use a single pipeline for data coming in multiple input format, we use multiple configuration files (folder config_files) containing the configuration files for multiple organisms and multiple ChIP-seq repositories (e.g., ReMap (9), GTRD (10), ChIP-atlas (11), CistromeDB (12)).

In order to run the pipeline, users must enter one of the following analysis IDs:

• ReMap 2020 Homo sapiens hg38 : ReMap2020_Human
• ChIP-exo S cerevisiae peaks from ChExMix : ChExMix_sacCer3

snakemake --config analysis_id=ReMap2022_Human [--cores X]

All the snakemake rules in the pipeline, except the first one (extract_peak_summits) are the same. The reason is because the input genomic coordinates could vary across the source and may require to parse the table.

In the previous chunk, the argument –config analysis_id=ReMap2022_Human indicates to the pipeline to use the config file for ReMap 2020 Homo sapiens hg38, and run its respective extract_peak_summits rule.

rule RSAT_peakmotifs_per_exp

• RSAT peak-motifs none motif found: verify that the input NarrowPeak is not empty.

• RSAT local-word-analysis python version: this program is written in python and called within RSAT peak-motifs. In case that the default python environment is 3.5, local-word-analysis will not run, therefore the solution is modify directly the line in the $RSAT/perl-scripts/peak-motifs script adding “python2.7” before the program is called.

• RSAT position-analysis founds highly repetitive/uninformative motifs: The default interval for position-analysis is 50nt, which is too large for the JASPAR pipeline (where the peaks have a size of 101), this may produce non-relevant artifact motifs. By changing the interval to 25 the program may found the expected TF motif. Note that this change was already set in the config file.

rule annotate_best_centrimo_experiment

• Verify that the variables in the awk command correspond to the correct fields in your experiment table.

rule best_centrimo_experiment_logo

• Input files not found: this may occur when the Experiment_TF ID contains certain characters as: (, ), [, ], ., /, " "(one or more spaces). Before launching the pipeline, verify that none of the folders in your data folder contains such characters, the simplest option is remove them or substitute them by ‘-’.