RSEM

Description

According to the page of RSEM, RSEM (RNA-Seq by Expectation Maximization) is an RNA-Seq transcriptome quantification program developed in 2009. It is widely used to estimate transcript abundance from RNA-Seq data.

Available Versions

  • rsem/1.2.25 (default)

Loading the Module

# Load RSEM
module load rsem/1.2.25

# Available variables
echo $RSEM_PATH      # Path to RSEM executables
echo $EBSEQ_PATH     # Path to EBSeq executables
echo $BOWTIE2_PATH   # Path to Bowtie2

Index Preparation

prepare_rsem_index.sh
#!/bin/bash
#SBATCH -J rsem_index
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 12:00:00
#SBATCH --mem=32G

export INPUT="genome.fa genes.gtf"
export OUTPUT="reference/"

module load rsem/1.2.25

mkdir -p reference

# Prepare RSEM reference
job-nanny rsem-prepare-reference -p $SLURM_CPUS_PER_TASK \
                 --gtf genes.gtf \
                 --bowtie2 --bowtie2-path $BOWTIE2_PATH \
                 genome.fa \
                 reference/mouse_ref

Single Sample Quantification

submit_rsem_single.sh
#!/bin/bash
#SBATCH -J rsem_single
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 24:00:00
#SBATCH --mem=16G

export INPUT="reads_1.fastq reads_2.fastq"
export OUTPUT="sample_quant/"

module load rsem/1.2.25

mkdir -p sample_quant

job-nanny rsem-calculate-expression -p $SLURM_CPUS_PER_TASK \
                 --paired-end \
                 --bowtie2 --bowtie2-path $BOWTIE2_PATH \
                 --estimate-rspd \
                 --append-names \
                 --output-genome-bam \
                 reads_1.fastq reads_2.fastq \
                 reference/mouse_ref \
                 sample_quant/sample

Batch Processing with Script

submit_rsem_batch.sh
#!/bin/bash
#SBATCH -J rsem_batch
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 72:00:00
#SBATCH --mem=32G

export INPUT="samples.txt"
export OUTPUT="all_samples/"

module load rsem/1.2.25

mkdir -p all_samples

# Create execution script
cat > run_rsem.sh << 'EOF'
#!/bin/bash

while IFS= read -r sample; do
    echo "Processing sample: $sample"

    rsem-calculate-expression -p $SLURM_CPUS_PER_TASK \
        --paired-end \
        --bowtie2 --bowtie2-path $BOWTIE2_PATH \
        --estimate-rspd \
        --append-names \
        ${sample}_1.fastq ${sample}_2.fastq \
        reference/mouse_ref \
        all_samples/${sample}
done < samples.txt
EOF

chmod +x run_rsem.sh
job-nanny ./run_rsem.sh
samples.txt
sample1
sample2
sample3
sample4
sample5

Job Array for Multiple Samples

submit_rsem_array.sh
#!/bin/bash
#SBATCH -J rsem_array
#SBATCH --array=1-10
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 36:00:00
#SBATCH --mem=16G

SAMPLES=(
    "sample1"
    "sample2"
    "sample3"
    "sample4"
    "sample5"
    "sample6"
    "sample7"
    "sample8"
    "sample9"
    "sample10"
)

SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
export INPUT="${SAMPLE}_1.fastq ${SAMPLE}_2.fastq"
export OUTPUT="results/${SAMPLE}/"

module load rsem/1.2.25

mkdir -p results/${SAMPLE}

job-nanny rsem-calculate-expression -p $SLURM_CPUS_PER_TASK \
                 --paired-end \
                 --bowtie2 --bowtie2-path $BOWTIE2_PATH \
                 --estimate-rspd \
                 --append-names \
                 --output-genome-bam \
                 ${SAMPLE}_1.fastq ${SAMPLE}_2.fastq \
                 reference/mouse_ref \
                 results/${SAMPLE}/${SAMPLE}

Analysis with EBSeq

submit_ebseq.sh
#!/bin/bash
#SBATCH -J ebseq
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 12:00:00
#SBATCH --mem=8G

export INPUT="conditions.txt"
export OUTPUT="ebseq_results/"

module load rsem/1.2.25

mkdir -p ebseq_results

# Prepare count matrix
cat > prepare_ebseq.R << 'EOF'
# Combine results from multiple samples
sample_files <- list.files(path="results", pattern="*.genes.results",
                           recursive=TRUE, full.names=TRUE)

# Extract expected counts
counts <- data.frame()
for (f in sample_files) {
    data <- read.table(f, header=TRUE, row.names=1)
    sample_name <- basename(dirname(f))
    counts[[sample_name]] <- data$expected_count
}

write.table(counts, "ebseq_results/counts.txt", sep="\t", quote=FALSE)
EOF

Rscript prepare_ebseq.R

# Run EBSeq
cat > run_ebseq.R << 'EOF'
library(EBSeq)

# Read counts
counts <- read.table("ebseq_results/counts.txt", header=TRUE, row.names=1)

# Experimental conditions
conditions <- read.table("conditions.txt", header=FALSE)[,1]

# Normalize
sizes <- MedianNorm(counts)

# Estimate parameters
EBOut <- EBTest(Data=counts, Conditions=conditions, sizeFactors=sizes)

# Differentially expressed genes
PP <- GetPPMat(EBOut)
DE <- rownames(PP)[PP[,2] > 0.95]

write.table(DE, "ebseq_results/DE_genes.txt",
            row.names=FALSE, col.names=FALSE, quote=FALSE)
write.table(PP, "ebseq_results/PP_matrix.txt", sep="\t", quote=FALSE)
EOF

Rscript run_ebseq.R

Results Extraction

extract_rsem_results.sh
#!/bin/bash
#SBATCH -J extract_rsem
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 04:00:00
#SBATCH --mem=4G

module load rsem/1.2.25

# Create count matrix for all genes
cat > extract_counts.R << 'EOF'
# Combine all results
files <- list.files(path="results", pattern="*.genes.results",
                    recursive=TRUE, full.names=TRUE)

counts <- data.frame()
tpm <- data.frame()

for (f in files) {
    data <- read.table(f, header=TRUE, row.names=1)
    sample_name <- basename(dirname(f))

    if (nrow(counts) == 0) {
        counts <- data.frame(row.names=rownames(data))
        tpm <- data.frame(row.names=rownames(data))
    }

    counts[[sample_name]] <- data$expected_count
    tpm[[sample_name]] <- data$TPM
}

write.table(counts, "gene_counts.txt", sep="\t", quote=FALSE)
write.table(tpm, "gene_tpm.txt", sep="\t", quote=FALSE)

# Summary statistics
summary_stats <- data.frame(
    Sample = colnames(counts),
    Total_Counts = colSums(counts),
    Detected_Genes = colSums(counts > 0)
)

write.table(summary_stats, "summary_stats.txt",
            sep="\t", quote=FALSE, row.names=FALSE)
EOF

Rscript extract_counts.R

References

See also