SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
clustering RNA seq data gfmgfm Bioinformatics 6 05-22-2013 01:55 AM
RNA-Seq: Full-length transcriptome assembly from RNA-Seq data without a reference gen Newsbot! Literature Watch 7 10-26-2011 05:37 AM
Clustering RPKM data anjulka Bioinformatics 1 10-17-2011 09:27 PM
clustering RNASeq data PFS Bioinformatics 0 06-21-2011 09:14 AM
Clustering gene expression pattern of next generation seq data zhlyang Bioinformatics 0 12-02-2009 02:03 PM

Reply
 
Thread Tools
Old 04-06-2012, 04:28 AM   #1
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default Hierarchical clustering of RNA-seq data

I'm clearly out of my comfort zone here but it needs to be done. So I've taken a cursory look into clustering some RNA-seq data and I am more confused then when I started. Does anyone know of any good reviews that could help my mathematically limited knowledge-base understand some of the complexities of hierarchical clustering, especially with respect to RNA-seq data.

I have taken a look at a couple things I have found.
1. DESeq has some clustering abilities but seems pretty limited.
2. MBCluster.Seq ( http://cran.r-project.org/web/packag...luster.Seq.pdf ) looks interesting but seems to lack human friendly documentation.
3. hopach ( http://bioc.ism.ac.jp/2.6/bioc/html/hopach.html ) looks interesting but there is no documentation on how to feed it RNA-seq data (or any data for that matter).

Data preprocessing steps I've considered but may or may not be correct.
1. Seems best to select genes with high variance.
2. Perhaps values should be log2 transformed??
3. Filtering genes with low read counts seems like a good idea.
4. Perhaps depending on what you want to see, divide by lowest expressed condition for each gene so you are looking a fold change.
5. Are there tools that take into account variance within a condition???

It all looks pretty complicated so any insights would be greatly appreciated.

Thanks!!
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 04-06-2012, 06:24 AM   #2
themerlin
Member
 
Location: Flagstaff, AZ

Join Date: Feb 2010
Posts: 49
Default

I would try out the multiple experiment viewer.

http://www.tm4.org/mev/
themerlin is offline   Reply With Quote
Old 04-06-2012, 06:39 AM   #3
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default

That looks like a great resource. Thanks!!!!!
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 04-07-2012, 02:51 AM   #4
Wolfgang Huber
Senior Member
 
Location: Heidelberg, Germany

Join Date: Aug 2009
Posts: 100
Default

Dear ETHANol,

regarding the question about log-transformation, that is not a good idea since then you will create missing values in those cases where a gene had 0 counts in some of the samples. Genes that have 0 counts in some samples but not in others may be among the most interesting ones.

Now you can of course use log2(n+c), with some suitable choice of c. If one explores the question 'what is the best choice of c?' a bit deeper, one place to end up is variance stabilising transformation. This is explained in the DESeq vignette.

You said "DESeq has some clustering abilities but seems pretty limited". In fact, DESeq is a package for the R language and environment, which probably offers the most wide-ranging spectrum of clustering and machine learning methods on the planet. See the CRAN task view for an overview. DESeq itself does not do clustering, it 'just' preprocesses the data (and offers statistical tests and regression modelling).

Best wishes
Wolfgang
__________________
Wolfgang Huber
EMBL
Wolfgang Huber is offline   Reply With Quote
Old 04-07-2012, 07:18 AM   #5
kopi-o
Senior Member
 
Location: Stockholm, Sweden

Join Date: Feb 2008
Posts: 312
Default

Excellent points by Wolfgang and even more to the point with regard to the original question: the CRAN Cluster Analysis Task View
kopi-o is offline   Reply With Quote
Old 04-07-2012, 10:06 AM   #6
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default

Thanks for the comments.

I have already replaced all the 1's in my data set with 0's. Seems like a reasonable thing to do with RNA-seq data. The log2(n+c) sounds like a better option then just cutting off the low expressed genes at some arbitrary quantile. I'll have to read the DESeq vignette a little more carefully, but when it got to clustering it was pretty much over my head and I don't really see a way to take it further (not saying there isn't - I just don't have the background).

The links look like a jungle and sort of reinforce the idea that this is really difficult. But with anything there has to be a way in.
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 04-07-2012, 10:43 AM   #7
kopi-o
Senior Member
 
Location: Stockholm, Sweden

Join Date: Feb 2008
Posts: 312
Default

If you want a few R functions/packages to get started, you could try "agnes()" or good old "heatmap()".
kopi-o is offline   Reply With Quote
Old 04-07-2012, 10:32 PM   #8
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default

Quote:
Originally Posted by Wolfgang Huber View Post
DESeq itself does not do clustering, it 'just' preprocesses the data (and offers statistical tests and regression modelling).
I completely missed that point when looking at the vignette. So that is really nice as you can take the preprocessed data from DESeq and use it how you wish. Thanks for the heads up.

Code:
library(DESeq)

countsTable <- read.delim("countstable.txt", header=TRUE, stringsAsFactors=TRUE)
rownames( countsTable ) <- countsTable$gene
countsTable <- countsTable[ , -1 ]
conds <- factor(c("pos","pos","pos","pos","neg","neg","neg","neg"))
cds <- newCountDataSet( countsTable, conds )
cds <- estimateSizeFactors( cds )
normcds <- counts( cds, normalized=TRUE )

cds <- estimateDispersions( cds )
res <- nbinomTest( cds, "pso", "neg" )
write.table(res, file="DEG_list.txt", quote=FALSE, sep="\t", row.names=FALSE)



cdsBlind <- estimateDispersions( cds, method="blind" )
vsd <- getVarianceStabilizedData( cdsBlind )

select <- order(res$pval)[1:40]
colors <- colorRampPalette(c("white","darkblue"))(100)
heatmap( vsd[select,], col = colors, scale = "none")
A question on this code, which is taken directly out of the DESeq vigentte.
The object 'cds' is a counts data set in which the conditions are ascribed to each column. If I have several conditions, can I just write this code:
Code:
conds <- factor(c("pos","pos","het","het","het","neg","neg","neg"))
The the top 40 genes are then selected using this line:
Code:
select <- order(res$pval)[1:40]
In the case of an experiment that has several conditions and thus not being able to identify differentially expressed genes by one simple pairwise comparison, would there be some way to select the best genes, say using the genes that show the highest variance?

Thanks again for the help.
__________________
--------------
Ethan

Last edited by ETHANol; 04-07-2012 at 10:34 PM.
ETHANol is offline   Reply With Quote
Old 04-07-2012, 11:21 PM   #9
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 976
Default

Quote:
each column. If I have several conditions, can I just write this code:
Code:
conds <- factor(c("pos","pos","het","het","het","neg","neg","neg"))
Yes, that shoul work. Try it out.

Quote:
In the case of an experiment that has several conditions and thus not being able to identify differentially expressed genes by one simple pairwise comparison, would there be some way to select the best genes, say using the genes that show the highest variance?
Exactly. You could try
Code:
library(genefilter)
select <- order( -rowVars(vsd) )[1:40]
The '-' is to get the highest, not the lowest, variances.

And maybe use a few more than just 40 genes.

Simon
Simon Anders is offline   Reply With Quote
Old 04-09-2012, 05:17 AM   #10
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default

Thanks Simon and everyone else. This has all been a lot of help!!!
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 05-15-2012, 10:44 AM   #11
ETHANol
Senior Member
 
Location: Greece

Join Date: Feb 2010
Posts: 309
Default

I'm back to clustering. I wish there was someone else here more qualified to do this but no such luck.

I'll describe an experiment. I have four conditions which I will call:
Control
Treatment-A
Treatment-B
Treatment-C

I have three replicates for each condition.
What I would like to find out is the genes that are differentially expressed comparing the control to Treatment-A, Treatment-B and Treatment C. These binary comparisons are easy to do with DESeq.

I would also like to know which genes are using the control as a baseline up/down regulated in Treatment-A AND Treatment-B, Treatment-B AND Treatment-C, Treatment-A AND Treatment-C, and all three.

I believe throwing the gene lists from each pairwise comparison into Venn diagrams would produce a lot of false differences.

Is clustering the correct approach to getting these lists of genes? Anybody have any suggestion on how to best do this? Should I use the replicates or just the averages for each condition?
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 08-01-2013, 03:12 PM   #12
vyellapa
Member
 
Location: phoenix

Join Date: Oct 2011
Posts: 58
Default

After the normalizing step

Code:
data <-counts(cds, normalized=TRUE)
My samples are all tumor only.Then, is it appropriate to select genes based on a Standard deviation cutoff or a list of top 'n' genes with highest SD?

Code:
dataLog2=log2(data+1)
for (i in (1:94137)) {if(sd(dataLog2[i,])>3.5 {geneList=union(geneList, rownames (dataLog2)[i]);}}

Last edited by vyellapa; 08-01-2013 at 03:15 PM.
vyellapa is offline   Reply With Quote
Old 02-18-2014, 12:35 PM   #13
vyellapa
Member
 
Location: phoenix

Join Date: Oct 2011
Posts: 58
Default

I'm trying to identify subtypes in my tumor RNA-Seq data through unsupervised hierarchical clustering. For this, I'm curious if counts should normalized (using counts( cds, normalized=TRUE )) before estimateDispersions and Variance Stabilization or should counts not be normalized before these steps.

Thank you,
Teja

Code:
countsTable <- read.delim("countstable.txt", header=TRUE, stringsAsFactors=TRUE) 
rownames( countsTable ) <- countsTable$gene 
countsTable <- countsTable[ , -1 ] 
conds <- factor(c("tumor","tumor","tumor","tumor","tumor","tumor","tumor","tumor")) 
cds <- newCountDataSet( countsTable, conds ) 
cds <- estimateSizeFactors( cds ) 
cds <- counts( cds, normalized=TRUE )  
cdsBlind <- estimateDispersions( cds, method="blind" ) 
vsd <- getVarianceStabilizedData( cdsBlind )  
select <- order(res$pval)[1:250] 
colors <- colorRampPalette(c("white","darkblue"))(100) heatmap( vsd[select,], col = colors, scale = "none")
vyellapa is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 09:50 AM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.