Hi Mohamed - If you are familiar with R, it is all very easy: use DESeq2 package to produce normalized log-transformed data (using functions rlog or getVarianceStabilisedData), and then use packages pheatmap for clustering and vegan for principal coordinate analysis. See attached R script, should be easily modifiable for any dataset.

I think I used HC for RNA-seq 

Hi Mohamed - If you are familiar with R, it is all very easy: use DESeq2 package to produce normalized log-transformed data (using functions rlog or getVarianceStabilisedData), and then use packages pheatmap for clustering and vegan for principal coordinate analysis. See attached R script, should be easily modifiable for any dataset.

1) it has to be normalized, using counts is not proper

2) normalized data is the same as other readout, like array.

3) you should already normalized your counts. If you have multiple treatment groups, you may need to normalize the data.

Some good answers already, the only thing I would add is that clustering over all genes may fail due to the size of the problem, so you may want to limit the analysis to a subset of genes, for instance the ones which are differentially expressed, e.g. as assessed by DESeq2.

To add to Jinsong and Jonathan's answers above: depending on the read depth of your data, you might encounter genes that do not have reads associated with them. This poses a difficulty, because you are likely to find lowly expressed genes with no associated reads in one condition vs detectable, but still low detection in the other condition. This results in -Inf or +Inf fold changes, which mess up your heatmaps and clustering analysis. The problem is that the 'no expression -> low expression' situation might be quite relevant for your question, so just throwing out the genes with an average of zero reads in one condition is probably not what you want. DeSeq provides a model to statistically assess whether the changes are relevant, but I think the clustering problem remains.

Thank you for all your useful answers,

I'll try DESeq2.

I have an other question, how to choose the good fitType in DESeq2 ? 

Dear Mohamed,

I would be more in favor of performing EdgeR analysis on your count data if aim is to compare the groups for differentially expressed region. Many methods has been tried  for different count data including RNA-seq, NGS, Chip-Seq and found that EdgeR as more appropriate, although its more personal choice of bioinformaticians/analyst.

Actually in DESeq, sometime you found massive fold change but not significant of particular region, mainly because DESeq works through the logic of independent filtering and then you need to filter the data based on quantile (this can be done using DESeq vignette )

****************************************************************************************

For normalization of your count data using EdgeR, following scripts would be helpful:

## edgeR run script

library(edgeR)

setwd("C:/Documents and Settings/user/My Documents/")

# define input Expression file name to read

input

Dear Tushar - Edger is more appropriate, for counts data, really? Can you please provide the reference for that benchmarking study you mentioned? As far as I am aware, EdgeR is the older-fashion package orignially designed not for counts but for microarrays (color intensity values), and DESeq2 is the latest development for counts-based data. Please tell me I am wrong...

EdgeR was developed for microarray data years ago, but similar statistics are used in DESeq2 as well. I don't think there is much difference between EdgeR and DESeq2 (DESeq 1 is a bit of less sensitive though) in terms of RNAseq data analysis. You may see a difference of several genes, which are usually ambiguous calls anyway, when you apply the both methods to a same dataset. If you want to be comprehensive, you can use differentially expressed genes called from both packages as well.

Dear Prof. Mikhail,

You are absolutely right regarding the points of EdgeR and DeSEq2, but actually that there is not much difference in EdgeR and DESeq statistics of negative binomial dfistribution. (here I meant by DeSeq1 more based on personal experience)

ref.: http://www-huber.embl.de/pub/pdf/nprot.2013.099.pdf.

http://biorxiv.org/content/early/2014/05/28/005611

http://bib.oxfordjournals.org/content/16/1/59.long

Usually, I performed analysis by both methods to find the deferentially expressed genes.

Dear Mikhail,

you are wrong :-) -- edgeR has been developed for count data, please see http://www.ncbi.nlm.nih.gov/pubmed/19910308

It is limma, not edgeR, that has been developed for microarrays originally and can also be used to analyze RNA-Seq data via the voom() transformation. 

One benchmarking study is http://www.biomedcentral.com/1471-2105/14/91 

although I'm not sure that's the one mentioned before...

Cheers

DESeq and EdgeR differ in their normalization techniques: edgeR uses the trimmed mean of M values whereas DESeq uses a relative log expression approach by creating a virtual library that every sample is compared against. They also differ in the choice they make to estimate the dispersion. EdgeR uses dispersion estimates toward a trended mean according to the dispersion mean relationship. DEseq takes the maximum of the individual dispersion estimates and the dispersion mean trend. In practice, this means that EdgeR is more sensitive to outliers. However, looking at the literature, there might be only a few differences in the output of the two packages.

If your aim is to be really careful with your outputs, run both and keep only the genes showing up in both analyses...

 Oh yes, this backs up my feeling that DESeq, DESeq2 and EdgeR are all acceptable tools. Indeed I was wrong about EdgeR being not a counts analysis method - confused with limma.

I am a big fan of DESeq2 personally, with only one exception - the "independence filtering" procedure (discarding low-abundant genes to boost FDR-passing numbers). Feels wrong to me for some reason. I wrote an alternative simulation-based procedure that seems to be no less powerful than independence filtering but does not require actually discarding genes from the analysis. It is the R package empiricalFDR.DESeq2 - please check it out (you can get it from CRAN).

I agree with Mikhail V Matz. For my RNA-Seq data DESeq2 has provided the most meaningful results based on downstream analyses.  DESeq and EdgeR are comparable.  Once you have the transformed counts (i.e. rlog) then you can use the tool of your choice to cluster the data.  I also agree that I am not 100% comfortable with the independence filtering procedure, but otherwise the package has been very useful. 

For clustering data you may want to try mbcluster.seq, which has been developed specifically for count data. 

http://cran.r-project.org/web/packages/MBCluster.Seq/index.html

Anna

I thought I must mention Trinity. Basically it normalize the data and use both (users choice) to generate the RPKM/FPKM value. Nicely written perl scripts are included in this to call for differentially expressed genes. Please have a look.

Quite interesting discussion over clustering, clustering methods and finding differentially expressed genes and the tools one could use to find DE genes and for clustering. I would also like to add a few comments.

Although, this is already mentioned I would like to strengthen it by stating that RPKM/FPKM approach to determine DE genes is useful, when one compares expression levels of mRNA across a single library.

When it is needed to find DE genes between two libraries (two conditions) it does not really matter, whether one uses count or RPKM/FPKM, simply because the mRNA length does not play a role anymore, as the comparison is between expression levels of the same gene between two libraries.

Selecting a clustering method that would help us to extract clusters has been a problem for a very long time. Many in the field tried many different approaches to fine tune clustering approach to find clusters representing whole or a part of molecular networks. In recent years clustering has become quite important in the context of identifying novel players of molecular networks and novel molecular networks. Principal component analysis provides a way to reduce multidimensional data sets into two or three dimensions and to see whether the dimensionality reduction help us to identify gene sets that behave the same manner.

The following could provide a better understanding of clustering, data preparation and methods to use

http://www.biostat.jhsph.edu/GenomeCAFE/ExpressionistSeminarSlides/ClusterClassTopTenCarlo.ppt

http://zsi.tech.us.edu.pl/~nowak/bien/MODA3.pdf

http://www.biomedcentral.com/1471-2105/15/S2/S2

There are many clustering algorithms implemented in R and published in Bioconductor e.g.

Mfuzz (it also has a gui to make your life easy)

In recent years machine learning was getting popular in identifying clusters by making decisions using classifiers. One could make an attempt to use a machine learning algorithm available online with some tweaking. One could also use training sets (a positive and negative sets) to train the algorithm.

This has been going on for maybe too long.. but if anyone is still following: I was assuming that the original question by Mohamed was about how to cluster samples. To cluster genes into correlated modules, this is another issue. For gene clustering, we are big fans of weighted gene correlation network analysis (WGCNA) - it is really great, but you must have 15+ samples (it does not make much sense to try correlation-based clustering with fewer). See excellent WGCNA tutorials here: http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/

Correct me if I am wrong. Isn't it that clustering samples and cluster genes into correlated modules are two neighboring mile stones on the same road (or journey)? I have been using WGCNA as well and yes it is great a tool. Thanks Mikhail for bringing it up here, in fact it was in my mind when I was writing my previous post, but forgot to mention.

Hi Mohamed,

DeSeq, EdgeR or EasyRNAseq start form raw-count data. If you have all your tuxedo pipeline setup (tophat, cufflinks, cuffdiff), the easiest might be to use cumeRbund to explore your clusters.

If instead you have an .xls or tabular file with raw counts you can use one of the methods above. Services like InSilicoDB, Maverix, or StationX allow you to import and analyze your data easily. You can do this easily and free in InSilico DB (https://insilicodb.com).

To clarify the difference between Raw counts, FPKMs and RPKMs I recommend Lior Patcher's talk: https://www.youtube.com/watch?v=5NiFibnbE8o

I propose you NMF to perform cluster, is a robust method to perform cluster and obtain features (genes) that caracterize each cluster, them, you can perform an enrichmwnt analysis (for example with EnrichR web tool) on the lists of features (one per cluster). If you get too few genes for conventional enrichmwnt analysis. You can sort the list and perform GSEA to get functional pathways over and under regulated in each cluster

Where is the data?