Normalization is essential for RNAseq analysis. However, one needs to understand the underlining assumptions for each methods. Most methods assume there is no global changes between conditions (e.g. TMM normalization). However, this may not be true when global effect occurs. For example, if you delete a gene that controls transcription, you expect to see global gene expression reduction. In that case, other normalization methods need to be considered. (e.g. spike-in controls). The same principle applies to other high-throughput sequencing data such as ChIPseq.
read this very important paper by Rafael A Irizarry: Genome-wide repressive capacity of promoter DNA methylation is revealed through epigenomic manipulation
DESseq2 normalization by Simon Anders:
To estimate the library size, simply taking the total number of (mapped or unmapped) reads is, in our experience, not a good idea. Sometimes, a few very strongly expressed genes are differentially expressed, and as they make up a good part of the total counts, they skew this number. After you divide by total counts, these few strongly expressed genes become equal, and the whole rest looks differentially expressed.
The following simple alternative works much better:
- Construct a "reference sample" by taking, for each gene, the geometric mean of the counts in all samples.
- To get the sequencing depth of a sample relative to the reference, calculate for each gene the quotient of the counts in your sample divided by the counts of the reference sample. Now you have, for each gene, an estimate of the depth ratio.
- Simply take the median of all the quotients to get the relative depth of the library.
This is what the
estimateSizeFactorsfunction of our DESeq package doese.
If one wants to use a set of genes that are not affected by the global change, do
dds = newCountDataSet(CountTable, Design$condition )
dds <- estimateSizeFactors(dds,
controlGenes = rownames(dds) %in% norm_genes)
dds_global <- estimateSizeFactors(dds)
dds_global <- DESeq(dds_global)
res_global <- results(dds_global)
or give self-defined size factors.
sizeFactors(dds) = c(my_Values)
Two nature protocols for RNA-seq analysis
Count-based differential expression analysis of RNA sequencing data using R and Bioconductor Based on DESeq and EdgeR.
Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks
A nice tutorial from f1000 research RNA-Seq workflow: gene-level exploratory analysis and differential expression from Michael Love who is the author of DESeq2.
f1000 bioconductor workflow: RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR
f1000 From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline by Gordon Smith.
A post from Nextgeneseek
The three papers kind of replaces earlier tools from Salzberg’s group (Bowtie/TopHat,Cufflinks, and Cuffmerge)
they offer a totally new way to go from raw RNA-seq reads to differential expression analysis:
align RNA-seq reads to genome (HISATinstead of Bowtie/TopHat, STAR),
assemble transcripts and estimate expression (StringTie instead of Cufflinks), and
perform differential expression analysis (Ballgown instead of Cuffmerge).
Simulation-based comprehensive benchmarking of RNA-seq aligners A nature method paper.
We found that performance varied by genome complexity, and accuracy and popularity were poorly correlated. The most widely cited tool underperforms for most metrics, particularly when using default settings
RapMap: A Rapid, Sensitive and Accurate Tool for Mapping RNA-seq Reads to Transcriptomes. From Sailfish group.
For mapping based methods, usually the raw reads are mapped to transcriptome or genome (need to model gaps by exon-exon junction), and then a gene/transcript level counts are obtained by:
featureCounts, integrate HTseq.Finally, differential expression is carried out by
EBseq An R package for gene and isoform differential expression analysis of RNA-seq data
JunctionSeq differential usage of exons and splice junctions in High-Throughput, Next-Generation RNA-Seq datasets. The methodology is heavily based on the DEXSeq bioconductor package.The core advantage of JunctionSeq over other similar tools is that it provides a powerful automated tools for generating readable and interpretable plots and tables to facilitate the interpretation of the results. An example results report is available here.
MetaSeq Meta-analysis of RNA-Seq count data in multiple studies
derfinder Annotation-agnostic differential expression analysis of RNA-seq data at base-pair resolution
DGEclust is a program for clustering and differential expression analysis of expression data generated by next-generation sequencing assays, such as RNA-seq, CAGE and others
Degust: Perform RNA-seq analysis and visualisation. Simply upload a CSV file of read counts for each replicate; then view your DGE data.
Vennt Dynamic Venn diagrams for Differential Gene Expression.
GlimmaInteractive HTML graphics for RNA-seq data.
bcbio.rnaseq
RNAseqGUI. I have used several times. looks good.
compcodeR
paper: Benchmark Analysis of Algorithms for Determining and Quantifying Full-length mRNA Splice Forms from RNA-Seq Data
paper: Comparative evaluation of isoform-level gene expression estimation algorithms for RNA-seq and exon-array platforms
paper:A benchmark for RNA-seq quantification pipelines
It provides a method to detect changes in the relative abundance of the alternative transcripts (isoforms) of genes. This is called Differential Transcript Usage (DTU).
Detecting DTU is supplementary to the quantification of transcripts by tools like Salmon, Sailfish and Kallisto and the detection of Differential Transcript Expression (DTE) by tools such as Sleuth.
I particularly like the figure in the tutorial showing the differences among DTU, DTE and DEG. The paper transcript-level estimates improve gene-level inferences above also talks about the differences:
! MATS is a computational tool to detect differential alternative splicing events from RNA-Seq data. The statistical model of MATS calculates the P-value and false discovery rate that the difference in the isoform ratio of a gene between two conditions exceeds a given user-defined threshold.
A biostar post: Do not feed rounded estimates of gene counts from kallisto into DESeq2 (please make sure you read through all the comments, and now there is a suggested workflow for feeding rounded estimates of gene counts to DESeq etc)
There is some confusion in the answers to this question that hopefully I can clarify with the three comments below:
- kallisto produces estimates of transcript level counts, and therefore to obtain an estimate of the number of reads from a gene the correct thing to do is to sum the estimated counts from the constituent transcripts of that gene. Of note in the language above is the word "estimate", which is necessary because in many cases reads cannot be mapped uniquesly to genes. However insofar as obtaining a good estimate, the approach of kallisto (and before it Cufflinks, RSEM, eXpress and other "transcript level quantification tools") is superior to naïve "counting" approaches for estimating the number of reads originating from a gene. This point has been argued in many papers; among my own papers it is most clearly explained and demonstrated in Trapnell et al. 2013.
- Although estimated counts for a gene can be obtained by summing the estimated counts of the constituent transcripts from tools such as kallisto, and the resulting numbers can be rounded to produce integers that are of the correct format for tools such as DESeq, the numbers produced by such an approach do not satisfy the distributional assumptions made in DESeq and related tools. For example, in DESeq2, counts are modeled "as following a negative binomial distribution". This assumption is not valid when summing estimated counts of transcripts to obtain gene level counts, hence the justified concern of Michael Love that plugging in sums of estimated transcript counts could be problematic for DESeq2. In fact, even the estimated transcript counts themselves are not negative binomial distributed, and therefore also those are not appropriate for plugging into DESeq2. His concern is equally valid with many other "count based" differential expression tools.
- Fortunately there is a solution for performing valid statistical testing of differential abundance of individual transcripts, namely the method implemented in sleuth. The approach is described here. To test for differential abundance of genes, one must first address the question of what that means. E.g. is a gene differential if at least one isoform is? or if all the isoforms are? The tests of sleuth are performed at the granularity of transcripts, allowing for downstream analysis that can capture the varied questions that might make biological sense in specific contexts.
In summary, please do not plug in rounded estimates of gene counts from kallisto into DESeq2 and other tools. While it is technically possible, it is not statistically advisable. Instead, you should use tools that make valid distributional assumptions about the estimates.
However, Charlotte Soneson, Mike Love and Mark Robinson showed in a f1000 paper: Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences that rounded values from transcript level can be fed into DESeq2 etc for gene-level differential expression, and it is valid and preferable in many ways.
Thanks Rob Patro for pointing it out!
Kallisto and caching everything we possibly can.Isolator has a particular focus on producing stable, consistent estimates. Maximum likelihood approaches produce unstable point estimates: small changes in the data can result in drastically different results, conflating downstream analysis like clustering or PCA. Isolator produces estimates that are in general, simultaneously more stable and more accurate other methods
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.
hg38 and gencode v21 as annotation.Considerable differences are found between the methods in terms of the number and characteristics of the genes that are called differentially expressed. Pre-filtering of lowly expressed genes can have important effects on the results, particularly for some of the methods originally developed for analysis of bulk RNA-seq data. Generally, however, methods developed for bulk RNA-seq analysis do not perform notably worse than those developed specifically for scRNA-seq.
325 followers · starred May 2018
283 followers · starred Nov 2019
Normalization is essential for RNAseq analysis. However, one needs to understand the underlining assumptions for each methods. Most methods assume there is no global changes between conditions (e.g. TMM normalization). However, this may not be true when global effect occurs. For example, if you delete a gene that controls transcription, you expect to see global gene expression reduction. In that case, other normalization methods need to be considered. (e.g. spike-in controls). The same principle applies to other high-throughput sequencing data such as ChIPseq.
read this very important paper by Rafael A Irizarry: Genome-wide repressive capacity of promoter DNA methylation is revealed through epigenomic manipulation
DESseq2 normalization by Simon Anders:
To estimate the library size, simply taking the total number of (mapped or unmapped) reads is, in our experience, not a good idea. Sometimes, a few very strongly expressed genes are differentially expressed, and as they make up a good part of the total counts, they skew this number. After you divide by total counts, these few strongly expressed genes become equal, and the whole rest looks differentially expressed.
The following simple alternative works much better:
- Construct a "reference sample" by taking, for each gene, the geometric mean of the counts in all samples.
- To get the sequencing depth of a sample relative to the reference, calculate for each gene the quotient of the counts in your sample divided by the counts of the reference sample. Now you have, for each gene, an estimate of the depth ratio.
- Simply take the median of all the quotients to get the relative depth of the library.
This is what the
estimateSizeFactorsfunction of our DESeq package doese.
If one wants to use a set of genes that are not affected by the global change, do
dds = newCountDataSet(CountTable, Design$condition )
dds <- estimateSizeFactors(dds,
controlGenes = rownames(dds) %in% norm_genes)
dds_global <- estimateSizeFactors(dds)
dds_global <- DESeq(dds_global)
res_global <- results(dds_global)
or give self-defined size factors.
sizeFactors(dds) = c(my_Values)
Two nature protocols for RNA-seq analysis
Count-based differential expression analysis of RNA sequencing data using R and Bioconductor Based on DESeq and EdgeR.
Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks
A nice tutorial from f1000 research RNA-Seq workflow: gene-level exploratory analysis and differential expression from Michael Love who is the author of DESeq2.
f1000 bioconductor workflow: RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR
f1000 From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline by Gordon Smith.
A post from Nextgeneseek
The three papers kind of replaces earlier tools from Salzberg’s group (Bowtie/TopHat,Cufflinks, and Cuffmerge)
they offer a totally new way to go from raw RNA-seq reads to differential expression analysis:
align RNA-seq reads to genome (HISATinstead of Bowtie/TopHat, STAR),
assemble transcripts and estimate expression (StringTie instead of Cufflinks), and
perform differential expression analysis (Ballgown instead of Cuffmerge).
Simulation-based comprehensive benchmarking of RNA-seq aligners A nature method paper.
We found that performance varied by genome complexity, and accuracy and popularity were poorly correlated. The most widely cited tool underperforms for most metrics, particularly when using default settings
RapMap: A Rapid, Sensitive and Accurate Tool for Mapping RNA-seq Reads to Transcriptomes. From Sailfish group.
For mapping based methods, usually the raw reads are mapped to transcriptome or genome (need to model gaps by exon-exon junction), and then a gene/transcript level counts are obtained by:
featureCounts, integrate HTseq.Finally, differential expression is carried out by
EBseq An R package for gene and isoform differential expression analysis of RNA-seq data
JunctionSeq differential usage of exons and splice junctions in High-Throughput, Next-Generation RNA-Seq datasets. The methodology is heavily based on the DEXSeq bioconductor package.The core advantage of JunctionSeq over other similar tools is that it provides a powerful automated tools for generating readable and interpretable plots and tables to facilitate the interpretation of the results. An example results report is available here.
MetaSeq Meta-analysis of RNA-Seq count data in multiple studies
derfinder Annotation-agnostic differential expression analysis of RNA-seq data at base-pair resolution
DGEclust is a program for clustering and differential expression analysis of expression data generated by next-generation sequencing assays, such as RNA-seq, CAGE and others
Degust: Perform RNA-seq analysis and visualisation. Simply upload a CSV file of read counts for each replicate; then view your DGE data.
Vennt Dynamic Venn diagrams for Differential Gene Expression.
GlimmaInteractive HTML graphics for RNA-seq data.
bcbio.rnaseq
RNAseqGUI. I have used several times. looks good.
compcodeR
paper: Benchmark Analysis of Algorithms for Determining and Quantifying Full-length mRNA Splice Forms from RNA-Seq Data
paper: Comparative evaluation of isoform-level gene expression estimation algorithms for RNA-seq and exon-array platforms
paper:A benchmark for RNA-seq quantification pipelines
It provides a method to detect changes in the relative abundance of the alternative transcripts (isoforms) of genes. This is called Differential Transcript Usage (DTU).
Detecting DTU is supplementary to the quantification of transcripts by tools like Salmon, Sailfish and Kallisto and the detection of Differential Transcript Expression (DTE) by tools such as Sleuth.
I particularly like the figure in the tutorial showing the differences among DTU, DTE and DEG. The paper transcript-level estimates improve gene-level inferences above also talks about the differences:
! MATS is a computational tool to detect differential alternative splicing events from RNA-Seq data. The statistical model of MATS calculates the P-value and false discovery rate that the difference in the isoform ratio of a gene between two conditions exceeds a given user-defined threshold.
A biostar post: Do not feed rounded estimates of gene counts from kallisto into DESeq2 (please make sure you read through all the comments, and now there is a suggested workflow for feeding rounded estimates of gene counts to DESeq etc)
There is some confusion in the answers to this question that hopefully I can clarify with the three comments below:
- kallisto produces estimates of transcript level counts, and therefore to obtain an estimate of the number of reads from a gene the correct thing to do is to sum the estimated counts from the constituent transcripts of that gene. Of note in the language above is the word "estimate", which is necessary because in many cases reads cannot be mapped uniquesly to genes. However insofar as obtaining a good estimate, the approach of kallisto (and before it Cufflinks, RSEM, eXpress and other "transcript level quantification tools") is superior to naïve "counting" approaches for estimating the number of reads originating from a gene. This point has been argued in many papers; among my own papers it is most clearly explained and demonstrated in Trapnell et al. 2013.
- Although estimated counts for a gene can be obtained by summing the estimated counts of the constituent transcripts from tools such as kallisto, and the resulting numbers can be rounded to produce integers that are of the correct format for tools such as DESeq, the numbers produced by such an approach do not satisfy the distributional assumptions made in DESeq and related tools. For example, in DESeq2, counts are modeled "as following a negative binomial distribution". This assumption is not valid when summing estimated counts of transcripts to obtain gene level counts, hence the justified concern of Michael Love that plugging in sums of estimated transcript counts could be problematic for DESeq2. In fact, even the estimated transcript counts themselves are not negative binomial distributed, and therefore also those are not appropriate for plugging into DESeq2. His concern is equally valid with many other "count based" differential expression tools.
- Fortunately there is a solution for performing valid statistical testing of differential abundance of individual transcripts, namely the method implemented in sleuth. The approach is described here. To test for differential abundance of genes, one must first address the question of what that means. E.g. is a gene differential if at least one isoform is? or if all the isoforms are? The tests of sleuth are performed at the granularity of transcripts, allowing for downstream analysis that can capture the varied questions that might make biological sense in specific contexts.
In summary, please do not plug in rounded estimates of gene counts from kallisto into DESeq2 and other tools. While it is technically possible, it is not statistically advisable. Instead, you should use tools that make valid distributional assumptions about the estimates.
However, Charlotte Soneson, Mike Love and Mark Robinson showed in a f1000 paper: Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences that rounded values from transcript level can be fed into DESeq2 etc for gene-level differential expression, and it is valid and preferable in many ways.
Thanks Rob Patro for pointing it out!
Kallisto and caching everything we possibly can.Isolator has a particular focus on producing stable, consistent estimates. Maximum likelihood approaches produce unstable point estimates: small changes in the data can result in drastically different results, conflating downstream analysis like clustering or PCA. Isolator produces estimates that are in general, simultaneously more stable and more accurate other methods
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.
hg38 and gencode v21 as annotation.Considerable differences are found between the methods in terms of the number and characteristics of the genes that are called differentially expressed. Pre-filtering of lowly expressed genes can have important effects on the results, particularly for some of the methods originally developed for analysis of bulk RNA-seq data. Generally, however, methods developed for bulk RNA-seq analysis do not perform notably worse than those developed specifically for scRNA-seq.
325 followers · starred May 2018
283 followers · starred Nov 2019