A wet-dry hybrid biologist's take on genetics and genomics. Mostly is about Linux, R, python, reproducible research, open science and NGS. Grab my book to transform yourself to a computational biologist https://divingintogeneticsandgenomics.ck.page/
Recently, there are several papers focusing on batch recomputing the public RNA-seq data sets either by local clusters or remote cloud computing. I will put the links below: 1. A cloud-based workflow to quantify transcript-expression levels in public cancer compendia. one can download the TCGA and CCLE RNA-seq data in transcript-level in counts and TPM units. used gencode v24 as reference and kallisto to calculate. 2. Rapid and efficient analysis of 20,000 RNA-seq samples with Toil. Access the data here. Note that the data are mapped to hg38 reference genome and gencode v23 as annotation. TCGA/ICGC /GTEx and CCLE data. 3. intropolis is a list of exon-exon junctions found across 21,504 human RNA-seq samples on the Sequence Read Archive (SRA) from spliced read alignment to hg19 with Rail-RNA. 4. RESTful RNA-seq Analysis APIA simple RESTful API to access analysis results of all public RNAseq data for nearly 200 species in European Nucleotide Archive. 5. ExpressionAtlas bioconductor package: This package is for searching for datasets in EMBL-EBI Expression Atlas, and downloading them into R for further analysis. Each Expression Atlas dataset is represented as a SimpleList object with one element per platform. Sequencing data is contained in a SummarizedExperiment object, while microarray data is contained in an ExpressionSet or MAList object. 6. Bgeeas suggested in the comment. 7. The Digital Expression ExplorerThe Digital Expression Explorer (DEE) is a repository of digital gene expression profiles mined from public RNA-seq data sets. These data are obtained from NCBI Short Read Archive.
However, I just got several RNA-seq data to play with, and I think it is a good time-point for me to get my hands wet on those RNA-seq quantification tools (especially those alignment-free ones) and get a personal idea of how different tools perform. I am not doing bench-marking, as one should simulate the RNA-seq reads by e.g. polyester to have the ground truth.
Choosing alignment based tools (such as tophat, STAR, bowtie, HISAT) or alignment free ones depends on the purpose of your study. Salmon and kallisto requires the reads "pesudo-map" to the transcriptome, so one has to provide a fasta file containing all the transcripts you want to quantify. Therefore, if you want to find novel transcripts, you probably should go with the alignment based methods. It is also shown recently that Widespread intron retention diversifies most cancer transcriptomes. If you want to do similar things, you need to use mappers with the genome (not the transcriptiome) as a reference.
kallisto can output a pseudo-bam which can be useful for some people. Salmon will have the same functionality in the next release according to Rob.
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.6.0/SalmonBeta-0.6.0_DebianSqueeze.tar.gz
tar xvzf SalmonBeta-0.6.0_DebianSqueeze.tar.gz
To use Salmon in quasi-mapping-based mode, then you first have to build an Salmon index for your transcriptome.Ensemble release 75 is the latest version for GRCh37 (hg19). starting from release 76 it is GRCh38 (hg38)
Download the hg19 version of cDNA and non-coding RNA fasta:
Default index --- The quasi index has been made the default type. This means that it is no longer necessary to provide the --type option to the index command. The fmd index remains enabled, but may be removed in a future version.
salmon index -t Homo_sapiens.GRCh37.75.cdna.ncrna.fa -i Homo_sapiens.GRCh37.75_quasi_index
It finished in a little over 10 mins.
quantify the transcript:
The RNA-seq library I am analyzing is single-end stranded, reads from reverse strand. so I specified -l SR. readsalmon doc for different library types. I especially like the figure representation below:
Salmon can also give gene-level quantification as long as feed a gtf file. In addtion, I specify the fragment length to be 200 bp (default) and standard deviation of 20 (default 80). Details here. For paired-end RNA-seq, the fragment length distribution can be infered from the fastq files, but for single-end data, it needs to be specified.
It is recommended using tximport to get the gene-level quantification. I asked the difference between tximport andsalmon quant -g. The developer of salmon @Rob Patro answered:
Main diffs I can think of (1) in R
(2) integrated with DESeq2
(3) Can derive multi-sample effective gene lengths
kallisto index -i Homo_sapiens.GRCh37.75.cdna.ncrna.kalisto.idx Homo_sapiens.GRCh37.75.cdna.ncrna.fa
It took me around 15 mins to build the index.
quantification:
-l, --fragment-length=DOUBLE Estimated average fragment length -s, --sd=DOUBLE Estimated standard deviation of fragment length(default: value is estimated from the input data)
In the case of single-end reads, the -l option must be used to specify the average fragment length. Typical Illumina libraries produce fragment lengths ranging from 180–200 bp but it’s best to determine this from a library quantification with an instrument such as an Agilent Bioanalyzer.
Common values for single end reads are insert length 200 and sd 20. If you have any better information, like the person who prepped the library or better yet, data from bioanalyzer that will of course be better.
Finished in ~6 mins. again, blazing fast as Salmon does.
counts versus TPM/RPKM/FPKM
I will need to quote from this blog post on explaning technical differences among RPKM, FPKM and TPM.
These three metrics attempt to normalize for sequencing depth and gene length. Here’s how you do it for RPKM:
Count up the total reads in a sample and divide that number by 1,000,000 – this is our “per million” scaling factor.
Divide the read counts by the “per million” scaling factor. This normalizes for sequencing depth, giving you reads per million (RPM)
Divide the RPM values by the length of the gene, in kilobases. This gives you RPKM.
FPKM is very similar to RPKM. RPKM was made for single-end RNA-seq, where every read corresponded to a single fragment that was sequenced. FPKM was made for paired-end RNA-seq
TPM is very similar to RPKM and FPKM. The only difference is the order of operations. Here’s how you calculate TPM:
Divide the read counts by the length of each gene in kilobases. This gives you reads per kilobase (RPK).
Count up all the RPK values in a sample and divide this number by 1,000,000. This is your “per million” scaling factor.
Divide the RPK values by the “per million” scaling factor. This gives you TPM.
So you see, when calculating TPM, the only difference is that you normalize for gene length first, and then normalize for sequencing depth second. However, the effects of this difference are quite profound.
When you use TPM, the sum of all TPMs in each sample are the same. This makes it easier to compare the proportion of reads that mapped to a gene in each sample. In contrast, with RPKM and FPKM, the sum of the normalized reads in each sample may be different, and this makes it harder to compare samples directly.
My colleage @samir processed the same data using STAR-HTseq piepline, and it gives the raw counts for each gene. He is very skepitcal about transcript level quantification and would focus on gene-level for now. What adds the complexity a bit is that he used gencode v19 as annotation and the gene name has a dot + digits in the end of each gene. e.g. ENSG00000000003.10 vs ENSG00000000003 in gtf files downloaded from ensemble. I will need to get rid of the digits in the end.
I will need to convert the raw counts from the STAR-HTseq pipeline to TPM for comparison as Salmon and kallistooutput TPM and estimated counts. Read the post: convert counts to TPM. Kamil Slowikowski wrote a function to convert counts to TPM, and the function involves an effective length of the features.
If we take the fragment length to be fixed, then the effective length is how many fragments can occur in the transcript. This turns out to be length - frag_len +1. The number of fragments coming from a transcript will be proportional to this number, regardless of whether you sequenced one or both ends of the fragment. In turn, when it comes to probabilistically assigning reads to transcripts the effective length plays a similar role again. Thus for short transcripts, there can be quite a difference between two fragment lengths. To go back to your example if you have transcript of length 310, your effective length is 10 (if fragment length is 300) or 160 (if fragment length is 150) in either case, which explains the discrepancy you see.
From @Rob
The effective length is computed by using the fragment length distribution to determine the effective number of positions that can be sampled on each transcript. You can think of this as convolving the fragment length distribution with the characteristic function (the function that simply takes a value of 1) over the transcript. For example if we observe fragments of length 50 --- 1000, a position more than 1000 bases from the end of the transcript will contribute a value of 1 to the effective length, while a position 150 bases will contribute a value of F(150), where F is the cumulative distribution function of the fragment length distribution. For single end data, where we can't learn an empirical FLD, we use a gaussian whose mean and standard deviation can be set with --fldMean and --fldSD respectively.
From Harold Pimentel's post above. He is in Lior Pachter's group.
Effective length refers to the number of possible start sites a feature could have generated a fragment of that particular length. In practice, the effective length is usually computed as:
where uFDL is the mean of the fragment length distribution which was learned from the aligned read. If the abundance estimation method you’re using incorporates sequence bias modeling (such as eXpress or Cufflinks), the bias is often incorporated into the effective length by making the feature shorter or longer depending on the effect of the bias.
R scripts to convert HTseq counts to TPM
library(dplyr)
setwd("/Volumes/mdarisngc03.mdanderson.edu/scratch/RNA-seq-26357338/")
list.files()
counts.3R<- read.table("STAR_3R-30390482_htseq.cnt", sep="\t", header=F, stringsAsFactors=F)
counts.50R<- read.table("STAR_50R-30393469_htseq.cnt", sep="\t", header=F, stringsAsFactors=F)
counts.WT<- read.table("STAR_WT-30393468_htseq.cnt", sep="\t", header=F, stringsAsFactors=F)
counts.df<- cbind(counts.WT, counts.3R, counts.50R)
counts.df<-counts.df[, c(1,2,4,6)]
names(counts.df)<- c("gene_name", "HTseq.WT", "HTseq.3R", "HTseq.50R")
## get rid of the digits in the end of the gene_namecounts.df<-counts.df %>% mutate(gene_name= gsub("\\.[0-9]+", "", gene_name))
To convert the raw counts from HTSeq (gencode v19 as annotation), I will need the length of each gene. Although one can compute the gene length from the gtf files, the gene-level output of Salmon has already computed it for me. I will just need to "borrow" from there.
counts_to_tpm<-function(counts, featureLength, meanFragmentLength) {
# Ensure valid arguments.
stopifnot(length(featureLength) == nrow(counts))
stopifnot(length(meanFragmentLength) == ncol(counts))
# Compute effective lengths of features in each library.effLen<- do.call(cbind, lapply(1:ncol(counts), function(i) {
featureLength-meanFragmentLength[i] +1
}))
# Exclude genes with length less than the mean fragment length.idx<- apply(effLen, 1, function(x) min(x) >1)
counts<-counts[idx,]
effLen<-effLen[idx,]
featureLength<-featureLength[idx]
# Process one column at a time.tpm<- do.call(cbind, lapply(1:ncol(counts), function(i) {
rate= log(counts[,i]) - log(effLen[,i])
denom= log(sum(exp(rate)))
exp(rate-denom+ log(1e6))
}))
# Copy the column names from the original matrix.
colnames(tpm) <- colnames(counts)
return(tpm)
}
featureLength<-counts.df$Length## fragment length is 200bp meanFragmentLength<- c(200, 200, 200)
counts<- as.matrix(counts.df[,c(2,3,4)])
rownames(counts)<-counts.df$gene_nameTPM.from.HTSeq<- counts_to_tpm(counts, featureLength, meanFragmentLength)
Import salmon and kallisto transcript level
Use tximport to import salmon and kallisto transcript level quantification
Kallisto and Salmon seems to have very tight correlation. STAR-HTSeq based RNA-seq pipeline is a bit off when correlated with the other two. Note that I did not do adaptor and low qualit base trimming, STAR may discard some informative reads. In contrast, salmon (and kallisto?) is very robust to quality and adapter trimming. If you have other thoughs, please comment :) I am sorry that my analysis is not fully reproducible as I can not share the data now.
One of the benefits of the quasi-mapping approach taken by sailfish is that it is rather robust to quality and adapter trimming. One way you can assess this is by looking at the mapping rate (i.e. the fraction of reads that are mapped) --- which appears in the comments of the main output file "quant.sf". You may be able to recover a small fraction of extra reads with quality / adapter trimming, but generally this is not necessary for Sailfish to map reads accurately for quantification purposes.
UPDATE 07/24/16 Robs suggested that "taking the aggregate feature length (i.e. the gene length) and then correcting it may not always yield the same result as correcting the feature length (i.e. the transcript lengths) and then aggregating them to the gene level." see comment below.
I then re-calculated TPM from counts using the method he suggested, and the correlation looks better, but spearman correlation does not improve.
get rid of the digits (gene version) in the end for the gene names (gencode v19)
cat STAR_WT-30393468_htseq.cnt| sed -E 's/\.[0-9]+//'> WT_htseq.cnt
add #! /bin/env python to the head of the python script Rob wrote.
#! /bin/env pythonimport argparse
import pandas as pd
import numpy as np
import sys
defmain(args):
gtable = pd.read_table(args.ginput).set_index('Name')
ttable = pd.read_table(args.tinput).set_index('Name')
tgmap = pd.read_table(args.tgmap, names=['t', 'g']).set_index('t')
gene_lengths = {}
j =0# Map over all gene groups (a gene and its associated transcripts)for g, txps in tgmap.groupby('g').groups.iteritems():
if j %500==1:
print("Processed {} genes".format(j))
j +=1# The set of transcripts present in our salmon index
tset = []
for t in txps:
if t in ttable.index:
tset.append(t)
# If at least one of the transcripts was presentiflen(tset) >0:
# The denominator is the sum of all TPMs
totlen = ttable.loc[tset,'TPM'].sum()
# Turn the relative TPMs into a proper partition of unityif totlen >0:
tpm_fracs = ttable.loc[tset, 'TPM'].values / ttable.loc[tset,'TPM'].sum()
else:
tpm_fracs = np.ones(len(tset)) /float(len(tset))
# Compute the gene's effective length as the abundance-weight# sum of the transcript lengths
elen =0.0for i,t inenumerate(tset):
elen += tpm_fracs[i] * ttable.loc[t, 'EffectiveLength']
gene_lengths[g] = elen
# Give the table an effective length field
gtable['EffectiveLength'] = gtable.apply(lambdar : gene_lengths[r.name] if r.name in gene_lengths else1.0, axis=1)
# Write it to the output file
gtable.to_csv(args.output, sep='\t', index_label='Name')
if__name__=="__main__":
parser = argparse.ArgumentParser(description="Compute gene-level effective lengths")
parser.add_argument('--ginput', type=str, help='gene level input table')
parser.add_argument('--tinput', type=str, help='transcript level input table')
parser.add_argument('--tgmap', type=str, help='transcript -> gene mapping')
parser.add_argument('--output', type=str, help='output table with extra column')
args = parser.parse_args()
main(args)
chmod u + x glef.py
/glef.py --ginput WT_htseq.cnt --tinput WT_quant.sf --tgmap tx2gene.tsv --output WT_htseq_effective_length.tsv
compare effecitve length derived from transcript level and salmon output gene level
WT.count.df<- read.table("~/playground/gene-level-effective-length/WT_htseq_effective_length.tsv",
sep="\t", header=T, stringsAsFactors=F)
head(counts.df)
counts.effective.df<- left_join(counts.df, WT.count.df, by=c("gene_name"="Name"))
## scatter plot of effective length from gene level output of salmon and effective length derived from transcript level salmon output
plot(log2(counts.effective.df$EffectiveLength.x+1), log2(counts.effective.df$EffectiveLength.y+1), xlab="log2 effective gene length from salmon gene level output", ylab="log2 effective gene length derived from transcript level")
abline(a=0, b=1, col="red")
