<p>Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.</p>

CC BY-SA 3.0
2018-04-18T13:41:03.480Z
Yes, you are right. If some thing went wrong what all should I check anything like mapped reads percentage or anything else?
stack_learner
CC BY-SA 3.0
2018-04-18T13:42:52.657Z
Mapped percentages, size factors, and perhaps sparsity (percentage of genes with no counts). Basically you want to ensure that they're all similar. Also, try redoing the counts with featureCounts. Then you can eliminate an oddity in the counting as an issue.
Devon Ryan
CC BY-SA 3.0
2018-04-18T18:28:45.197Z
Which is the mapping percentage in the above comments? It gives overall alignment rate.
stack_learner
CC BY-SA 3.0
2018-04-18T18:45:07.507Z
The alignment rate is the mapping percentage.
Devon Ryan
CC BY-SA 3.0
2018-04-19T07:55:05.093Z
Counts with featurecounts should I use bam or sorted.bam?
stack_learner
CC BY-SA 3.0
2018-04-19T07:55:36.877Z
It doesn't matter which you use.
Devon Ryan
CC BY-SA 3.0
2018-04-19T09:24:09.820Z
I initially checked on one bam file. Geneid, Chr, Start, End, Strand, Length, /path/to/A.sorted.bam. These are the columns in the featureCounts output. Is the last column read counts? Do you think this is fine? featureCounts -a annotation.gtf -F GTF -p -s 2 -T 8 -P -o countMatrix_all.txt /pathto/*.sorted.bam
stack_learner
CC BY-SA 3.0
2018-04-19T09:42:15.177Z
Yes, the last column is the counts (I wish there were an option to simplify that output to just geneid and count). You probably don't need the `-P` option, not that it's likely hurting anything. BTW, that could be its own post here.
Devon Ryan
CC BY-SA 3.0
2018-04-19T11:36:27.040Z
One last question. When I did samtools flagstat on one of the bam file i see this. 202671876 + 0 in total (QC-passed reads + QC-failed reads)
128036784 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
179886195 + 0 mapped (88.76% : N/A)
74635092 + 0 paired in sequencing
37317546 + 0 read1
37317546 + 0 read2
45092910 + 0 properly paired (60.42% : N/A)
48583902 + 0 with itself and mate mapped
3265509 + 0 singletons (4.38% : N/A)
621892 + 0 with mate mapped to a different chr
377778 + 0 with mate mapped to a different chr (mapQ>=5) In this what is mapped and properly paired?
stack_learner
CC BY-SA 3.0
2018-04-19T13:20:07.770Z
Please post things like this as stand alone questions (make sure to google them first, this is answered elsewhere).
Devon Ryan
CC BY-SA 3.0
2018-04-18T11:54:43.307Z
2018-04-18T11:54:43.307Z
Devon Ryan
<p>The largest variation differences in total expression will be seen in cell population differences. It's a good thing that your conditions aren't clustering on the first and second principal components when all genes are considered, because it means that cell population is probably not a factor in the variation of expression.</p>

<p>There may be some other principal component that <em>does</em> separate the conditions, but you might need to drill down fairly deep into the data. This will depend on how much the treatment conditions contribute to the total expression variation. We got lucky with our <a href="http://jem.rupress.org/content/early/2016/12/01/jem.20160470#F3" rel="nofollow noreferrer">Th2 RNASeq study</a> in that cell population was mostly on the first and second principal components, and treatment condition appeared on the third principal component.</p>

<p>If you want to force condition clustering, identify the genes that are differentially expressed, and only use those genes for the principal component analysis.</p>

CC BY-SA 3.0
2018-04-19T07:08:12.440Z
I'm not sure how useful doing a PCA on DE genes is. You already know your samples differ by group in those genes, after all.
Devon Ryan
CC BY-SA 3.0
2018-04-18T22:46:06.663Z
2018-04-18T22:46:06.663Z
gringer
<p>I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in <a href="http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual" rel="nofollow noreferrer">Stringtie manual</a>. With the counts data I have used edgeR for Differential analysis.</p>

<p>Before that I have seen MDS plot showing in the following way:</p>

<p><a href="https://i.sstatic.net/peGiD.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/peGiD.png" alt="enter image description here"></a></p>

<p>I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.</p>

<p><a href="https://i.sstatic.net/2Dujc.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/2Dujc.png" alt="enter image description here"></a></p>

<p>Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?</p>

<rna-seq><r><rna><clustering><edger>
What could be the reason for the samples not clustering?
# What could be the reason for the samples not clustering?

(Asked by: username_0 on 2018-04-18)

I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in [Stringtie manual](http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual). With the counts data I have used edgeR for Differential analysis.

Before that I have seen MDS plot showing in the following way:

[![enter image description here](https://i.sstatic.net/peGiD.png)](https://i.sstatic.net/peGiD.png)

I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.

[![enter image description here](https://i.sstatic.net/2Dujc.png)](https://i.sstatic.net/2Dujc.png)

Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?

**Tags:** rna-seq, r, rna, clustering, edger

---

## 2 answers

### Answer by username_1 on 2018-04-18. Score: 2

Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.

> **Comments:**
> **username_0**: Yes, you are right. If some thing went wrong what all should I check anything like mapped reads percentage or anything else?  
> **username_1**: Mapped percentages, size factors, and perhaps sparsity (percentage of genes with no counts). Basically you want to ensure that they're all similar. Also, try redoing the counts with featureCounts. Then you can eliminate an oddity in the counting as an issue.  
> **username_0**: Which is the mapping percentage in the above comments? It gives overall alignment rate.  
> **username_1**: The alignment rate is the mapping percentage.  
> **username_0**: Counts with featurecounts should I use bam or sorted.bam?  
> **username_1**: It doesn't matter which you use.  
> **username_0**: I initially checked on one bam file. Geneid, Chr, Start, End, Strand, Length, /path/to/A.sorted.bam. These are the columns in the featureCounts output. Is the last column read counts? Do you think this is fine? featureCounts -a annotation.gtf -F GTF -p -s 2 -T 8 -P -o countMatrix_all.txt /pathto/*.sorted.bam  
> **username_1**: Yes, the last column is the counts (I wish there were an option to simplify that output to just geneid and count). You probably don't need the `-P` option, not that it's likely hurting anything. BTW, that could be its own post here.  
> **username_0**: One last question. When I did samtools flagstat on one of the bam file i see this. 202671876 + 0 in total (QC-passed reads + QC-failed reads) 128036784 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 179886195 + 0 mapped (88.76% : N/A) 74635092 + 0 paired in sequencing 37317546 + 0 read1 37317546 + 0 read2 45092910 + 0 properly paired (60.42% : N/A) 48583902 + 0 with itself and mate mapped 3265509 + 0 singletons (4.38% : N/A) 621892 + 0 with mate mapped to a different chr 377778 + 0 with mate mapped to a different chr (mapQ>=5) In this what is mapped and properly paired?  
> **username_1**: Please post things like this as stand alone questions (make sure to google them first, this is answered elsewhere).  

---

### Answer by username_2 on 2018-04-18. Score: 0

The largest variation differences in total expression will be seen in cell population differences. It's a good thing that your conditions aren't clustering on the first and second principal components when all genes are considered, because it means that cell population is probably not a factor in the variation of expression.

There may be some other principal component that *does* separate the conditions, but you might need to drill down fairly deep into the data. This will depend on how much the treatment conditions contribute to the total expression variation. We got lucky with our [Th2 RNASeq study](http://jem.rupress.org/content/early/2016/12/01/jem.20160470#F3) in that cell population was mostly on the first and second principal components, and treatment condition appeared on the third principal component.

If you want to force condition clustering, identify the genes that are differentially expressed, and only use those genes for the principal component analysis.

> **Comments:**
> **username_1**: I'm not sure how useful doing a PCA on DE genes is. You already know your samples differ by group in those genes, after all.  

---

