RNASeq: Obtaining read counts

Read counting implies counting the number of reads that map inside a specific annotation feature. The tutorials listed here demonstrate read counting as part of differential gene expression using the R library DESeq/DESeq2. Alternatively, reads may be counted with the python program HTSeq-count, see the manual for instructions (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html).

Read counting may be CPU-intensive, depending on the size of the BAM file(s) used. It is thus recommended to run this process as a job script on Abel. Such a job script must first load the R module on Abel, subsequently executing an R script containing the read-counting R code. Such a job script may look like:

!/bin/bash

SBATCH --job-name=my_R_script_name

SBATCH --account=myAccountName

SBATCH --time=48:00:00

SBATCH --mem-per-cpu=3500M

SBATCH --nodes=1

SBATCH --ntasks-per-node=1

source /cluster/bin/jobsetup

module purge module load R R CMD BATCH /path/to/Rscript.R