Hi everyone! I'm working on a gene expression file from an RNA seq experiment and try to compare between two groups of samples. I'm pretty new to this world of gene expressions and need some help. I have a strong reason to believe that there should be a biological difference between the two groups (from the biological aspect). When I do different gene expression I get many differentially expressed genes and when I cluster the samples according to the genses, I don't see two clear separation between the groups. Is there a way to filter the genes I get from my differential expression algorithm so I can see a better clustering effect? I am interested in selecting genes from the list of differentially expressed genes that will separate the two groups in the best way, however, there are so many genes that are differentially expressed that I don't know how to effectively do it. Also, I have been trying two different differential gene expression methods (ttest and limma) and the genes I get from each limma don't appear in the gene list I get from the ttest. Which one should I use? Thanks a lot!
If you already know the two groups and want to find the genes that discriminate between them, try linear discriminant analysis or any other suitable supervised method.
I gave an upvote for this as it is quite a powerful method, i.e., linear discriminant analysis (LDA).
Thanks! have you used LDA on a high throughput data? I was looking for an R code to run it, but could not find one suitable for large scale data.
It depends on what you call large scale data. In R, one typically uses the lda() function in the MASS package.
Consider that you probably don't need to use all your data to build the model. Presumably a subset is enough to use as training set.
If for some reason the lda() function can't handle your data, you could try to do the computation yourself. The most intensive operation is the eigenvector-eigenvalue computation and you could try using solvers for large eigenvalue problems (e.g. the RSpectra package)
Sorry about the confusion. By large scale data I mean the gene expression matrix. It has more than 30,000 genes. Will lda() handle that large amount of features? Also, should I give lda() the raw count as an input? I guess not, because of the difference in library depth and variance between the samples, so should I give it the normalized expression count generated by Deseq2? Thanks
Did you use the LDA on the log transformed expression matrix?
Log in to answer this question.
How are you clustering the genes? Have you tried any sort of pathway or functional enrichment analyses (GSEA, DAVID, etc) of your differentially expressed genes? What is the ttest method, literally just doing t-tests between your two groups?
Personally, I've had good success with
limmaand was able to make biologically interesting (and reasonable) conclusions from its output.How did you normalize the samples in the ttest method? I usually use DESeq2 with great results. You can use PCA to see if the samples can be clustered.
Thanks! I work with DESeq2 too. Do you run the pca on the rlog output of DESeq2?
Thank you all very much! Interestingly, when I do pca with the DE genes I can see a pretty nice separation between my two groups. For some reason though, in the heatmap the samples don't cluster as nicely. I used Euclidean and weighted for building the heatmap. My data is actually genetically heterogeneous. I have few driving translocation that each have a different gene expression signature and I am trying to find a difference between two groups that each group contains samples with different types of translocations. That is, I am trying to find the difference between groups where there is a lot of genetic variability within each group to begin with.
I initially wrote my response as a comment but was asked to write it as a reply to your comment. Sorry for the confusion. When I look at my data using PCA the variance explained by PC1 is 14% and 10% by PC2. I know the samples are biologically different regardless of my grouping as some of them have known mutations while other don't (it may not directly affect the gene expression, but previous studies have shown that within these samples are biologically different from each other. Also, each sample is from different patient, so there is another level of biological variability). Looking at these samples I divided them into two groups that I think should have distinct gene expression signature and colored the pca plot by the these two groups, where it showed that they do separate. The PCA plot was obtained using the expression matrix for all of the samples (i.e. the samples in the two groups) and more than 30,000 genes. I then did differential analysis between the two groups using Deseq2 and used the genes I got to draw a heatmap, but in the heatmap the separation was not that obvious.
Thank you so much for your help! They don't cluster so nicely in the dendogram. I'll try to describe in more detail what I did: I created a Deseq2 object using my raw count matrix and metadata. I then applied variance stabilization transformation on the data and plot PCA. I see two distinct groups in my plot. The groups on the pca correlate with the groups in the design matrix. Here is the code I was using:
dds <-DESeqDataSetFromMatrix(countData = ep,colData = cp,design = ~Risk)
dds <- estimateSizeFactors(dds)
vsd <- varianceStabilizingTransformation(dds)
plotPCA(vsd, intgroup="Risk")
Now I ran Deseq on dds and got 150 genes which are differentially expressed between the two groups in my design matrix. I used those genes to create the heatmap and the clustering wasn't as clear.
Is there a way to know which genes separate the two groups I see in the pca plot? If so, I can try giving those genes to the clustering algorithm and see I get a better separation in the heatmap.
Thank you for the code! Is it supposed to be project.pca <- prcomp(assay(vsd))? The column in my vsd are the samples and the row are the genes. If I go by t(assay(vsd)) I got PCA for the samples rather than the genes.
Thanks! I'm now stuck with another related problem on a different dataset. I run Deseq2 on the dataset, I get a handful of DE genes between two groups of interests in my dataset. I run pca on the expression matrix with only the DE genes selected, and I don't see the two groups separating. How can it be that when I run pca on the dataset using the DE genes (with adj p value less than 0.05 and log fold on 2) the two groups don't separate? Did you run into a similar issue in the past? Which visualization method is best to show two groups clearly separated by the DE genes? Thanks for your help again!
Thank you all very much! Interestingly, when I do pca with the DE genes I can see a pretty nice separation between my two groups. For some reason though, in the heatmap the samples don't cluster as nicely. I used Euclidean and weighted for building the heatmap. My data is actually genetically heterogeneous. I have few driving translocation that each have a different gene expression signature and I am trying to find a difference between two groups that each group contains samples with different types of translocations. That is, I am trying to find the difference between groups where there is a lot of genetic variability within each group to begin with.