Genomic Analyses

Basic tools for processing and analyzing genomic sequencing (next-generation sequencing, NGS) data are available preinstalled in a container environment.

The path to the container on Komondor is: /opt/software/packages/containers/datascience_genomics.sif. The definition file of the container (datascience_genomics.def), together with the conda environment specifications needed to build it, is available in the “HPC Bioinformatics / Containers” public git repository.

Note

The previous version of the container is still available in the same directory as /opt/software/packages/containers/datascience_genomics_20250121.sif (e.g. for reproducing earlier analyses). Tool versions differ between the two containers.

Quick start

$ module load apptainer
$ SIF=/opt/software/packages/containers/datascience_genomics.sif

# run a tool from the container as a Slurm job
$ srun apptainer exec $SIF samtools --version

# short description of the container (tool groups, special commands)
$ srun apptainer run-help $SIF

# exact tool versions (one file per tool group)
$ srun apptainer exec $SIF ls /opt/genomics/lock/
$ srun apptainer exec $SIF cat /opt/genomics/lock/core.yml

Always run the tools of the container as Slurm jobs, not directly on the login node: short commands with srun as above, longer analyses with sbatch scripts (see the examples below).

Note

The short commands above run with the default resources (1 task, 1 CPU core). For real analyses always set the resources you need with the --ntasks, --cpus-per-task, --mem-per-cpu and --time (and --partition) options, as in the examples below. See the list of the options in Preparing Job Scripts (sbatch).

Only the $HOME directory is mounted in the container automatically: use the -B option to mount the PROJECT and SCRATCH file systems (e.g. apptainer exec -B /project,/scratch $SIF ...).

Available tools

Open-source tools that are currently available in the container:

  • Quality control of sequencing data: FastQC, MultiQC

  • Trimming and cleaning sequences: Cutadapt, Trimmomatic, Trim Galore

  • Sequence read mapping tools: BWA, Bowtie, Bowtie 2, STAR, BBMap, HISAT2, Bismark, minimap2

  • Multiple sequence alignment: MUSCLE, Clustal Omega, ClustalW, MAFFT

  • Manipulation of sequences and alignment data analysis: Samtools, Picard, MACS2, deepTools, HOMER, HTSeq, tabix/bgzip, SeqKit

  • Genome arithmetics, genomic features: Bedtools, GenomeTools

  • Gene expression analysis: Salmon, kallisto, DESeq2

  • Genetic variant identification and analysis: GATK4, BCFtools, VCFtools, freebayes

  • Population genetics and phylogenetics: ANGSD (with NGSadmix), PCAngsd, ngsDist, IQ-TREE 3, TASSEL

  • Structure analysis and prediction: ViennaRNA

  • Chromosome conformation capture analysis: HiCUP, HiCExplorer, HiC-Pro, hic2cool, cooler

  • Genome assembly: ABySS, MEGAHIT, Velvet, SPAdes, Canu, Trinity, MetaPlatanus, Tigmint, TGS-GapCloser, Racon, MetaBAT 2

  • Sequence analysis and identification: BLAST+, HMMER, TRF, RepeatMasker, RepeatModeler, RepeatScout, CD-HIT

  • Bioconductor packages: DESeq2, Gviz, ggbio, trackViewer (and ggpubr), available through the R-bioc / Rscript-bioc commands and the R (Bioconductor) Jupyter kernel

  • Data analysis and scripting: R, Python (Jupyter datascience stack with Biopython and py3Dmol), Julia, Perl, Ruby

See the “HPC Bioinformatics / Containers public git repository” for detailed software list and versions and more info about the container.

If you need a tool or package that is not in the container, you can build your own environment (see Building your own environment below).

Important

Note that authors usually request a citation in the publication if their software tool was used in the published work. Please always check the terms of use and the citation request on the tool’s website.

Scheduled analyses (sbatch)

Using the container, bioinformatics pipeline scripts can be run as scheduled Slurm jobs.

Example: Bowtie2 alignment

In this example, suppose we have paired-end NGS data from an experiment that we want to align to the human reference genome.

We have the reference genome prepared (indexed) in the reference directory (on the PROJECT filesystem). We have the two NGS fastq (read1 and read2) files in the data directory and we created an output directory for the alignment files (on the faster SCRATCH filesystem).

$ ls -1 /project/my_project/reference
hs_grch38.fa.gz
hs_grch38.fa.gz.1.bt2
hs_grch38.fa.gz.2.bt2
hs_grch38.fa.gz.3.bt2
hs_grch38.fa.gz.4.bt2
hs_grch38.fa.gz.rev.1.bt2
hs_grch38.fa.gz.rev.2.bt2

$ ls -1 /scratch/my_project/fastq
my_sample_R1.fastq.gz
my_sample_R2.fastq.gz

$ ls -1 /scratch/my_project/bam
$

Now we prepare an sbatch script to do the alignment using bowtie2 that is already installed in the genomics container (/opt/software/packages/containers/datascience_genomics.sif). We allocate 16 processor cores for the job and set bowtie2 to use 16 threads accordingly. The output of bowtie2 is passed directly to samtools sort, so we get a sorted BAM file instead of a large SAM file; finally, samtools index creates the index of the BAM file. The pipe (|) has to run inside the container, so the commands are given to bash -c in the container.

Note, that we mount both the PROJECT and SCRATCH file systems for the container since we have input data in both locations. As configured in the script, bowtie2’s log will be in the “slurm-%j__bowtie2.out” file (%j is replaced by the Slurm job ID).

$ cat ./alignment.sbatch
#!/usr/bin/env bash

#SBATCH --job-name=alignment
#SBATCH --partition=cpu
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --mem-per-cpu=2000
#SBATCH --time=01:00:00
#SBATCH --output=slurm-%j__bowtie2.out

# input data (fastq files)
READ1=/scratch/my_project/fastq/my_sample_R1.fastq.gz
READ2=/scratch/my_project/fastq/my_sample_R2.fastq.gz
# output data (sorted bam file)
BAM_OUT=/scratch/my_project/bam/my_sample.bam
# reference genome index location
REF_IDX=/project/my_project/reference/hs_grch38.fa.gz

# loading Apptainer
module load apptainer

# Apptainer options
#   the genomic analysis tools are installed in the datascience_genomics.sif container
#   -B option specifies the file systems (with our data) we want to mount inside the container
CONTAINER_CALL="apptainer exec -B /project,/scratch /opt/software/packages/containers/datascience_genomics.sif"

# bowtie2 alignment (on 16 processor cores, as set in the SBATCH options above),
# the output is sorted by samtools (4 threads) into a bam file
srun ${CONTAINER_CALL} bash -c "bowtie2 --threads ${SLURM_CPUS_PER_TASK} -x ${REF_IDX} -1 ${READ1} -2 ${READ2} \
                                | samtools sort -@ 4 -o ${BAM_OUT}"

# index of the bam file
srun ${CONTAINER_CALL} samtools index ${BAM_OUT}

We can submit the sbatch script to Slurm using the following command:

$ sbatch ./alignment.sbatch
Submitted batch job 1234567

If the job completes successfully, the sorted alignment (bam) file and its index have been created and we see a bowtie2 log with the alignment stats:

$ ls -1 /scratch/my_project/bam
my_sample.bam
my_sample.bam.bai

$ cat ./slurm-1234567__bowtie2.out
50000 reads; of these:
  50000 (100.00%) were paired; of these:
    10646 (21.29%) aligned concordantly 0 times
    33390 (66.78%) aligned concordantly exactly 1 time
    5964 (11.93%) aligned concordantly >1 times
    ----
    10646 pairs aligned concordantly 0 times; of these:
      5243 (49.25%) aligned discordantly 1 time
    ----
    5403 pairs aligned 0 times concordantly or discordantly; of these:
      10806 mates make up the pairs; of these:
        6544 (60.56%) aligned 0 times
        1890 (17.49%) aligned exactly 1 time
        2372 (21.95%) aligned >1 times
93.46% overall alignment rate

Example: Multiple alignments in parallel (job array)

In this example, the task is similar to the previous example, but now we have three samples instead of one that we want to analyze in parallel.

We have the reference (indexed), the input files (three paired-end NGS data files) and the output directory:

$ ls -1 /project/my_project/reference
hs_grch38.fa.gz
hs_grch38.fa.gz.1.bt2
hs_grch38.fa.gz.2.bt2
hs_grch38.fa.gz.3.bt2
hs_grch38.fa.gz.4.bt2
hs_grch38.fa.gz.rev.1.bt2
hs_grch38.fa.gz.rev.2.bt2

$ ls -1 /scratch/my_project/fastq
sample1_R1.fastq.gz
sample1_R2.fastq.gz
sample2_R1.fastq.gz
sample2_R2.fastq.gz
sample3_R1.fastq.gz
sample3_R2.fastq.gz

$ ls -1 /scratch/my_project/bam
$

The simplest method to do the alignments in parallel is using a job array. With the following settings, the sbatch script will be run in parallel for all the three samples as separate jobs. The variable $SLURM_ARRAY_TASK_ID can be used to distinguish the runs and specify the sample name individually for each run. The range of the array task IDs is set with the --array sbatch option. In the name of the output file, %A is replaced by the job ID of the array and %a by the array task ID. (See also the “Job Array” section in “Basic Usage / Advanced Job Definitions”.)

$ cat ./alignment_array.sbatch
#!/usr/bin/env bash

#SBATCH --job-name=alignment
#SBATCH --partition=cpu
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --mem-per-cpu=2000
#SBATCH --array=0-2
#SBATCH --time=01:00:00
#SBATCH --output=slurm-%A_%a__bowtie2.out

# loading Apptainer
module load apptainer

# Apptainer options
#   the genomic analysis tools are installed in the datascience_genomics.sif container
#   -B option specifies the file systems (with our data) we want to mount inside the container
CONTAINER_CALL="apptainer exec -B /project,/scratch /opt/software/packages/containers/datascience_genomics.sif"

# we want to run the alignment for three samples in parallel
samples=(sample1 sample2 sample3)

# current run (sample)
#   we use the index of the current task within the job array to select the sample name
run=${samples[$SLURM_ARRAY_TASK_ID]}

# input data (fastq files)
READ1=/scratch/my_project/fastq/${run}_R1.fastq.gz
READ2=/scratch/my_project/fastq/${run}_R2.fastq.gz
# output data (sorted bam file)
BAM_OUT=/scratch/my_project/bam/${run}.bam
# reference genome index location
REF_IDX=/project/my_project/reference/hs_grch38.fa.gz

# bowtie2 alignment (on 16 processor cores, as set in the SBATCH options above),
# the output is sorted by samtools (4 threads) into a bam file
srun ${CONTAINER_CALL} bash -c "bowtie2 --threads ${SLURM_CPUS_PER_TASK} -x ${REF_IDX} -1 ${READ1} -2 ${READ2} \
                                | samtools sort -@ 4 -o ${BAM_OUT}"

# index of the bam file
srun ${CONTAINER_CALL} samtools index ${BAM_OUT}

Submitting the sbatch script to Slurm:

$ sbatch ./alignment_array.sbatch
Submitted batch job 1234567

On a successful run, we will have the output (bam) file for each sample, and a separate output log for each alignment:

$ ls -1 /scratch/my_project/bam
sample1.bam
sample1.bam.bai
sample2.bam
sample2.bam.bai
sample3.bam
sample3.bam.bai

$ ls -1 ./slurm-*__bowtie2.out
slurm-1234567_0__bowtie2.out
slurm-1234567_1__bowtie2.out
slurm-1234567_2__bowtie2.out

Useful configurations

Links to Slurm configurations that can help organize and efficiently run your bioinfo pipelines that consist of many steps and/or parallel tasks:

Job arrays (mentioned above): Submitting and managing collections of similar jobs quickly and easily - https://slurm.schedmd.com/job_array.html

Packed jobs: Mechanism for resource management to the job within its allocation regarding multiple job steps - https://slurm.schedmd.com/srun.html#OPT_exclusive

Job dependencies: Organizing sequential compute jobs that have to be run in order - https://slurm.schedmd.com/sbatch.html#OPT_dependency

Heterogeneous jobs: Submitting and managing jobs with different resources for different pieces of the job - https://slurm.schedmd.com/heterogeneous_jobs.html

Multiple programs configuration: Run a job with different programs and different arguments for each task - https://slurm.schedmd.com/srun.html#OPT_multi-prog

Interactive analyses

The genomics container can also be launched interactively on the command line or on the JupyterHub.

Execute the analysis interactively

Instead of preparing and submitting an sbatch script, we can also run a short test alignment interactively on the command line using the container (provided that the task is really short and does not use too many resources). We use the same setup (reference, data files, task) as described in the “Example: Bowtie2 alignment” section above:

$ REF_IDX=/project/my_project/reference/hs_grch38.fa.gz
$ READ1=/scratch/my_project/fastq/my_sample_R1.fastq.gz
$ READ2=/scratch/my_project/fastq/my_sample_R2.fastq.gz
$ BAM_OUT=/scratch/my_project/bam/my_sample.bam

$ module load apptainer

$ CONTAINER_CALL="apptainer exec -B /project,/scratch /opt/software/packages/containers/datascience_genomics.sif"
$ srun --ntasks=1 --cpus-per-task=16 --mem-per-cpu=2000 ${CONTAINER_CALL} \
      bash -c "bowtie2 --threads 16 -x ${REF_IDX} -1 ${READ1} -2 ${READ2} | samtools sort -@ 4 -o ${BAM_OUT}"

Execute the analysis in the container shell

We can also load the genomics container for full interactive use and issue commands directly in the container’s shell. The environment variables set on the host (here REF_IDX, READ1, READ2 and BAM_OUT from the previous section) are also available inside the container:

$ module load apptainer

$ CONTAINER_SHELL="apptainer shell -B /project,/scratch /opt/software/packages/containers/datascience_genomics.sif"
$ srun --ntasks=1 --cpus-per-task=16 --mem-per-cpu=2000 --pty ${CONTAINER_SHELL}

Apptainer> bowtie2 --version
/opt/genomics/envs/core/bin/bowtie2-align-s version 2.5.5
...
Apptainer> realpath $READ1
/scratch/my_project/fastq/my_sample_R1.fastq.gz

Apptainer> bowtie2 --threads 16 -x ${REF_IDX} -1 ${READ1} -2 ${READ2} | samtools sort -@ 4 -o ${BAM_OUT}
50000 reads; of these:
  50000 (100.00%) were paired; of these:
...

Apptainer> ls /scratch/my_project/bam
my_sample.bam

Apptainer> exit
$

Interactive use on the JupyterHub

Login to the JupyterHub and select datascience-genomics from the container list when starting a new server. You can use the bash kernel to create a notebook for the analysis or use the bash terminal to run the tools directly. For the Bioconductor packages (e.g. DESeq2), use the R (Bioconductor) kernel.

Known limitations

  • Two R installations: the R and Rscript commands start the R of the Jupyter stack, without the Bioconductor packages. For DESeq2, Gviz, ggbio, trackViewer and ggpubr, use the R-bioc / Rscript-bioc commands or the R (Bioconductor) Jupyter kernel.

  • Python libraries of the tools: the tools are installed in separate environments, so their Python libraries (e.g. deepTools, pysam, cooler, MACS2, HiCExplorer) cannot be imported in Python scripts or notebooks; use them as command-line tools. From Python, only the Jupyter datascience stack, Biopython and py3Dmol are available.

  • Read-only container: new packages cannot be installed into the container (conda install or pip install does not work); see Building your own environment.

  • Repeat libraries: the Dfam and RepBase repeat libraries are not included for RepeatMasker / RepeatModeler; you have to provide your own repeat library (e.g. with the -lib option of RepeatMasker).

  • No GPU support: the container contains no CUDA/GPU libraries, run it on the CPU partition.

  • Java tools: start Picard, Trimmomatic and TASSEL with their commands (picard, trimmomatic, run_pipeline.pl), not with java -jar.

  • Host environment: Apptainer passes the environment variables of the host into the container. If you get unexpected Python or R errors, check that no conda environment is activated and PYTHONPATH / R_LIBS_USER are not set when you start the container.

Building your own environment

If you need tools or package versions that are not available in the genomics container, you can build your own environment. See Conda environments with Apptainer for building a container from a conda environment (conda-forge and bioconda packages), and Installing Python packages for the options of installing Python packages.