Conda environments

Install the miniconda environment manager :

## https://www.anaconda.com/docs/getting-started/miniconda/install

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

bash ~/Miniconda3-latest-Linux-x86_64.sh

source ~/.bashrc

conda update -n base -c conda-forge conda

conda install mamba -c conda-forge # faster packages manager

Install the workflow manager snakemake in the base environment

mamba install snakemake openpyxl # worflow manager 

Create the environment in your favorite path (-p) from the environment.yaml file:

cd GNT_GuideSeq/

mamba env create -f 01-envs/environment.yml

You should now have an environment named ‘guideseq’ containing all the required programs when running:

#list environments
mamba env list 

# programs version in guideseq environment
mamba list -n guideseq

Reference genomes

The pipeline uses the bowtie2 program to align reads on the reference genome (@langmead2018).

For efficient genome alignment, you can use a pre-built index for Bowtie2. These indices can be downloaded from the Bowtie2 website. Using a pre-built index saves computational resources and time, as it eliminates the need to build the index from scratch.

conda activate guideseq # load the 'guideseq' environment to access tools

# keep only chromosomes, dicard scaffolds
seqkit grep  -r -p '^[chr]?[0-9XYM]+' Homo_sapiens.GRCh38.dna.primary_assembly.fa > your_clean_fasta

bowtie2-build -@ {threads} {your_clean_fasta} {index_prefix}

conda activate           # return to "base" environment

Annotation

Prediction tool

The pipeline uses the SWOffinder program (@yaish2024) to predict gRNA off-targets on the fly and annotate OT accordingly as predicted or not.

• Install from the github repository : https://github.com/OrensteinLab/SWOffinder
• Set the path to SWOffinder in the config file (Prediction)

Prepare samples data sheet (SDS)

The sample data-sheet (SDS) is a simple file provided as delited ( ; ) or xlsx format that contains information about each sample to process in the run. An example is proposed in ./test_dataset/sampleInfo.csv and``./test_dataset/sampleInfo.xlsx.

• sampleName : Sample name to use.
  • Will be use to name output files and folders.
  • Samples with the same sampleName will be merged before processing (Merging samples)
  • Samples name must not include ‘-’ in their name as this symbol is used in the pipeline for a special purpose. Instead, use “_” or any other separator of your choice.
• Genome: Uses to define which reference genome to use for each sample.
  • This must be one of the values present in config file genome key (see Reference genome).
• gRNA_name : name of the gRNA used.
• gRNA_sequence: Sequence of the gRNA without the PAM sequence.
• orientation : which PCR orientation was chosen [“positive”, “negative”, “mix” if they share the same indexes]
  • (this info is only used for metadata purposes, both PCR orientations are automatically processes by the pipeline).
• Cas: name of the Cas used
• PAM_sequence: Sequence of the PAM
• PAM_side: [5 or 3] indicating if the PAM is 5’ or 3’
• Cut_Offset: Distance from gRNA end where the cut occurs [default : -4]
• type: Type of experiment [“guideseq”, “iguideseq”] or other. This value will define which sequence to trim (Reads adapter & ODN trimming sequences). Type of library in the SDS file must match one identical type in the configuration file.
• index1: Sequence of index 1 for demultiplexing
• index2: Sequence of index 2 for demultiplexing

Additional columns can be added for metadata annotation purpose.

SDS format is automatically validated using the snakemake.utils - validate function.

Prepare the configuration file

The configuration file is a yaml formatted file with key-value dictionary used to fine-tune the pipeline behavior. Settings will apply to all samples of the run.

An example of configuration file is proposed in ./test_dataset/guideSeq_GNT.yml

config file must be present in the working directory when starting the pipeline.

author: "Guillaume CORRE, PhD"
affiliation: "Therapeutic Gene Editing - GENETHON, INSERM U951, Paris-Saclay University, Evry, France"
contact: "gcorre@genethon.fr"
version: 'V1.0'

clean_intermediates_files: 'TRUE'

## Library informations
sampleInfo_path: "sampleInfo.csv"
read_path: "."  # path to reads if not current folder
R1: "Undetermined_S0_L001_R1_001.fastq.gz"
R2: "Undetermined_S0_L001_R2_001.fastq.gz"
I1: "Undetermined_S0_L001_I1_001.fastq.gz"
I2: "Undetermined_S0_L001_I2_001.fastq.gz"

## path to references
genome:
  GRCh37:
    fasta: "/PATH_TO_REFERENCE/GRCh37/Sequences/GRCh37.primary_assembly.genome.fa"
    index: "/PATH_TO_REFERENCE/GRCh37/Indexes/Bowtie2/GRCh37.primary_assembly.genome" 
    annotation: "/PATH_TO_REFERENCE/GRCh37/Annotations/gencode.v19.annotation.gtf.gz"
  GRCh38:
    fasta: "/PATH_TO_REFERENCE/GRCh38/Sequences/GRCh38.primary_assembly.genome.fa"
    index: "/PATH_TO_REFERENCE/GRch38/Indexes/Bowtie2/GRCh38.primary_assembly.genome" 
    annotation: "/PATH_TO_REFERENCE/GRCh38/Annotations/gencode.v46.annotation.gtf.gz"
    oncogene_list: "/PATH_TO_REFERENCE/GRCh38/Annotations/OncoList_OncoKB_GRCh38_2025-07-04.tsv"
  GRCm39:
    fasta: "/PATH_TO_REFERENCE/GRCm39/Sequences/GRCm39.primary_assembly.genome.fa"
    index: "/PATH_TO_REFERENCE/GRCm39/Indexes/Bowtie2/GRCm39.primary_assembly.genome" 
    annotation: "/PATH_TO_REFERENCE/GRCm39/Annotations/gencode.vM36.annotation.gtf.gz"
    oncogene_list: ""

################################################
minLength: 25 ## Minimal read length after trimming, before alignment
################################################

## Alignement 
################################################
aligner: "bowtie2"   ## Aligner to use (bowtie2 only for now)
minFragLength: 100         # Minimal fragment length after alignment
maxFragLength: 1500        # Maximal fragment length after alignment 
rescue_R2: "TRUE" # keep R2 from unaligned pairs or pairs with too short R1 reads and align them as single-end reads.
################################################

################################################
## Off targets calling
tolerate_bulges: "TRUE"           # whether to include gaps in the gRNA alignment (this will change the gap penalty during SW pairwise alignment)
max_edits_crRNA: 6                # filter clusters with less or equal than n edits in the crRNA sequence (edits = substitutions + INDELs)
ISbinWindow: 100                  # insertion sites closer than 'ISbinWindow' will be clustered together
minReadsPerUMI: 4                 # Min number of reads per UMIs (>=)
minUMIPerIS: 4                      # Min number of UMI per IS (>=)
slopSize: 50                      # window size (bp) around IS (both directions) to identify gRNA sequence (ie 50bp = -50bp to +50bp)
min_predicted_distance: 100       # distance between cut site and predicted cut site to consider as predicted
################################################

################################################
# post alignment
minMAPQ: 20                      # Min MAPQ score to keep alignments 
UMI_hamming_distance: 1         # min distance to cluster UMI using network-based deduplication, use [0] to keep raw UMIs
UMI_deduplication: "Adjacency"  # method to correct UMI (cluster or Adjacency)
UMI_pattern: "NNWNNWNN"  
UMI_filter: "TRUE"              # If TRUE, remove UMIs that do no match the expected pattern [FALSE or TRUE]
################################################

################################################
# reporting
max_clusters: 100                 # max number of cluster alignments to report
minUMI_alignments_figure: 1       # filter clusters with more than n UMI in the report alignment figure (set to 0 to keep all clusters -> can be slow)
################################################

# Prediction
################################################
SWoffFinder:
  path: "/opt/SWOffinder" ## Path to SWoffinder on your server (downloaded from https://github.com/OrensteinLab/SWOffinder)
  maxE: 6                 # Max edits allowed (integer).
  maxM: 6                 # Max mismatches allowed without bulges (integer).
  maxMB: 6                # Max mismatches allowed with bulges (integer).
  maxB: 3                 # Max bulges allowed (integer).
  window_size: 100
################################################

################################################
# Sequences for the trimming steps

guideseq:
  positive:
    R1_trailing: "GTTTAATTGAGTTGTCATATGT"
    R2_leading: "ACATATGACAACTCAATTAAAC"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"
  negative:
    R1_trailing: "ATACCGTTATTAACATATGACAACTCAA"
    R2_leading: "TTGAGTTGTCATATGTTAATAACGGTAT"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"



iguideseq:
  positive:
    R2_leading: "ACATATGACAACTCAATTAAACGCGAGC"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"
    R1_trailing: "GCTCGCGTTTAATTGAGTTGTCATATGT"
  negative:
    R1_trailing: "TCGCGTATACCGTTATTAACATATGACAACTCAA"
    R2_leading: "TTGAGTTGTCATATGTTAATAACGGTATACGCGA"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"
    
    
    
olitagseq:
  positive:
    R1_trailing: "GGGGTTTAATTGAGTTGTCATATGTT"
    R2_leading: "AACATATGACAACTCAATTAAACCCC"
    R2_trailing: "TCCGCTCCCTCG"
  negative:
    R1_trailing: "CCCATACCGTTATTAACATATGAC"
    R2_leading: "GTCATATGTTAATAACGGTATGGG"
    R2_trailing: "TCCGCTCCCTCG"

tagseq:
  positive:
    R1_trailing: "TGCGATAACACGCATTTCGCATAA"
    R2_leading: "CTTATGCGAAATGCGTGTTATCGCA"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"
  negative:
    R1_trailing: "ATCTCTGAGCCTTATGCGAAATGC"
    R2_leading: "CGCATTTCGCATAAGGCTCAGAGAT"
    R2_trailing: "AGATCGGAAGAGCGTCGTGT"
    
################################################

Folder structure

In order to start a run:

• create a new directory in the installation folder and cd into.

• Then :

  • add the sample data sheet (Prepare samples data sheet (SDS))

  • add the configuration file (Prepare the configuration file)

  • add the sequencing undeterminded_R1/R2/I1/I2.fastq.gz files (undemultiplexed) or symblic link to the raw data.

Input fastq files should respect the following structure from original paper :

• R1 : contains fragment sequence starting in gDNA

• R2: Starts with ODN sequence followed by gDNA sequence and potential trailing adaptor sequence

• i1 : Contains barcode 1 (usually 8 nucleotides)

• i2 : Contains barcode 2 and UMI sequence (usually 8 + 8 nucleotides)

Your folder architecture should look similar to :

Path/to/guideseq/
├── 00-pipeline/
├── 01-envs/
├── 02-resources/
├── 03-misc/
├── test_dataset/
├── My/folder/
    ├── guideSeq_GNT.yml
    ├── sampleInfo.csv
    ├── undertermined_R1.fastq.gz
    ├── undertermined_R2.fastq.gz
    ├── undertermined_I1.fastq.gz
    └── undertermined_I2.fastq.gz

From inside your analysis folder, run the command below after adjusting for number of CPU (-j) :

snakemake -s ../00-pipeline/guideseq_GNT.smk \
  -j 24 \
  -k \
  --use-conda \ 
  --report-after-run --report runtime_report.html  \
  -n
    
# remove the -n argument to start the pipeline

Each rule takes a maximum of 6 threads (alignment, trimming) to speed up the data processing. Set -j as a multiple of 6 to process 2 (12threads) , 3 (18 threads) … samples in parallel.

## usefull snakemake arguments 
# see https://snakemake.readthedocs.io/en/stable/executing/cli.html#all-options

--quiet
--notemp # keep all intermediate files
--rerun-trigger {code,input,mtime,params,software-env} 
--rerun-incomplete

 --filegraph | sed -n '/digraph/,$p' | dot -Tpdf > dag_files.pdf #DAG representation of workflow
 --rulegraph | sed -n '/digraph/,$p' | dot -Tpdf > dag_rules.pdf 

Report

A general report is generated. It summarizes all main QC and key features obtained from the run using graphical representations and tables.

An example is available in the ./test_dataset/results folder.

Off-targets files

For each sample, an Excel file containing all the OT information is created. This file includes numerous columns, some of which are particularly noteworthy. It lists all detected positions, even those without any gRNA matches. For each position, the file reports the total number of UMIs and reads, providing the same information for both positive and negative PCRs. Insertion sites are annotated to genes and oncogenes if they are defined in the configuration file. For OT sites with gRNA matches, a summary of indel and mismatch positions is also provided.

Here is the data formatted into a markdown table:

Demultiplexing

Tool : cutadapt

Undetermined fastq files are demultiplexed to sampleName fastq files according to barcodes present in the sample data sheet.

• i1 and i2 fastq files are concatenated to a single fastq file (i3)

• R1 and R2 fastq are demultiplexed according to i1+i2 sequence present in i3 fastq files.

ODN trimming:

Tool : cutadapt

The leading ODN sequence is remove from R2 reads according to the method defined in the sample data sheet for each sample.

Reads that do not start with the ODN sequence are discarded.

PCR orientation is automatically detected and added to read name for future processing:

@M02111:194:000000000-LT722:1:1101:15810:6788_positive

UMI extraction

Tool : cutadapt

UMI is extracted from the I3 read generate before based on UMI pattern length defined in the configuration file.

UMI is added to reads name:

@M02111:194:000000000-LT722:1:1101:15810:6788_positive_GCTGTAGG

Adaptor trimming:

Tool : cutadapt

The adaptor trailing sequences are trimmed from R1 and R2 reads respectively if present.

Read filtering:

Tool : cutadapt

After trimming, only read pairs with both mates longer than a minLength defined in the configuration file are selected for alignment.

If the rescue_R2 variable is set to TRUE, then R2 reads that are longer that minLength is a discrarded pair are rescued and processed as single-end reads.

Genome alignment:

Tools: Bowtie2

Reference sequences are specified in the configuration file. Genome index is build if it does not already exist.

Trimmed reads pairs that passed the filtering steps are then aligned on the reference genome specified in the sample data sheet for each sample using the aligner specified in configuration file.

Cutting site calling:

Tool : bedtools and awk

Following alignment on the reference genome, nuclease cutting sites are extracted from the start position of R2 read alignment.

Reads that align at the same cutting site with the same UMI are aggregated together keeping track of total reads count.

UMI correction:

Tool : R script

UMI sequences are corrected for potential sequencing error. UMI with less than n edits (hamming distance defined in the configuration file) are clustered together.

Cut site clustering:

Tool: bedtools

Cut sites than fall in the same window of ISbinWindow defined in the configuration file are clustered together.

gRNA match:

Tool: R script

For each cluster of cutting sites, the gRNA sequence defined in the sample data sheet for each sample is looked up in a window of +/- slopSize bp using the Swith-Watterman algorithm.

Gap tolerance can be accepted to detect bulges in gDNA or gRNA if the tolerate_bulges variable is set.

Annotation of clusters:

Tool : R script

Clusters of cutting sites are annotated using the gtf file specified in the configuration file for each organism.

A first step prepare the annotation file to extract only gene and exons features.

A second step annotate clusters to gene and position relative to those genes (exon/intron). Multiple annotations may be present for each cluster and are all reported.

Off target prediction:

Tool : SWOffinder

For each gRNA sequence, PAM sequence and each organism specified in the sample data sheet, a prediction of OTS is realized using the SWOffinder tool with edits and mismatches tolerance defined in the configuration file.

Reporting:

A report is generated with main tables and graphical representations to better understand pipeline results.