Description

This pipeline aims to provide a simple alignment tool for RNAseq data (without UMI). The fastq files are aligned on a reference genome with STAR and counted with HTSeq-Count.

Prerequisites

• Git
• Conda

Input

The fastq.gz files need to be gathered in a directory. The pathway to this directory will be specified in the config.json. One fastq file per sample (or two if data are sequenced in paired end).

Installation

git clone https://github.com/DimitriMeistermann/rnaseq_align.git
cd rnaseq_align
conda env create -f virtualEnvs/rnaseq_align.yml

Configuration

• IS_PAIRED_END: true or false, is the data in paired end (forward and reverse reads ) ?
• PRELOAD_GENOME: true or false (experimental). STAR can preload the genome in the RAM, so the next job can use it. Currently it does not work on bird2cluster.
• STAR_ALIGN_MULTIPLE_FILE: true or false, align fastqs by batches, then resulting BAM is split. This is an alternative way to minimize the number of loading of reference genome if PRELOAD_GENOME is not possible.
• MAX_SIZE_PER_MULTIPLE_ALIGN: Size of alignment batches (if STAR_ALIGN_MULTIPLE_FILE=true). Example for 50 (means 50 gigabytes): if we have 10 samples of 10GB, the alignment will be performed in two batches, from two different STAR_ALIGN jobs. Alignment is also split every 1000 samples, even if MAX_SIZE_PER_MULTIPLE_ALIGN has not been reached.
• READ_LENGTH: String. Number of nucleotid in each reads, automatically picked if empty string by opening the first read of the first fastq. Warning: if cutadapt was already used on fastqs, you have to enter manually the read length before its use.
• STAR_GENOME_THREADS: Number of cores used by genome indexing, preload and align if STAR_ALIGN_MULTIPLE_FILE.
• THREAD_PER_SAMPLE: Number of cores used by most of the workflow jobs.
• FASTA_REFERENCE: Path of fasta file of the reference genome
• GTF_REFERENCE: Path of the GTF file (feature of the reference genome)
• FASTQ_PATH: Directory that contains the input fastqs.
• OUTPUT_PATH: Path where the results will be written.
• PAIR_END_FILE_PATTERN: characters between sample name and 1 or 2 for paired end data. Example: if forward reads are in sample_1.fastq.gz and reverse in sample_2.fastq.gz, then PAIR_END_FILE_PATTERN = "_"
• FORWARD_ADAPTATOR: Forward sequencing adaptator for cutadapt (adaptator trimming). If this argument is set to an empty string, cutadapt will not be used.
• REVERSE_ADAPTATOR: Reverse sequencing adaptator for cutadapt. This argument is not used for single end data.
• GENOME_INDEX_RAM: RAM allowed to be used by genome indexing.
• CUTADAPT_MIN_READ_LEN: Minimum read length after trimming. If the read is to short it will be removed from the fastq.
• DEDUP_UMI: true or false, deduplicate reads based on their UMI. UMI must be stored in read names. The UMI must be at the end of the read name, separated by an underscore, e.g. @[readname]_[UMI]. Please see the documentation of umi_tools extract to see more details on how to add UMI to read names if the sequencing method uses UMI.
• FEATURE_TYPE: Feature of the GTF (second column) were aligned reads are counted (example: Exon, Gene,...).
• FEATURE_ID: ID that will be the rows of the count table (gene_name for gene symbols).

Execution of the pipeline

Normal launch (with 8 cores), with default config.json

conda activate rnaseq_align
snakemake -rp -j 8

With specific config file

conda activate rnaseq_align
snakemake -rp -j 8 --configfile configFileSave/[yourConfig.json]

On SGE grid engine (example: bird2cluster).

You can find this code in clusterFiles/exeSnakemakeBird.sh. You can modify the parameter of jobs by editing jobscriptBird.sh.

conda activate rnaseq_align
snakemake -rp -j 20 --cluster 'qsub'  --jobscript clusterFiles/jobscriptBird.sh --latency-wait 10 --max-jobs-per-second 1 --configfile configFileSave/configTest.json

Test files can be found at https://gitlab.univ-nantes.fr/E114424Z/rnaseq_align/-/tree/master/testFolder/fastq

Output data

• log: log files from executed jobs, in the form of RULE_NAME[_sampleName].(err|out). Additional logs are maybe saved depending the tool that was used.
• FASTQ_CUTADAPT: optional, if cutadapt is used, trimmed fastqs are stored in this folder.
• fastQC: Zip and HTML results from the execution of FASTQC.
• BAM: this directory contains the results of alignment in the form of BAM files.
• DEDUP_BAM: if applicable, this directory contains the BAM deduplicated (based on their UMI)
• counts: this directory contains the resulting counts files obtained from htseqcount.
• results: this directory contains the raw counts table, a table of alignment stats, and the multiqc reports that present different quality scores from the data.