kallisto to tximport to DESeq2

Part 3 of the bulk RNA-Seq workflow

Summarise transcript abundances to gene level, fit the DESeq2 model, and write out an annotated differential expression table.

This walkthrough combines several package vignettes into one continuous workflow, picking up exactly where part 2 left off.

We start with:

a kallisto/ directory containing one folder per sample, each holding an abundance.h5
a samples.txt file describing which sample belongs to which condition
the .gtf annotation file that came in the same tarball as your kallisto index

If you have not installed the R packages yet, part 1 lists every one of them.

About the code on this page

The code blocks below are the commands you run; the grey blocks underneath are the real output from the worked example — six mouse samples, three treated with STLC and three DMSO controls. Your numbers will differ. The shapes of the objects should not.

Load your libraries

library("tximport")
library("readr")
library("rhdf5")
library("GenomicFeatures")
library("DESeq2")

Set up your working directory

Point R at the folder that contains both your kallisto/ directory and your samples.txt file:

setwd("path/to/your/project/")

Have a look at what is in there:

contents <- dir(file.path("./"))
head(contents)

## [1] "0h.csv"              "0h_STLC_vs_DMSO.csv" "0h_txi_kallisto.csv"
## [4] "5190-S1.fastq.gz"    "5190-S2.fastq.gz"    "5190-S3.fastq.gz"

Call this variable contents, not dir. dir() is a base R function, and naming a variable after it is a small trap you will eventually walk into.

Describe your experimental conditions

Read the sample sheet you made at the end of part 2:

samples <- read.table(file.path("./samples.txt"), header = TRUE)
samples

##     sample condition
## 1 S1output      DMSO
## 2 S2output      STLC
## 3 S3output      DMSO
## 4 S4output      STLC
## 5 S5output      DMSO
## 6 S6output      STLC

Name the rows after the samples:

rownames(samples) <- samples$sample
samples

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

Build the full path to every kallisto .h5 file:

files <- file.path("./kallisto", samples$sample, "abundance.h5")
files

## [1] "./kallisto/S1output/abundance.h5" "./kallisto/S2output/abundance.h5"
## [3] "./kallisto/S3output/abundance.h5" "./kallisto/S4output/abundance.h5"
## [5] "./kallisto/S5output/abundance.h5" "./kallisto/S6output/abundance.h5"

Give each path a name, so the columns of the count matrix are labelled:

names(files) <- paste0("sample", 1:nrow(samples))
files

##                            sample1                            sample2
## "./kallisto/S1output/abundance.h5" "./kallisto/S2output/abundance.h5"
##                            sample3                            sample4
## "./kallisto/S3output/abundance.h5" "./kallisto/S4output/abundance.h5"
##                            sample5                            sample6
## "./kallisto/S5output/abundance.h5" "./kallisto/S6output/abundance.h5"

Using 1:nrow(samples) rather than a hard-coded 1:6 means this line keeps working when your next experiment has eight samples.

You could equally use names(files) <- samples$sample, which labels the columns with your real sample names instead of sample1…sample6. That is often clearer when you come back to the results months later.

Before you go further, confirm that files and samples are in the same order. Everything downstream assumes that column n of the count matrix corresponds to row n of the sample table, and nothing will warn you if it does not.

Build a transcript-to-gene map

kallisto quantifies transcripts. Differential expression is nearly always asked at the level of genes. So you need a table mapping every transcript ID to its gene ID, and tximport will use it to sum transcripts up to genes.

Getting this mapping right is important, and slightly fiddly — several published vignettes suggest approaches that do not work in practice. The reliable route: the kallisto authors’ index repository ships a .gtf file alongside each prebuilt index, and that GTF is exactly what you want here, because it is guaranteed to be the annotation your index was built from.

txdb <- makeTxDbFromGFF("path/to/mus_musculus/Mus_musculus.GRCm38.96.gtf")

Have a look at what the object holds:

columns(txdb)

##  [1] "CDSCHROM"   "CDSEND"     "CDSID"      "CDSNAME"    "CDSPHASE"
##  [6] "CDSSTART"   "CDSSTRAND"  "EXONCHROM"  "EXONEND"    "EXONID"
## [11] "EXONNAME"   "EXONRANK"   "EXONSTART"  "EXONSTRAND" "GENEID"
## [16] "TXCHROM"    "TXEND"      "TXID"       "TXNAME"     "TXSTART"
## [21] "TXSTRAND"   "TXTYPE"

keytypes(txdb)

## [1] "CDSID"    "CDSNAME"  "EXONID"   "EXONNAME" "GENEID"   "TXID"     "TXNAME"

Pull out every transcript name, then map each one to its gene:

k <- keys(txdb, keytype = "TXNAME")

tx2gene <- AnnotationDbi::select(txdb, keys = k, columns = "GENEID", keytype = "TXNAME")
head(tx2gene)

## 'select()' returned 1:1 mapping between keys and columns
##               TXNAME             GENEID
## 1 ENSMUST00000193812 ENSMUSG00000102693
## 2 ENSMUST00000082908 ENSMUSG00000064842
## 3 ENSMUST00000192857 ENSMUSG00000102851
## 4 ENSMUST00000161581 ENSMUSG00000089699
## 5 ENSMUST00000192183 ENSMUSG00000103147
## 6 ENSMUST00000193244 ENSMUSG00000102348

Write AnnotationDbi::select, not select

If you have loaded the tidyverse — and sooner or later you will — dplyr::select masks AnnotationDbi::select, and this line fails with an error that does not obviously point at the cause. Naming the package explicitly makes the code work regardless of what else is loaded.

Run tximport

txi.kallisto <- tximport(files, type = "kallisto", tx2gene = tx2gene, ignoreTxVersion = TRUE)

## 1 2 3 4 5 6
## transcripts missing from tx2gene: 1673
## summarizing abundance
## summarizing counts
## summarizing length
## summarizing inferential replicates

Two things in that output are worth understanding

ignoreTxVersion = TRUE strips the version suffix from transcript IDs — ENSMUST00000193812.1 becomes ENSMUST00000193812 — so that IDs from your index match IDs from your GTF even when the two are numbered slightly differently. Without it you can get a mapping rate near zero.

transcripts missing from tx2gene: 1673 is the number of transcripts quantified by kallisto that had no entry in your map, and so were dropped. A number in the hundreds or low thousands out of ~120,000 is normal. A number in the tens of thousands means your index and your GTF do not match, and you should stop and fix that before going any further.

head(txi.kallisto$counts)

##                    sample1   sample2 sample3 sample4   sample5  sample6
## ENSMUSG00000000001    8299 14400.000    9423    9805 12225.000 9895.000
## ENSMUSG00000000003       0     0.000       0       0     0.000    0.000
## ENSMUSG00000000028     792  2803.000    1021    1244  1559.567 1514.813
## ENSMUSG00000000037     194   296.000     283     217   383.000  270.000
## ENSMUSG00000000049       0     0.000       0       0     0.000    0.000
## ENSMUSG00000000056    2401  3793.896    2382    2474  3247.000 2581.000

The counts are not whole numbers because kallisto assigns reads to transcripts probabilistically; a read compatible with three transcripts contributes a fraction to each. DESeq2 expects this and handles it.

Make a DESeq2 object

ddsTxi <- DESeqDataSetFromTximport(txi.kallisto,
                                   colData = samples,
                                   design = ~ condition)
ddsTxi

## class: DESeqDataSet
## dim: 36047 6
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(36047): ENSMUSG00000000001 ENSMUSG00000000003 ...
##   ENSMUSG00000118389 ENSMUSG00000118393
## rowData names(0):
## colnames(6): S1output S2output ... S5output S6output
## colData names(2): sample condition

Prefilter to remove low counts

This step is suggested in the DESeq2 vignette. It removes any gene with fewer than 10 reads across all samples — genes with almost no signal cost you statistical power in multiple-testing correction without ever being detectable.

keep <- rowSums(counts(ddsTxi)) >= 10
dds <- ddsTxi[keep, ]
dds

## class: DESeqDataSet
## dim: 20185 6
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(20185): ENSMUSG00000000001 ENSMUSG00000000028 ...
##   ENSMUSG00000118382 ENSMUSG00000118389
## rowData names(0):
## colnames(6): S1output S2output ... S5output S6output
## colData names(2): sample condition

Prefiltering removed a lot of genes here: 36,047 down to 20,185. That is expected — most annotated genes are not expressed in any given tissue.

Set the control condition

Make sure comparisons are made against the correct baseline. Set ref to whichever level is your control — here, DMSO:

dds$condition <- relevel(dds$condition, ref = "DMSO")

Skip this and R will pick the baseline alphabetically, which is very unlikely to be what you meant. Every fold change on the page would then be reported the wrong way round.

Run DESeq2

dds <- DESeq(dds)
res <- results(dds)
head(res)

## log2 fold change (MLE): condition STLC vs DMSO
## Wald test p-value: condition STLC vs DMSO
## DataFrame with 6 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE       stat    pvalue
##                    <numeric>      <numeric> <numeric>  <numeric> <numeric>
## ENSMUSG00000000001 10528.932     -0.0921937  0.136806 -0.6739009  0.500374
## ENSMUSG00000000028  1402.240      0.3068577  0.290703  1.0555719  0.291164
## ENSMUSG00000000037   293.507      0.0588162  0.630455  0.0932918  0.925672
## ENSMUSG00000000056  2862.085      0.3036562  0.289039  1.0505727  0.293455
## ENSMUSG00000000058   195.899      0.7890995  0.757948  1.0410992  0.297830

Confirm which comparison was actually made:

resultsNames(dds)

## [1] "Intercept"              "condition_STLC_vs_DMSO"

Read that line carefully every time. condition_STLC_vs_DMSO means a positive log2FoldChange is higher in STLC than in DMSO.

Optional: shrink effect sizes for plotting

Genes with low counts produce wildly noisy fold-change estimates. lfcShrink pulls those estimates towards zero in proportion to how uncertain they are, which makes plots far more honest without changing the p-values.

resFC <- lfcShrink(dds, coef = "condition_STLC_vs_DMSO", type = "apeglm")
head(resFC)

The coef argument must exactly match one of the names printed by resultsNames(dds) above. type = "apeglm" requires the apeglm package to be installed.

Shrinkage is not needed for ranking or for text-based analysis, so the result is written to a separate object and the unshrunken res is used below. Use resFC for volcano and MA plots; use res for the results table.

Examine the output

Order the results by p-value:

resOrdered <- res[order(res$pvalue), ]
summary(res)

##
## out of 20185 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 108, 0.54%
## LFC < 0 (down)     : 135, 0.67%
## outliers [1]       : 387, 1.9%
## low counts [2]     : 3884, 19%
## (mean count < 11)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

That summary is worth reading line by line. 243 genes pass an adjusted p-value threshold of 0.1. 387 genes were flagged as outliers by Cook’s distance — usually a single sample with an extreme count — and 3,884 were filtered out by independent filtering because their mean count was too low for the test to have any power.

Count the significant genes yourself:

sum(res$padj < 0.1, na.rm = TRUE)

## [1] 243

na.rm = TRUE is not optional. padj is NA for every gene that was filtered out, and without it the sum returns NA.

Check that nothing is duplicated:

any(duplicated(rownames(res)))
any(duplicated(colnames(res)))

## [1] FALSE
## [1] FALSE

Add gene symbols and descriptions

Ensembl IDs are unambiguous but unreadable. Map them to gene symbols and plain English descriptions so your collaborators can use the table.

library("AnnotationDbi")
library("org.Mm.eg.db")

resOrdered_df <- as.data.frame(resOrdered)
resOrdered_df$ensembl <- rownames(resOrdered_df)

resOrdered_df$symbol <- mapIds(org.Mm.eg.db,
                               keys = resOrdered_df$ensembl,
                               column = "SYMBOL",
                               keytype = "ENSEMBL",
                               multiVals = "first")

resOrdered_df$description <- mapIds(org.Mm.eg.db,
                                    keys = resOrdered_df$ensembl,
                                    column = "GENENAME",
                                    keytype = "ENSEMBL",
                                    multiVals = "first")

head(resOrdered_df)

Annotate with mapIds, not merge

An earlier version of this walkthrough joined the results to an annotation table with merge(). That is a trap, for three reasons, and all three are silent:

merge() performs an inner join, so every gene whose symbol did not map was dropped from the results without warning
where one symbol maps to several Entrez IDs, rows were duplicated
merge() re-sorts its output by the join column, which destroyed the p-value ordering you just created — the resulting file came out in alphabetical order by gene symbol

mapIds adds a column to the existing table instead of rebuilding it. Every gene is kept, unmapped genes get NA, the row count is unchanged, and the ordering survives. multiVals = "first" states explicitly what to do when an Ensembl ID maps to more than one symbol.

Use org.Hs.eg.db for human, org.Dm.eg.db for fly, and so on.

Write your results to a file

write.csv(resOrdered_df, file = "STLC_vs_DMSO_annotated.csv")

Your results are now readable in R, on the command line, in Python, and in Excel.

Record what you ran

Before you close R, capture the versions of everything you used. In six months, when a reviewer asks, this is the only thing that will let you answer:

sessionInfo()

Save that output alongside your results file.