Bulk RNAseq tutorial

1. Setup

Connect to HIVE, initiate a screen,

Initiate a screen to make sure if you disconnect accidentally from your terminal, the screen window will keep your sesssion going.

ssh {your_username}@hive.hpc.ucdavis.edu
screen

To detach or exit a screen

Detach from screen with Ctrl+A then Ctrl+D, this will keep the screen open and return you to the previous screen you were in. To exit and close a screen type exit. This will close any interactive node or modules you have loaded.

To list and reattach a screen

screen -ls
There are screens on:
    1707831.pts-134.login2  (09/16/2026 12:58:58 PM)  (Attached)
    3353424.pts-8.login2    (09/11/2026 09:40:16 AM)  (Detached)
screen -r 1707831.pts-134.login2

Request an interactive node

Important: Every time you connect to the HIVE cluster, you will initially have a session on the login node (login server). Because this computer is shared among all users and has limited computing resources, it’s bad practice to run even mild computations on this node. The first thing you should do after logging into the cluster is to ask for an interactive session on a worker node, by running the command:

#inside your screen
salloc --mem=40g -c 8 --partition=high --time=03:00:00 --account=luisccgrp

This opens up a 3 hour session with 8 CPUS/threads and up to 40GB of memory. We use 40GB memory mainly for the STAR alignment step as the genome we use is 38GB in size and must be loaded into the memory. You can allocate less memory if not doing STAR alignment. For full-sized samples, I would recommend allocating 16 cpus. Always ask for more time than you’ll actually need. After a few seconds (hopefully!), you should have been given an interactive session on a worker node. You can tell that you have moved between nodes by checking at the server name in your command prompt: the name of hive login node is loginX, while worker nodes have specific names e.g. hive-as-11-2-42.

Connect the HIVE to your laptop computer

sshfs -o allow_other,ro username@hive.hpc.ucdavis.edu:/quobyte/luisccgrp HIVE

In case you need to install sshfs on you LAPTOP WSL2 UBUNTU (not HIVE)

#from you laptop 
sudo apt update
sudo apt install sshfs -y
#you need to change your fuse.conf file to allow_other 
#remove the "#" before "allow_other" using vim or nano 
vim /etc/fuse.conf

You’ll be able to see HIVE files from your Windows folder by going to: \\wsl$\Ubuntu\home\yourusername

Quality control with FastQC

FastQC is a program that parses FASTQ read files and outputs a number of human-readable summary statistics and plots. Let’s use it to check the quality of the sequence data for our samples. The HIVE cluster has FastQC installed, you can load the program:

### list modules that are available on HIVE
module avail
module avail f  #shows only modules starting with the letter "f"

### load the fastqc module
module load fastqc

### get help to use fastqc
fastqc --help | less

Note: --help or -h are standard options, which many programs recognize as a request for information about the software.

So now, we can run FastQC on the reads of SampleF1_1.fq.gz like this:

cd /quobyte/luisccgrp/USERS/RNAseq_tutorial/<yourname>
mkdir fastqc_results  #making folder for fastqc outputs 
fastqc -t 4 -o fastqc_results SampleF1_1.fq.gz  
# Option (-t) sets number of threads/cpus
# Option (-o) sets the output folder otherwise will create in same folder by default

# create a "for loop" for all *fq.gz file in a directory
## First see what the loop actual outputs without executing the commands
cd /quobyte/luisccgrp/USERS/RNAseq_tutorial/<yourname>
for f in ./samples/*fq.gz; do
    echo "fastqc -t 4 -o fastqc_results $f"
    #fastqc -t 4 -o fastqc_results "$f"
done

# Actually output the loop
for f in ./samples/*fq.gz; do
    echo "fastqc -t 4 -o fastqc_results $f"
    fastqc -t 4 -o fastqc_results "$f"
done
### check files in the folder
ls fastqc_results
  • Run FastQC on both sequence files. Wait for the program to complete (this should only take about a minute).
  • How many output files did each run create? (Note: You can list the files in a directory in the order in which they were create using ls -ltr.)
  • How can we view the output files? (Does it require file transfer to local machine?) You can view html with results summary in Windows You can learn how to interpret results at: https://hbctraining.github.io/Intro-to-rnaseq-hpc-salmon/lessons/qc_fastqc_assessment.html or ask ChatGPT or Claude! If you have a severe QC issue, talk with Paul or I to see if sample is salvagable or talk with Paul or I.

Filtering low-quality reads with fastp

The first step in a sequence data analysis is usually to remove the subset of the data that has insufficient quality – keeping unreliable reads and base calls can introduce unnecessary noise in the analysis. This includes the Illumina adapter sequence that can contaminate the beginning of the reads.

TruSeq single index (previously LT) and TruSeq CD index (previously HT)-based kits:

Read 1: AGATCGGAAGAGCACACGTCTGAACTCCAGTCA
Read 2: AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT

To do so, we will use the program fastp (https://github.com/opengene/fastp). The HIVE cluster also provides this as a module ### Installing fastp

module load fastp
fastp --help

Tip–Using the time command

In addition to the fastp command, I’m adding the time command, this will output how long our fastp command take for one sample, which will help us decide how to run all the samples and if they need to be run in batches. Always pay attention to the real time output, this how much actual time passes.

Let’s use the time command in the fastp process we run below for 1 sample.

#make a directory for the new trimmed FASTQ
mkdir trimmed_fastq
time fastp \  #adding the time command will measure time for the command immediately after it
  -i ./samples/SampleF1_1.fq.gz -I ./samples/SampleF1_2.fq.gz \
  -o ./trimmed_fastq/SampleF1_1.trimmed.fastq.gz -O ./trimmed_fastq/SampleF1_2.trimmed.fastq.gz \
  --detect_adapter_for_pe \
  --cut_front \
  --cut_tail \
  --length_required 30 \  ##use 50 if concentrating on splicing 
  --thread 8 \
  --html ./trimmed_fastq/SampleF1.fastp.html --json /trimmed_fastq/SampleF1.json

Here we set the parameters for trimming fastq files for general differential analysis, with splice variant analysis you can be a little more stringent to get longer, better reads. However, this difference is likely neglible. alt text

  • Run fastp on both sequence files. Wait for the program to complete (this should only take about a minute).
  • How many output files did each run create? View the output files

Make a loop bash script and loop over all your files

vim run_fastp.sh
#Enter the full script below and change as necessary for your files
#!/bin/bash
set -euo pipefail

mkdir -p ./trimmed_fastq

for R1 in ./samples/*_1.fq.gz; do
    # Derive R2, sample base name, and output paths
    R2="${R1%_1.fq.gz}_2.fq.gz"
    SAMPLE=$(basename "${R1%_1.fq.gz}")

    echo "=== Processing ${SAMPLE} ==="

    fastp \
        -i "$R1" -I "$R2" \
        -o "./trimmed_fastq/${SAMPLE}_1.trimmed.fastq.gz" \
        -O "./trimmed_fastq/${SAMPLE}_2.trimmed.fastq.gz" \
        --detect_adapter_for_pe \
        --cut_front \
        --cut_tail \
        --length_required 30 \
        --thread 8 \
        --html "./trimmed_fastq/${SAMPLE}.fastp.html" \
        --json "./trimmed_fastq/${SAMPLE}.fastp.json"

    echo "--- Done: ${SAMPLE} ---"
done

echo "All samples processed."
#make your bash script executable
chmod +x ./run_fastp.sh
#run your script and output the log file of how the progress goes
./run_fastp.sh 2>&1 | tee ./trimmed_fastq/fastp_run.log
# 2>&1 will output the progress screens to the fastp_run.log, this is good to record what commands were run and if any errors were reported
#tee allows you to see the progress as well as writing the log file
less ./trimmed_fastq/fastp_run.log

Post-trimming FastQC

To be safe it’s good to look at the quality of your reads post trimming to make sure they are good and see effect of trimming

#Run loop as before but with new samples
##Make the directory for fastqc post trim results
mkdir fastqc_results/post_trim
for f in ./trimmed_fastq/*trimmed.fastq.gz; do
    echo "fastqc -t 4 -o fastqc_results/post_trim $f"
    fastqc -t 4 -o fastqc_results/post_trim "$f"
done
### check files in the folder
ls fastqc_results/fastqc_results/post_trim

Review a file pre and post trimmed to see how trimming affects the fastq files

Aligning reads to a reference genome

Here, we will align our RNA-seq Illumina reads to the reference human genome using the program STAR.

Creating a STAR genome index

The first step in a sequence data analysis is usually to remove the subset of the data that has insufficient quality – keeping unreliable reads and base calls can introduce unnecessary noise in the analysis.

Read alignment programs typically require the construction of a genome sequence database as a preliminary step before they are able to perform read alignments. This is necessary because aligning directly to a genome sequence in FASTA format is not computationally efficient, but only needs to be run once.

YOU DO NOT NEED TO RUN THIS SCRIPT, I’VE ALREADY DONE IT FOR YOU To create a genome database against which to map reads, STAR needs the genome sequence (FASTA file) and the gene/transcript/exon annotations (GTF file). Gencode https://www.gencodegenes.org/human/ is the website which has the standard gtf and fasta files that are used for RNAseq alignment. Here we will choose the comprehensive gene annotation on the reference chromosomes only gtf and the custom hg38 FASTA file based on the GDC file (https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files). This fasta contains viral sequences (including EBV) and I added the H.pylori genome. It also contains sequence decoys that can soak up reads that come from centromeres..etc (usually more important for DNA sequencing). Be sure to mark down which genome and gtf you use for you data.

AGAIN DON’T RUN THE SCRIPTS BELOW

Use the following pre-built index at: /quobyte/luisccgrp/REFERENCE_DATA/STAR/GRCh38.d1.vd1.Hpylori_gencode.v50/

##downloading the gtf file from gencode
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_50/gencode.v50.annotation.gtf.gz
gunzip gencode.v50.annotation.gtf.gz
#you will also need a FASTA file

#Generating a custom index
STAR
--runMode genomeGenerate \
--genomeDir  \
--genomeFastaFiles <reference> \
--sjdbOverhang 100 \
--sjdbGTFfile <gencode.v50.annotation.gtf> \
--runThreadN 16 \

If you generate any new references, it’s best to keep them in the /quobyte/luisccgrp/REFERENCE_DATA/ folder. ### Aligning reads to the reference genome HIVE has STAR installed as a module. Let’s load the STAR module: Be sure to have a interactive node allocation of 16 cpu and 64G RAM. The memory requirement is the most important as the human genome index is large, minimum is 32GB RAM.

module load star

##Running STAR 2-pass mode similar to GDC RNAseq pipeline DR32 https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/#rna-seq-alignment-command-line-parameters-dr32

These options set thresholds similar to TCGA processed RNAseq data and create outputs ready for differential analysis and STAR-Fusion calling

mkdir ./star_output

time STAR \
--readFilesIn ./trimmed_fastq/SampleF1_1.trimmed.fastq.gz ./trimmed_fastq/SampleF1_2.trimmed.fastq.gz \
--outFileNamePrefix ./star_output/SampleF1. \
--genomeDir /quobyte/luisccgrp/REFERENCE_DATA/STAR/GRCh38.d1.vd1.Hpylori_gencode.v50/ \
--readFilesCommand zcat \
--runThreadN 16 \
--twopassMode Basic \
--outFilterMultimapNmax 20 \
--alignSJoverhangMin 8 \
--alignSJDBoverhangMin 1 \
--outFilterMismatchNmax 999 \
--outFilterMismatchNoverLmax 0.1 \
--alignIntronMin 20 \
--alignIntronMax 1000000 \
--alignMatesGapMax 1000000 \
--outFilterType BySJout \
--outFilterScoreMinOverLread 0.33 \
--outFilterMatchNminOverLread 0.33 \
--limitSjdbInsertNsj 1200000 \
--outSAMstrandField intronMotif \
--outFilterIntronMotifs None \
--alignSoftClipAtReferenceEnds Yes \
--quantMode TranscriptomeSAM GeneCounts \
--outSAMtype BAM Unsorted \
--outSAMunmapped Within \
--genomeLoad NoSharedMemory \
--chimSegmentMin 15 \
--chimJunctionOverhangMin 15 \
--chimOutType Junctions SeparateSAMold WithinBAM SoftClip \
--chimOutJunctionFormat 1 \
--chimMainSegmentMultNmax 1 \
--outSAMattributes NH HI AS nM NM ch

Advanced application: If planning to do precise variant calling based on the RNAseq data, then it’s advised to add the read group information with the following options: –outSAMattrRGline.

After the program completes, check for errors in the terminal and in the STAR log file (especially toward the end), and for existing non-empty output files.

SampleF1.Aligned.out.bam

The primary genome-coordinate BAM file — your main alignment output. Contains all reads with genomic coordinates, including both uniquely and multi-mapped reads, spliced alignments (shown as N operations in CIGAR), and unmapped reads (retained because of –outSAMunmapped Within). This is the file you sort and index for IGV viewing, Picard QC, variant calling, and as input to STAR-Fusion.

SampleF1.Aligned.toTranscriptome.out.bam

Alignments in transcript coordinates rather than genomic coordinates. Each read is mapped to a specific transcript (ENST), with positions relative to the transcript start. Produced by –quantMode TranscriptomeSAM. This is the input for isoform-level quantification tools:

RSEM — estimates transcript abundance and isoform expression
Salmon (alignment mode) — transcript-level quantification
StringTie — transcript assembly and quantification

You do not need this for gene-level DESeq2 analysis. Keep it if isoform analysis is planned; otherwise it can be deleted to save disk (it’s roughly the same size as the genome BAM).

SampleF1.Chimeric.out.junction

A tab-delimited list of chimeric (fusion) junctions detected by STAR. Each row describes a breakpoint where reads span two genomic locations on different chromosomes, strands, or distant loci — the signature of a gene fusion or translocation. Produced by the –chimSegmentMin 15 and related chimeric options. This is the primary input for STAR-Fusion:

##Full bash script loop to run samples in a directory

#!/bin/bash
set -euo pipefail

STAR_DB="/quobyte/luisccgrp/REFERENCE_DATA/STAR/GRCh38.d1.vd1.Hpylori_gencode.v50"
OUTDIR="./star_output"
mkdir -p "$OUTDIR"

for R1 in ./trimmed_fastq/*_1.trimmed.fastq.gz; do
    R2="${R1%_1.trimmed.fastq.gz}_2.trimmed.fastq.gz"
    SAMPLE=$(basename "${R1%_1.trimmed.fastq.gz}")

    echo "=== Aligning ${SAMPLE} ==="

    STAR \
        --readFilesIn "$R1" "$R2" \
        --outFileNamePrefix "${OUTDIR}/${SAMPLE}." \
        #--outSAMattrRGline ID:"${SAMPLE}" SM:"${SAMPLE}" PL:ILLUMINA LB:lib1 PU:unit1 \ #only necesary if doing variant calling
        --genomeDir "$STAR_DB" \
        --readFilesCommand zcat \
        --runThreadN 16 \
        --twopassMode Basic \
        --outFilterMultimapNmax 20 \
        --alignSJoverhangMin 8 \
        --alignSJDBoverhangMin 1 \
        --outFilterMismatchNmax 999 \
        --outFilterMismatchNoverLmax 0.1 \
        --alignIntronMin 20 \
        --alignIntronMax 1000000 \
        --alignMatesGapMax 1000000 \
        --outFilterType BySJout \
        --outFilterScoreMinOverLread 0.33 \
        --outFilterMatchNminOverLread 0.33 \
        --limitSjdbInsertNsj 1200000 \
        --outSAMstrandField intronMotif \
        --outFilterIntronMotifs None \
        --alignSoftClipAtReferenceEnds Yes \
        --quantMode TranscriptomeSAM GeneCounts \
        --outSAMtype BAM Unsorted \
        --outSAMunmapped Within \
        --genomeLoad NoSharedMemory \
        --chimSegmentMin 15 \
        --chimJunctionOverhangMin 15 \
        --chimOutType Junctions SeparateSAMold WithinBAM SoftClip \
        --chimOutJunctionFormat 1 \
        --chimMainSegmentMultNmax 1 \
        --outSAMattributes NH HI AS nM NM ch

    echo "--- Done: ${SAMPLE} ---"
done

While the program runs, try to understand the above command.

  • Review the alignment statistics in the Log.final.out output file. What proportion of reads mapped unambiguously to the genome database?

Quality control summary

Now that we have gotten to the alignment BAM stage, let’s summarize all the QC outputs from each step: 1. Pre-trim Fastqc 2. Fastp trimming 3. Post-trim fastqc 4. STAR alignment

We could look through each file individually, but multiqc makes it much simpler to see multiple files at once. https://docs.seqera.io/multiqc/

Again the HIVE has a module for multiqc making installation easy

module load multiqc

Running multiqc

mkdir ./multiqc_report
multiqc ./fastqc_results ./trimmed_fastq ./fastqc_results/post_trim ./star_output -o ./multiqc_report

View the html files in your browser to evaluate you samples

Indexing and viewing the contents of you BAM file

The BAM file is a valuable resource that is the basis for many downstream applications. It contains all the alignments for every read in your sample and where they are mapped to on the genome.

Using samtools to view and manipulate BAM files

Samtools is a set of utilities (https://www.htslib.org/doc/samtools.html) that manipulate alignments in the SAM (Sequence Alignment/Map), BAM, and CRAM formats. It converts between the formats, does sorting, merging and indexing, and can retrieve reads in any regions swiftly.

Again, the HIVE has samtools preinstalled, all we need to do is load the module. The first step is to sort and index our BAM files which allows faster processing and search of the BAM file:

module load samtools
cd /quobyte/luisccgrp/USERS/RNAseq_tutorial/<yourname>
cd star_output
# Script to sort all Aligned.out.bam files in a directory
for bam in *.Aligned.out.bam; do
    # Create output filename
    sorted="${bam%.bam}.sorted.bam"
    # Run samtools sort
    samtools sort -o "$sorted" "$bam"
done

# Script to index all Aligned.out.bam files in a directory
for bam in *.Aligned.out.sorted.bam; do samtools index $bam; done

These two commands will produce Aligned.out.sorted.bam and Aligned.out.sorted.bam.bai files for each sample.

Next we can use samtools idxstats which will tell us how many reads are seen per chromosome (or virus/organism)

module load samtools
cd /quobyte/luisccgrp/USERS/RNAseq_tutorial/<yourname>
cd star_output
# Use idxstats to list, then pipe "|" results into the less command
samtools idxstats SampleF1.Aligned.out.sorted.bam | less

You will get an output with 4 columns.
column 1: chromosome or organism

column 2: length of chromosome or organism genome in BP

column 3: number of mapped reads

column 4: number of unmapped reads

Using Integrative Genomics Viewer (IGV) to visualize your alignments

The Integrative Genomics Viewer (IGV) https://igv.org/ is an app on your computer (not server) that you can use to view your BAM file visually.

To install, go to https://igv.org/ and download the IGV desktop application. Do this in your laptop or desktop computer.

It should download the hg38 genome, then you can “Load from File” in the file menu. Go to the HIVE from your laptop (must have run sshfs as above to have access to the HIVE from your laptop)

alt text
alt text