A fairly standard protocol could be:

1. use something like tophat (if you map on the genome sequence) or bowtie2 (if you map directly on the transcripts) to align the reads;

2. quantify the abundance of each transcript e.g. you can use cufflinks in conjunction with the genome sequence where you mapped the reads and the gff file containing the coordinates of all genes. Alternatively you can use coverageBed or similar to simply count how many reads are overlapped with the features in your gff.

3. this will give you one file for library (or all libraries in the case of cufflinks) so you have to concatenate them taking into account the gene identifiers. In case you use cufflinks, check carefully the command line you use, as it allows to apply multiple corrections (for instance with eukaryotes it is important to correct for compositional biases and multi-mapping something like -u -b parameters in cufflinks, parameters that, however consumes a lot of memory, so, if you get a segmentation fault, try at first to remove those two parameters). You can follow the pipeline suggested on the documentation of cufflinks.

4. let's assume you were able to quantify the abundance of your transcripts. cufflinks will give you both the fpkm information in one file and the (normalized for the above biases) counts. Depending on your taste, you can use both. Maybe check some theory about the two approaches (counts-rpkm vs fpkm). You have to parse the file with the abundances such that you can load them e.g. in R or Matlab.

5. use a clustering algorithm after standardization of the expression level. Better first to remove genes that are very lowly expressed (or even too much expressed, which may indicate e.g. multiple copies of the gene which indicates that establishing the right expression level for those transcripts is much more difficult than for the rest of the genes). I usually standardize by row assuming that the rows are genes and the columns are different libraries. Please note that standardization is important and in general it is better to apply it on log transformed counts or fpkm, but not always. It is a matter of trial and error. To start you can simply use an agglomerative clustering (there are many both in R and matlab). Please also notice I am assuming that at this point you put together the replicates (which is however done automatically by cufflinks if I remember well). Checking for this is also very important because it also allows to identify genes with very large variability which is another indication that something went wrong with that gene. better to remove them. IF you want to reduce even more your matrix, you can also think about removing genes with no or very little variability. At the extreme of perfectly constant gene expression level, you may end up with errors depending on the distance you use for the clustering e.g. correlation does not accept constant data and degrades consistently when the variability is very small as you are basically calculating the correlation in noise (if the small variations are indeed noise).

So I would start with a distance matrix to be input to the clustering algorithm calculated using 1-correlation or the euclidean distance (which are basically the same distance with standardized data). The algorithm in its simplest form will assign a cluster number to each gene, such that genes with the same number belong to the same cluster (i.e. they are most similar concerning their expression level).

6. It is very important to use a package like clValid (R) to test the ensuing clustering, as there are objective measures of its quality (connectivity, Silhouette Width and Dunn Index). If you have some information about what could be a biological validation of your clusters e.g. if you think that the clusters should be functionally homogeneous, that you can also integrate the functional annotation of genes and use some external validation measure again available in clValid. In this case the algorithm will run the clustering many times, each time calculating one of the external measures, to give you the best output possible (following the criterion you asked for).

You can also take a look at the following packages in R:

clValid

WGCNA

Petal

(simply look for "gene coexpression networks R-cran" in google)

Matlab has a few very interesting packages (in general also available in R), like the one developed by Faith et al. in 2008Article Many Microbe Microarrays Database: Uniformly normalized Affy...

But you can apply any clustering algorithm you are familiar with, maybe you'll need to adapt it, for instance the MCL (by Stijn van Dongen), SOM ecc ecc.

I think that the most important thing is to get a good quantification (so applying the right normalization, bias corrections) than stick with the most accurate clustering algorithm.

Hope I was clear enough.

You can use WGCNA package in R for extracting the co-expressed modules or clusters. Also, you have to use limma package for normalizing the data.

If you are starting with RNA-seq data, one standard workflow used in our group (and used by The Cancer Genome Atlas Project) is to map RNA-seq reads using MapSplice, then generate the gene expression matrix using RSEM followed by Upper Quartile Normalization (see: https://wiki.nci.nih.gov/display/TCGA/RNASeq+Version+2)

Another approach would be to use STAR followed by SALMON

Hope this helps,

Benjamin

A fairly standard protocol could be:

1. use something like tophat (if you map on the genome sequence) or bowtie2 (if you map directly on the transcripts) to align the reads;

2. quantify the abundance of each transcript e.g. you can use cufflinks in conjunction with the genome sequence where you mapped the reads and the gff file containing the coordinates of all genes. Alternatively you can use coverageBed or similar to simply count how many reads are overlapped with the features in your gff.

3. this will give you one file for library (or all libraries in the case of cufflinks) so you have to concatenate them taking into account the gene identifiers. In case you use cufflinks, check carefully the command line you use, as it allows to apply multiple corrections (for instance with eukaryotes it is important to correct for compositional biases and multi-mapping something like -u -b parameters in cufflinks, parameters that, however consumes a lot of memory, so, if you get a segmentation fault, try at first to remove those two parameters). You can follow the pipeline suggested on the documentation of cufflinks.

4. let's assume you were able to quantify the abundance of your transcripts. cufflinks will give you both the fpkm information in one file and the (normalized for the above biases) counts. Depending on your taste, you can use both. Maybe check some theory about the two approaches (counts-rpkm vs fpkm). You have to parse the file with the abundances such that you can load them e.g. in R or Matlab.

5. use a clustering algorithm after standardization of the expression level. Better first to remove genes that are very lowly expressed (or even too much expressed, which may indicate e.g. multiple copies of the gene which indicates that establishing the right expression level for those transcripts is much more difficult than for the rest of the genes). I usually standardize by row assuming that the rows are genes and the columns are different libraries. Please note that standardization is important and in general it is better to apply it on log transformed counts or fpkm, but not always. It is a matter of trial and error. To start you can simply use an agglomerative clustering (there are many both in R and matlab). Please also notice I am assuming that at this point you put together the replicates (which is however done automatically by cufflinks if I remember well). Checking for this is also very important because it also allows to identify genes with very large variability which is another indication that something went wrong with that gene. better to remove them. IF you want to reduce even more your matrix, you can also think about removing genes with no or very little variability. At the extreme of perfectly constant gene expression level, you may end up with errors depending on the distance you use for the clustering e.g. correlation does not accept constant data and degrades consistently when the variability is very small as you are basically calculating the correlation in noise (if the small variations are indeed noise).

So I would start with a distance matrix to be input to the clustering algorithm calculated using 1-correlation or the euclidean distance (which are basically the same distance with standardized data). The algorithm in its simplest form will assign a cluster number to each gene, such that genes with the same number belong to the same cluster (i.e. they are most similar concerning their expression level).

6. It is very important to use a package like clValid (R) to test the ensuing clustering, as there are objective measures of its quality (connectivity, Silhouette Width and Dunn Index). If you have some information about what could be a biological validation of your clusters e.g. if you think that the clusters should be functionally homogeneous, that you can also integrate the functional annotation of genes and use some external validation measure again available in clValid. In this case the algorithm will run the clustering many times, each time calculating one of the external measures, to give you the best output possible (following the criterion you asked for).

You can also take a look at the following packages in R:

clValid

WGCNA

Petal

(simply look for "gene coexpression networks R-cran" in google)

Matlab has a few very interesting packages (in general also available in R), like the one developed by Faith et al. in 2008Article Many Microbe Microarrays Database: Uniformly normalized Affy...

But you can apply any clustering algorithm you are familiar with, maybe you'll need to adapt it, for instance the MCL (by Stijn van Dongen), SOM ecc ecc.

I think that the most important thing is to get a good quantification (so applying the right normalization, bias corrections) than stick with the most accurate clustering algorithm.

Hope I was clear enough.