Next Generation Sequencing Analysis
Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
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
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
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
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)