Raw sequencing reads to kallisto pseudoalignment

Part 2 of the bulk RNA-Seq workflow

Quality control, adapter trimming, and pseudoalignment — from fastq.gz files to per-sample transcript abundances.

This is the first part of a workflow that takes you through a complete RNA-Seq analysis. Here we take raw sequencing reads, process them, and map them to gene transcripts. In part 3 we take the transcript abundances for each sample and each condition and combine them to perform differential expression analysis.

Everything is commented so that someone with little or no programming experience can run a robust RNA-Seq analysis. If you would rather not write code at all, Galaxy and UseGalaxy are an open-source, web-based alternative.

Install first, analyse second

Every program used below is installed in part 1: setting up your machine. Do that once, in one sitting, before you start here. Stopping halfway through an analysis to install something is how mistakes get made.

Make sure you can run bash

macOS and Linux — you already have a bash-compatible shell. Open Terminal and you are ready to go.

Windows — you need WSL (Windows Subsystem for Linux). Enable WSL and install Ubuntu by following Canonical’s instructions, then type bash at the command line to open the Linux shell.

The WSL install can take up to an hour on a laptop. Time for a coffee break.

Set up conda and bioconda

An effective way to install and run bioinformatics packages is conda, using the bioconda channel. If you do not already have conda, install it as described in part 1.

Configure the channels. The order matters, because --add puts each new channel at the top of the priority list — so run them in this order and you end up with conda-forge highest, then bioconda, then defaults, which is what bioconda expects:

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

Then update:

conda update --all

conda update --all needs both dashes. conda update all is interpreted as “update a package literally named all” and will simply tell you no such package exists.

Create a conda environment for this project

Notice that your prompt begins with (base) or something like it. That tells you which environment you are in. Make a new environment for each workflow, so that the packages you install are specific to that workflow and do not conflict with versions used elsewhere.

conda create --name RNA
conda activate RNA

You will need to run conda activate RNA every time you open a new terminal session to work on this analysis. To check where you are, run conda env list — the * marks the active environment.

Download your sequencing files and run FastQC

Download your sequencing files — they are almost certainly fastq.gz — and ask your sequencing core or service what processing they already performed. They have very likely demultiplexed them, and they may have trimmed them. If they did nothing, you will need to run quality checks and trim sequencing adapters yourself.

Run FastQC on each of your samples. The output tells you about sequencing quality and the presence or absence of adapter sequences, and from it you can tell what you need to do next.

If any samples are of poor quality you will need to decide whether to keep them or throw them out. Then you can move forward with processing.

Trim adapters with Trim Galore

Trim Galore depends on cutadapt, and optionally on FastQC so it can run a quality report after trimming. Install all three from bioconda, inside your active RNA environment:

conda install -c bioconda cutadapt fastqc trim-galore

Check the versions to confirm they installed:

cutadapt --version
fastqc -v
trim_galore --version

Run Trim Galore with default options, and run FastQC afterwards:

trim_galore --fastqc path/to/your/sample.fastq.gz

Once you are comfortable, a short for loop will run Trim Galore across every file in a folder. Depending on file size and sample number this can take several hours, so start it before you go home.

Concatenate samples run across multiple lanes

If you mitigated batch effects by splitting samples across multiple lanes, you need to concatenate those runs into a single file per sample:

cat 5190-S17* > S17_DMSO_6h.fastq.gz

You can also concatenate from one folder and write the output to another:

cat ./Aaron_data_copy/Aaron_5190_181108B1/5190-S5_S105_L00* > ./Untrimmed/5190-S5.fastq.gz

The * is a wildcard: it matches any sequence of characters at that position. So 5190-S5_S105_L00* means “every file whose name starts with 5190-S5_S105_L00, whatever comes after”. This is a compact way to say “all the lane files belonging to sample S5”.

Keep your files gzipped, and keep the .fastq.gz extension. cat works correctly on gzipped files — concatenating two gzip streams produces a valid gzip file — but downstream tools decide how to read a file partly from its name.

Build or download a kallisto index

kallisto runs in two steps: first you generate an index from a transcriptome, then you pseudoalign each sample’s reads against that index.

Install kallisto if it is not already present:

conda install -c bioconda kallisto
kallisto

Running kallisto with no arguments prints the version and the list of available commands, which is a quick way to confirm the install worked.

For common organisms you can skip building an index entirely, because the authors maintain a repository of prebuilt ones. Download the tarball for your species.

Use the prebuilt tarball if your organism is on the list

The tarball contains both the .idx index and the .gtf annotation the index was built from. In part 2 you need a transcript-to-gene map, and it must come from the same annotation that built the index. Taking both files from the same tarball guarantees that. Mixing an index from one annotation release with a GTF from another is the single most common way this workflow goes quietly wrong — you get results, they are just missing genes.

If your organism is not on the list, build the index yourself following the kallisto manual, and keep the GTF you used.

Extract the tarball, naming the destination folder:

tar -xvzf path/to/downloaded/FILE-NAME.tar.gz -C path/to/destination/folder

The extracted folder now contains your .idx file, which you can point kallisto at directly, and the .gtf file you will need in part 2.

Run kallisto

Now run kallisto on each sample. It is worth reading the manual before you start; what follows is enough to get you going.

kallisto assumes paired-end reads by default, taking two files per sample and computing the fragment length distribution itself. For single-end reads you must pass --single and supply the mean fragment length -l and its standard deviation -s.

-l and -s are fragment length, not read length

These two values describe the sequencing library — the distribution of cDNA fragment sizes that went into the sequencer — not the length of the reads you got out. FastQC cannot tell you them. Get them from your sequencing core, or from the Bioanalyzer or TapeStation trace of your library, remembering to subtract the adapter length from the reported fragment size.

Getting these wrong biases the effective length correction, which biases every abundance estimate. If you genuinely cannot obtain them, say so in your methods and treat small fold changes with suspicion.

Navigate to the directory containing your samples. Here we run kallisto on single-end reads with 100 bootstrap replicates (-b):

kallisto quant -i mouse.idx -o ./kallisto/S1output -b 100 --single -l 200 -s 30 ./S1.fastq.gz

For paired-end reads the command is simpler — no --single, no -l, no -s, and both read files given in order:

kallisto quant -i mouse.idx -o ./kallisto/S1output -b 100 ./S1_R1.fastq.gz ./S1_R2.fastq.gz

Give every sample its own output directory

kallisto always writes files called abundance.h5, abundance.tsv and run_info.json. If you point two samples at the same -o directory, the second silently overwrites the first.

If you write a for loop to iterate over samples, make sure the output path changes each time:

-o ./kallisto/S1output
-o ./kallisto/S2output
and so on

The naming used here — a kallisto/ folder containing one directory per sample, named S1output, S2output … — is what part 3 expects. Use it and the next part will work without modification.

Describe your experiment before you leave the shell

Part 3 needs a plain text file that says which sample is which condition. Make it now, while the experiment is fresh in your mind, and save it as samples.txt alongside your kallisto/ folder:

sample  condition
S1output    DMSO
S2output    STLC
S3output    DMSO
S4output    STLC
S5output    DMSO
S6output    STLC

The columns are whitespace-separated, the header row is required, and the values in the sample column must exactly match your kallisto output directory names.

Once you have a directory per sample and a samples.txt that describes them, you are ready for part 3.