RSEM

Descrição

De acordo com a página do RSEM, RSEM (RNA-Seq by Expectation Maximization) é um programa de quantificação de transcriptomas RNA-Seq desenvolvido em 2009. É amplamente utilizado para estimar abundância de transcritos a partir de dados de RNA-Seq.

Versões Disponíveis

  • rsem/1.2.25 (default)

Carregando o Módulo

# Carregar RSEM
module load rsem/1.2.25

# Variáveis disponíveis
echo $RSEM_PATH      # Caminho para executáveis RSEM
echo $EBSEQ_PATH     # Caminho para executáveis EBSeq
echo $BOWTIE2_PATH   # Caminho para Bowtie2

Preparação do Índice

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

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

Quantificação de Amostra Única

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

Processamento em Lote com 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

# Criar script de execução
cat > run_rsem.sh << 'EOF'
#!/bin/bash

while IFS= read -r sample; do
    echo "Processando amostra: $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 para Múltiplas Amostras

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}

Análise com 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

# Preparar matriz de contagem
cat > prepare_ebseq.R << 'EOF'
# Combinar resultados de múltiplas amostras
sample_files <- list.files(path="results", pattern="*.genes.results",
                           recursive=TRUE, full.names=TRUE)

# Extrair contagens esperadas
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

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

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

# Condições experimentais
conditions <- read.table("conditions.txt", header=FALSE)[,1]

# Normalizar
sizes <- MedianNorm(counts)

# Estimar parâmetros
EBOut <- EBTest(Data=counts, Conditions=conditions, sizeFactors=sizes)

# Genes diferencialmente expressos
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

Extração de Resultados

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

# Criar matriz de contagem para todos os genes
cat > extract_counts.R << 'EOF'
# Combinar todos os resultados
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)

# Estatísticas resumidas
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

Referências

Ver também