distal-seq

Bash script for aligning reads, calling peaks and running the necessary scripts for distal-seq. Takes paired-end fastq files and all necessary options as input. Output consists of alignment stats and distal-seq blender output figures and tables.

setup of conda environment

Create the conda environment with the following command:

conda env create -f environment.yml

usage

conda activate distal-seq
./distal-seq_pipeline.sh <path/to/R1.fastq.gz> <path/to/R2.fastq.gz> <NAME> <ARM1> <ARM2> <CUTOFF> <FILTERCUTOFF> <BOWTIETHREADS> <SAMTOOLSTHREADS> <FORCEBLENDER> <FIGUREONLY> <MDOE> <WINDOW> <PATH/TO/BOWTIEREF> <PATH/TO/GENOME.fa>

The arguments are: R1: Path to R1 fastq.gz file R2: Path to R2 fastq.gz file NAME: Unique sample/run name (used as output folder name) ARM1: Arm1 sequence ARM2: Arm2 sequence CUTOFF: Number of blunt end reads used as threshold for raw data (normally 5) FILTERCUTOFF: Disco score cutoff to define high confidence hits in the processed filtered output (normally 15) BOWTIETHREADS: Number of threads for bowtie2 SAMTOOLSTHREADS: Number of threads for each samtools sub-command FORCEBLENDER: 0 or 1. Re-run blender even if output files exist (1) or not (0) FIGUREONLY: 0 or 1. Re-run figure generation if forceblender is 0 and FIGUREONLY is 1 MODE: Mode for calculating arm alignment differences. Using substitution score ("substitution") or total count of mismatches ("mismatches") WINDOW: Search window for arms around cut site (normally 40) BOWTIE_REF: Path to bowtie2 reference ("/path/to/Bowtie2Index/genome") FASTA_REF: Path to genome fasta file ("/path/to/genome.fa")

example command for test data

To run the sub-sampled AAVS1 fastq files from the data folder you need your own hg38 bowtie2 index and fasta file. Then use the following command:

conda activate distal-seq
./distal-seq_pipeline.sh data/AAVS1_sub_R1.fastq.gz data/AAVS1_sub_R2.fastq.gz AAVS1 ACCCCCCTCCTCACCACA GCTGGCTACTGGCCTTAT 5 15 40 10 1 0 substitution 40 <bowtie2 index> <genome fasta>