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=32g -c 16 --partition=high --time=03:00:00 --account=luisccgrp

This opens up a 3 hour session with 16 CPUS/threads and up to 32GB of memory. 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

Learning to use symbolic links with raw data

In the LCC lab we normally store all raw FASTQ data in the /quobyte/luisccgrp/RAW_DATA/ folder (based on whether it is RNAseq, Exome, WGS etc). When aligning and analyzing our data we put it into a corresponding fold in the /quobyte/luisccgrp/SEQ_DATA folder. It’s most convenient to have all our sample FASTQs in one folder for alignment, but would use excessive disk space and time to copy all the FASTQ into a new folder. Therefore, we often use symbolic links that create a symbolic file that links back to the original file in the RAW_DATA folder. We’ll do a similar symbolic link into our sample folder.

#Go to your folder
cd /quobyte/luisccgrp/USERS/RNAseq_tutorial/<yourfolder>
ls
#I moved your files into "full_sample" folder and create a short version for us to use today. Look at these and the sizes of the files.  
ls full_sample
du -h full_sample/*.gz
ls full_sample/shortsample
du -h full_sample/shortsample/*.gz

Now let’s create the symbolic links.

#create a "sample" folder and go into it
mkdir sample
cd sample
#Create symbolic links to all *.gz files in a particular folder to our current one
ln -s ../full_sample/shortsample/*.gz ./
#take a look at how the link works
ls
ls -l *.gz
du -h *.gz

Note: Some algorithms don’t like to take in symbolic links and need the actual file. ### Working with FASTQ files

Let’s look at one of your FASTQ files: I’ll be looking at my file SampleF1_1.fq.gz , which comprise the sequence reads for SampleF1 in FASTQ format.

Note: This file has been compressed with GZIP, so you may need to decompress it using gzip or zcat to access the FASTQ text contained. Some programs (such as the zless text viewer, and many bioinformatic tools) know how to deal with gzipped files.

### view the contents in a file
zless Data/SampleF1_1.fq.gz

### check the total number of lines
zless Data/SampleF1_1.fq.gz | wc -l

#Count the number of lines in one of your shortsamples
zcat full_sample/shortsample/SampleF1_1.fq.gz | wc -l 

# This is how you would unzip the gzipped file, DON'T NEED TO DO THIS
#gunzip SampleF1_1.fq.gz

How many reads are in your file? - How many lines does each file have? How many reads do we have for each sample? - How long is the first read for SampleF1_1.fq.gz ? - How to tell if this is paired Read 1 or paired Read 2? - What is the ReadGroupID, FlowCell and Lane number?

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 l  #shows only modules starting with the letter "L"

### 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

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. The HIVE cluster also provides this as a module ### Installing fastp

module load fastp
fastp --help

To filtering low-quality reads, we use the following command:

#make a directory for the new trimmed FASTQ
mkdir trimmed_fastq
fastp \
  -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

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

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/ \
#--outSAMattrRGline ID:"${SAMPLE}" SM:"${SAMPLE}" PL:ILLUMINA LB:lib1 PU:unit1 \ #only necesary if doing variant calling
--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

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