Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown

Article reproduced here

System setup

The tutorial presumes that the analysis will be run on a 64 bit Linux system. This particular work flow has been done on Ubuntu 22.04.

Software required

• Hisat2
• Stringtie
• samtools

For installation, unpack each of the software above (individually), then do the following in the terminal for each of them:
sudo export PATH=$PATH:path/to/extracted/folderORexecutibleFiles
It may be necessary to add the line export PATH=$PATH:path/to/extracted/folderORexecutibleFiles to .bashrc file in the home directory.

Setting up R environment

R and Rstudio setup sudo apt install r-base

Download and install Rstudio from here: www.rstudio.com/ide/download/desktop
For Ubuntu 22.04, download and install: https://download1.rstudio.org/electron/jammy/amd64/rstudio-2023.12.1-402-amd64.deb

Before opening RStudio, open the terminal and install the following applications in Ubuntu:
sudo apt install build-essential libcurl4-openssl-dev libxml2-dev libssl-dev gfortran

Start Rstudio and enter the following commands:

>source("http://www.bioconductor.org/biocLite.R")
>biocLite(c("alyssafrazee/RskittleBrewer","ballgown","genefilter","dplyr","devtools"))

Note: If there is any error in the installation, pay close attention to the error message and install any unmet dependency

Files to be downloaded

In this paper the reference genome is the Human X chromosome. Only X chromosome is being used as a reference genome in order to reduce the download and analysis time. The user may also use the entire genome of the relevant organism. The reference genome may also be downloaded from the Hisat2 website since it is pre-indexed.

https://daehwankimlab.github.io/hisat2/download/

This link contains the reference genome as well as the sample fastq files. The following workflow includes files from the above archive. The data may be replaced with the relevant data as well.

Protocol

HISAT2: Alignment of RNA-Seq reads to the genome

Map the reads from each sample to the reference genome:

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188044_chrX_1.fastq.gz -2 chrX_data/samples/ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam

The explanation for each part of the command is as follows:

hisat2 : calls the hisat2 program. This will work only if the hisat2 folder has been added to the PATH properly
-p 8 : "-p" specifies how many threads should be allocated to the command. Here, 8 threads have been allocated. Should always be equal or less than the available threads on the system
-dta : reports alignments tailored for transcript assemblers
-x : Name (or path to name). hisat2 index file comes right after this chrX_data/indexes/chrX_tran : Path to the index files. "chrX_tran" is the common name of all the reference file in the index folder. e.g., the files in the folder are chrX_tran.1.ht2, chrX_tran.2.ht2, chrX_tran.3.ht2, etc
-1 : Has to be put before the forward read of the sample
chrX_data/samples/ERR188044_chrX_1.fastq.gz : Forward read of sample
-2 : Has to be put before the reverse read of the sample
chrX_data/samples/ERR188044_chrX_2.fastq.gz : Reverse read of sample
-S : Default sam output
ERR188044_chrX.sam : File for SAM output
Likewise, the same commands have to be issued for all the input files, as follows:

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188044_chrX_1.fastq.gz -2 chrX_data/samples/ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188104_chrX_1.fastq.gz -2 chrX_data/samples/ERR188104_chrX_2.fastq.gz -S ERR188104_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188234_chrX_1.fastq.gz -2 chrX_data/samples/ERR188234_chrX_2.fastq.gz -S ERR188234_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188245_chrX_1.fastq.gz -2 chrX_data/samples/ERR188245_chrX_2.fastq.gz -S ERR188245_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188257_chrX_1.fastq.gz -2 chrX_data/samples/ERR188257_chrX_2.fastq.gz -S ERR188257_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188273_chrX_1.fastq.gz -2 chrX_data/samples/ERR188273_chrX_2.fastq.gz -S ERR188273_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188337_chrX_1.fastq.gz -2 chrX_data/samples/ERR188337_chrX_2.fastq.gz -S ERR188337_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188383_chrX_1.fastq.gz -2 chrX_data/samples/ERR188383_chrX_2.fastq.gz -S ERR188383_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188401_chrX_1.fastq.gz -2 chrX_data/samples/ERR188401_chrX_2.fastq.gz -S ERR188401_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188428_chrX_1.fastq.gz -2 chrX_data/samples/ERR188428_chrX_2.fastq.gz -S ERR188428_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188454_chrX_1.fastq.gz -2 chrX_data/samples/ERR188454_chrX_2.fastq.gz -S ERR188454_chrX.sam

hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR204916_chrX_1.fastq.gz -2 chrX_data/samples/ERR204916_chrX_2.fastq.gz -S ERR204916_chrX.sam

Samtools: Convert the output sam files to bam

samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam

Usage: samtools [options]

samtools: invokes the program
-@ : Set number of sorting and compression threads
8 : Here we are using 8 threads since this is being done on an 8 thread system
-o : output to file following this
ERR188044_chrX.bam : Output file
ERR188044_chrX.sam : Input file
Likewise, the same commands have to be issued for all the input files, as follows:

samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
samtools sort -@ 8 -o ERR188104_chrX.bam ERR188104_chrX.sam
samtools sort -@ 8 -o ERR188234_chrX.bam ERR188234_chrX.sam
samtools sort -@ 8 -o ERR188245_chrX.bam ERR188245_chrX.sam
samtools sort -@ 8 -o ERR188257_chrX.bam ERR188257_chrX.sam
samtools sort -@ 8 -o ERR188273_chrX.bam ERR188273_chrX.sam
samtools sort -@ 8 -o ERR188337_chrX.bam ERR188337_chrX.sam
samtools sort -@ 8 -o ERR188383_chrX.bam ERR188383_chrX.sam
samtools sort -@ 8 -o ERR188401_chrX.bam ERR188401_chrX.sam
samtools sort -@ 8 -o ERR188428_chrX.bam ERR188428_chrX.sam
samtools sort -@ 8 -o ERR188454_chrX.bam ERR188454_chrX.sam
samtools sort -@ 8 -o ERR204916_chrX.bam ERR204916_chrX.sam

Stringtie: Assemble and quantify expressed genes and transcripts

As per the Stringtie website:

StringTie is a fast and highly efficient assembler of RNA-Seq alignments into potential transcripts. It uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus. Its input can include not only the alignments of raw reads used by other transcript assemblers, but also alignments longer sequences that have been assembled from those reads.In order to identify differentially expressed genes between experiments, StringTie's output can be processed by specialized software like Ballgown, Cuffdiff or other programs (DESeq2, edgeR, etc.)

Usage: stringtie <input.bam ..> [-G ] [-l ] [-o ] [-p ]

stringtie -p 8 -G chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf ERR188044_chrX.bam

Stringtie : Invokes the program
-p : the integer following this specifies the number of threads to assign to the program
8 : 8 threads have been assigned here
-G : Specifies the reference annotation to use for guiding the assembly process (GTF/GFF3)
chrX_data/genes/chrX.gtf : Path to the reference genome here
-o : output file name for the merged transcripts GTF
ERR188044_chrX.bam : This is the input file here. As can be observed in the usage syntax, these arguments may be rearranged.
Rest of the files are processed as follows:

stringtie ERR188044_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf
stringtie ERR188104_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188104_chrX.gtf
stringtie ERR188234_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188234_chrX.gtf
stringtie ERR188245_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188245_chrX.gtf
stringtie ERR188257_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188257_chrX.gtf
stringtie ERR188273_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188273_chrX.gtf
stringtie ERR188337_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188337_chrX.gtf
stringtie ERR188383_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188383_chrX.gtf
stringtie ERR188401_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188401_chrX.gtf
stringtie ERR188428_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188428_chrX.gtf
stringtie ERR188454_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR188454_chrX.gtf
stringtie ERR204916_chrX.bam -p 8 -G chrX_data/genes/chrX.gtf -o ERR204916_chrX.gtf

Stringtie: Merge transcripts from all samples

stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o stringtie_merged.gtf chrX_data/mergelist.txt

Here, stringtie creates a file stringtie_merged.gtf using the mergelist.txt file. mergelist.txt has the list of all the sample.gtf files. sample.gtf files were generated in the previous step by the --merge argument. Paths are necessary if the gtf files are not in the working directory. mergelist.txt looks like this in our experiment:

ERR188044_chrX.gtf
ERR188104_chrX.gtf
ERR188234_chrX.gtf
ERR188245_chrX.gtf
ERR188257_chrX.gtf
ERR188273_chrX.gtf
ERR188337_chrX.gtf
ERR188383_chrX.gtf
ERR188401_chrX.gtf
ERR188428_chrX.gtf
ERR188454_chrX.gtf
ERR204916_chrX.gtf

R- For Differential expression analysis

Load the required libraries by entering the following commands at the R prompt

>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)

You should also have a phenotype file in csv format which specifies information about your sample. The example data includes a phenotype file (geuvadis_phenodata.csv) with the following content:	ids	sex	population
ERR188044	male	YRI
ERR188104	male	YRI
ERR188234	female	YRI
ERR188245	female	GBR
ERR188257	male	GBR
ERR188273	female	YRI
ERR188337	female	GBR
ERR188383	male	GBR
ERR188401	male	GBR
ERR188428	female	GBR
ERR188454	male	YRI
ERR204916	female	YRI

Since we already have the expression values of various genes for these samples, we will be using R to see which genes are differentially regulated based on several parameters. Here we may choose to compare them on the basis of "population" or "sex".

pheno_data <- read.csv("geuvadis_phenodata.csv")

Next, we read in the expression data that has been calculated by StringTie bg_chrX = ballgown(dataDir <- "ballgown", samplePattern = "ERR", pData=pheno_data)

Removing low-abundance genes

Often with RNA-seq data, we may get genes which have too few or even zero counts. Here we may remove all transcripts with a variance across sample less than one.
bg_chrX_filt = subset(bg_chrX, "rowVars(texpr(bg_chrX)) >1", genomesubset=TRUE)

Finding genes with significantly different expression

Here we identify transcripts that have statistically significant difference in expression between groups. Here we see if there is an expression difference on the basis of sex. We can do so using stattest function of Ballgown.

results_transcripts <- stattest(bg_chrX_filt, feature= "transcript",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")

We use "getFC=TRUE" parameter so that we can look at the confounder-adjusted fold change between the two groups. Next we identity genes which show statistical significant difference groups Identify genes that show statistically significant differences between groups. For this we can run the same function that we used to identify differentially expressed transcripts, but here we set feature="gene" in the following command:

results_genes = stattest(bg_chrX_filt, feature="gene", covariate="sex", adjustvars = c("population"), getFC=TRUE, meas="FPKM")

Add gene names and gene IDs to the results_transcripts data frame

results_transcripts←data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)

Sort the results from the smallest P value to the largest

results_transcripts = arrange(results_transcripts, pval) results_genes = arrange(results_genes, pval)