Analysis pipeline for RNA-seq data

This pipeline is specifically designed for the analysis of RNA-seq data to fit our B-ALL subtypes prediction pipeline MD-ALL. Upon execution, the pipeline will produce outputs detailing gene read counts, mutations, fusions, and chromosomal-level copy number variations (gains/losses) derived from the RNA-seq data. You can choose to install all the packages yourself or use the Singularity containers we have created, which include all the essential software. The use of Singularity containers is highly recommended, as some of the software can be challenging to install correctly. This pipeline is still under active development, and new analyses regarding B-ALL subtype prediction will be released.

The pipe written in Snakemake. The workflow of this pipeline:

For singularity container users

Please download the singularity containers and the reference database.

Configuration

Users will need to edit the config.yaml file for configurations. Please use the full path for all the reference data as required by Singularity. After editing, please add the pull path of config.yaml in the rnaseq.smk file.

Parameters

‘ref_fa’, the fasta file of reference genome of human GRCh38.

‘gtf’, gtf annotation file of the reference genome.

‘bed_DUX4’, bed file of DUX4 genes. This file is used in the read counts patching process for DUX4 genes.

‘ref_star’, the directory of reference used by STAR to do alignment.

‘ref_fusioncatcher’, the directory of reference used by FusionCatcher to call gene fusions.

‘ref_RNApeg_flat’, the refFlat file used by RNApeg.

‘cores_star’, ‘cores_samtoolsSort’, ‘cores_fusioncatcher’ and ‘cores_RNApeg’ are the number of threads used by the the corresponding software.

Pipeline running

Before running, making sure the Singularity is correctly installed and loaded. Note that all working directories including reference, pipeline and project must be added to singularity using –bind. Otherwise singularity will not be able to use them. The files in current folder does not need to bind. Files Code for running the pipeline:

singularity exec --bind /full_path_to/ref,/full_path_to/pipeline,/full_path_to/project /full_path_to/app_gulab_rnaseq-20231121.sif \
snakemake \
--config dir_in=/full_path_to/project/fq \
dir_out=/full_path_to/project/out \
sample=testId \
-s /full_path_to/pipeline/rnaseq.smk --rerun-incomplete -p --keep-going -j16

The code will initate the pipeline to runing the analysis for sample id testId with 16 cores.

‘dir_in’, the directory of input fastq files. Only gz compressed paired-end fastq files are supported currently. The file names should follow the pattern {sample}.R1.fq.gz and {sample}.R2.fq.gz. If a sample id is COH000456_D1, then the fastq file names should be COH000456_D1.R1.fq.gz and COH000456_D1.R2.fq.gz.

‘dir_out’, the output directory. Results will be stored in sub-directories within this folder, each named according to the respective sample ID.

‘sample’, the sample ID list that will be processed for analysis. The corresponding fastq files need to be stored in the directory ‘dir_in’.

output

‘bam’, sorted and duplication-marked BAM files generated by STAR, Samtools, and GATK MarkDuplicates, along with duplication metrics and STAR junctions.

‘FusionCatcher’, output folder of FusionCather. The file final-list_candidate-fusion-genes.txt will be used for MD-ALL subtype prediction.

‘HTSeq’, HTSeq output folder containing the read counts for each gene, along with the DUX4 gene-patched gene count file.

‘Mutation’, variants called by GATK HapolotyperCaller.

‘RNApeg’, output folder of RNApeg with counts of junctions and can be used as input for Cicero.

‘RNAseqCNV’, Output folder of RNAseqCNV containing the chromosomal-level CNV changes.

‘log’, log files

Cicero

Since fusion calling based on Cicero sometimes runs very slowly, users may choose whether to use it or not. Therefore, we have developed a separate Singularity container for Cicero to facilitate fusion calling. Please download the singularity Cicero containers and the reference database for Cicero. Code for running the Cicero:

singularity exec --bind /full_path_to_Cicero_ref/ref \
/full_path_to/cicero_0.3.0p2.sif \
Cicero.sh -n 8 -b input.bam -g GRCh38_no_alt -r ref -j input.junctions -s 2 -c 10 -o output

It is recommended to use the full path for all input files. The input.junctions file will be the .bam.junctions.tab.shifted.tab file generated by RNApeg.

For non-singularity users

Dependencies

Please install all the necessary software before running this pipeline. The software versions may be interchangeable, but users should proceed with caution.

perl-v5.26.2

Snakemake-7.16.0

STAR-2.7.2b

samtools-1.12

GATK-4.3.0.0

HTSeq-count-0.13.5

FusionCatcher-1.33

RNApeg

Cicero-0.3.0p2

RNAseqCNV-1.2.2

Configuration

Configuration please refer to the singularity part.

Run the pipeline

After all the necessary software is installed, the following code can be used to run the pipeline:

snakemake \
--config dir_in=/full_path_to/fq \
dir_out=/full_path_to/out \
sample=testId \
-s /full_path_to/run_rnaseq.smk --rerun-incomplete -p --keep-going -j16

Contact:

Zunsong Hu: zuhu@coh.org

Zhaohui Gu: zgu@coh.org