S3 : Solutions

Exercise 3.1 - RNA-Seq Mapping and Count Data 😺

Copy the text below and save it as a bash script rnaseq.sh

#!/bin/bash -e

# Jane Doe
# 16 April 2026
#This script uses relative paths. It has to be executed from rna_seq parent directory OR edit paths accordingly. 

# Load required modules and setup the environment
module purge
module load HISAT2/2.2.1-gompi-2023a
module load SAMtools/1.22-GCC-12.3.0
module load Subread/2.0.7-GCC-12.3.0

#Print the current working directory
echo "${PWD}"

#Create results directories. We are moving away from "Mapping" and "Counts"
mkdir -p results/{sam,bam,counts}

# Index genome
genome=ref_genome/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa
index=ref_genome/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel

hisat2-build -p 2 -f $genome $index

# Align to indexed genome
for file in trimmed_reads/*.fastq
do
      # Extract base name
      base=$(basename ${file} .fastq)

      # Set variables
      sam=results/sam/${base}.sam
      summary=results/sam/${base}.summary.txt
      sortedbam=results/bam/${base}_sorted.bam
      mapstat=results/bam/${base}_mapstat.txt

      # Run commands
      hisat2 -p 2 -x $index -U $file -S $sam --summary-file $summary
      samtools view -S -b $sam | samtools sort -o $sortedbam
      samtools flagstat $sortedbam > $mapstat
done

# count how many reads aligned to each genome feature (exon)  
featureCounts -a ref_genome/Saccharomyces_cerevisiae.R64-1-1.99.gtf \
              -o results/counts/yeast_counts.txt -T 2 -t exon -g gene_id \
                results/bam/*_sorted.bam

What is the \ in the featureCounts command line?

You can use the backslash to continue your code over multiple lines. This increases the ease of visibility across very long lines of code. You'll frequently see scripts with each flag on a new line to make it easy to see and edit flag arguments.

Exercise 5.4 - Variant calling workflow slurm script 📜

Copy the text below and save it as a slurm script variant-calling.sl

#!/bin/bash -e

#SBATCH --account       nesi02659
#SBATCH --job-name      variant_calling_workflow
#SBATCH --cpus-per-task 4
#SBATCH --time          00:15:00
#SBATCH --mem           4G
#SBATCH --output        logs/%x-%j.out
#SBATCH --error         logs/%x-%j.err
#SBATCH --mail-type     END
#SBATCH --mail-user     myemail@email.co.nz # remember to update with your own email!

# Optional: create and cd into new working dir   
mkdir ex_5.4 && cd ex_5.4

# Optional: echo back some variables for help debugging later
echo "The date is $(date)"
echo "My current working directory path is $PWD"

# Load all the required modules
module purge
module load BWA/0.7.18-GCC-12.3.0
module load SAMtools/1.22-GCC-12.3.0
module load BCFtools/1.22-GCC-12.3.0


# create the results directories
mkdir -p results/{sam,bam,bcf,vcf}

# indexing the genome
genome=~/scripting_workshop/variant_calling/ref_genome/ecoli_rel606.fasta 
trimmed=~/scripting_workshop/variant_calling/trimmed_reads 

# Optional: re-index here, or re-use the index files we already created
# bwa index $genome

# create a loop that map reads to the genome, sort the bam files and call variants
for fq1 in ${trimmed}/*_1.trim.sub.fastq
 do
    echo "working with file $fq1"

    base=$(basename $fq1 _1.trim.sub.fastq)
    echo "base name is $base"

    # setting the variables
    fq1=${trimmed}/${base}_1.trim.sub.fastq
    fq2=${trimmed}/${base}_2.trim.sub.fastq
    sam=results/sam/${base}.aligned.sam
    bam=results/bam/${base}.aligned.bam
    sorted_bam=results/bam/${base}.aligned.sorted.bam
    raw_bcf=results/bcf/${base}_raw.bcf
    variants=results/vcf/${base}_variants.vcf
    final_variants=results/vcf/${base}_final_variants.vcf

    # running the analysis steps
    bwa mem $genome $fq1 $fq2 > $sam
    samtools view -S -b $sam > $bam
    samtools sort -o $sorted_bam $bam
    samtools index $sorted_bam
    bcftools mpileup -O b -o $raw_bcf -f $genome $sorted_bam
    bcftools call --ploidy 1 -m -v -o $variants $raw_bcf
    vcfutils.pl varFilter $variants > $final_variants

done

echo "ALL DONE!"

Exercise 5.5 - RNA-seq workflow slurm script 📜

Copy the text below and save it as a slurm script rnaseq.sl

#!/bin/bash -e

#SBATCH --account       nesi02659
#SBATCH --job-name      rna-seq_workflow
#SBATCH --cpus-per-task 4
#SBATCH --time          00:15:00
#SBATCH --mem           4G
#SBATCH --output        logs/%x-%j.out
#SBATCH --error         logs/%x-%j.err
#SBATCH --mail-type     END
#SBATCH --mail-user     myemail@email.org.nz

# Optional: create and cd into new working dir   
mkdir ex_5.5 && cd ex_5.5

# Optional: echo back some variables for help debugging later
echo "The date is $(date)"
echo "My current working directory path is $PWD"

# load modules
module purge
module load HISAT2/2.2.1-gompi-2023a
module load SAMtools/1.22-GCC-12.3.0
module load Subread/2.0.7-GCC-12.3.0


# Create results directories. 
mkdir -p results/{sam,bam,counts}

# Set variables.
genomedir=~/scripting_workshop/rna_seq/ref_genome
genome=${genomedir}/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa 
index=${genomedir}/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel
trimmeddir=~/scripting_workshop/rna_seq/trimmed_reads  

# Optional: re-index here, or re-use the index files we already created. 
# hisat2-build -p 4 -f $genome $index

# map reads, convert to BAM and sort
for filename in ${trimmeddir}/*
  do
    base=$(basename ${filename} .fastq)

    # Set variables
    sam=results/sam/${base}.sam
    summary=results/sam/${base}.summary.txt
    sortedbam=results/bam/${base}_sorted.bam
    mapstat=results/bam/${base}_mapstat.txt

    # Run commands
    hisat2 -p 4 -x $index -U $filename -S $sam --summary-file $summary
    samtools view -S -b $sam | samtools sort -o $sortedbam
    samtools flagstat $sortedbam > $mapstat

  done

# count reads per feature
featureCounts \
  -a ${genomedir}/Saccharomyces_cerevisiae.R64-1-1.99.gtf \
  -o results/counts/yeast_counts.txt \
  -T 4 \
  -t exon \
  -g gene_id \
  results/bam/*sorted.bam

echo "All done!"

Back to homepage