<p>bioinformatics scientists.
I hope your work goes smoothly well and be safe.</p>
<p>I'm stuck in trouble at the step of clustering.</p>
<p>Briefly, I downloaded fastq files from a published study (<a href="https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861</a>). There are 24 samples that come from 6 different mouse tumors.
Each tumor consists of two biological replicates where one biological replicate includes two technical replicates. So, I ran the command &quot;cellranger count&quot;, and generated 24 &quot;filtered_feature_bc_matrix&quot; files. For example, for the MC38, bio1-tech1, bio1-tech2, bio2-tech1, bio2-tech2 samples exist. FYI, These samples contain tumor cells as well as tumor-infiltrating immune cells.</p>
<p>And then, I carried out the command &quot;cellranger aggr&quot; to combine all matrix files into one matrix file.
As suggested Seurat tutorial (<a href="https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html" rel="nofollow noreferrer">https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html</a>), Standard pre-processing workflow, linear dimensional reduction, determining dimensionality, clustering, and differential expression were done. Below are the commands I put.</p>
<blockquote>
<p>tumors &lt;- NormalizeData(tumors)</p>
<p>tumors &lt;- FindVariableFeatures(tumors, selection.method = &quot;vst&quot;, nfeatures = 2000)</p>
<p>all.genes &lt;- rownames(tumors)</p>
<p>tumors &lt;- ScaleData(tumors, features = all.genes)</p>
<p>tumors &lt;- RunPCA(tumors, features = VariableFeatures(object = tumors))</p>
<p>tumors &lt;- JackStraw(tumors, num.replicate = 100)</p>
<p>tumors &lt;- ScoreJackStraw(tumors, dims = 1:20)</p>
<p>ElbowPlot(tumors)</p>
<p>tumors &lt;- FindNeighbors(tumors, dims = 1:10) # I chose dims from 5, 6, 7, 8, 9, 10 or 20 as a cutoff.</p>
<p>tumors &lt;- FindClusters(tumors, resolution = 1.2)</p>
<p>tumors &lt;- RunUMAP(tumors, dims = 1:10)</p>
<p>VlnPlot(tumors, features = c(&quot;Cd79a&quot;, &quot;Foxp3&quot;, &quot;Cd3e&quot;, &quot;Cd4&quot;, &quot;Zap70&quot;, &quot;Cd8a&quot;, &quot;Cd14&quot;, &quot;Il2ra&quot;))</p>
</blockquote>
<p><a href="https://i.sstatic.net/rhRyA.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/rhRyA.png" alt="enter image description here" /></a></p>
<p>Interestingly, as you can see here, Foxp3, Cd4, Cd8a are found to be overlapped. This phenomenon was the same when I drew vlnplot using different dimensionality from dims=1:5 to 1:10 or 1:20.</p>
<p>Next, I performed this single-cell RNA sequencing analysis using each tumor; mc38, b16, sa1n, ll2, emt6, and ct26. Representatively, MC38, SA1N, EMT6, and CT26 are presented as below.</p>
<p>For mc38, CD4, CD8, and Foxp3 are overlapped.
<a href="https://i.sstatic.net/Uji2V.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Uji2V.png" alt="enter image description here" /></a>
For sa1n, CD4, CD8, and Foxp3 are separated.
<a href="https://i.sstatic.net/kd7cg.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/kd7cg.png" alt="enter image description here" /></a>
For emt6, CD4, CD8, and Foxp3 are overlapped.
<a href="https://i.sstatic.net/N11NV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/N11NV.png" alt="enter image description here" /></a>
For ct26, CD4, CD8, and Foxp3 are overlapped.
<a href="https://i.sstatic.net/TWm3P.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/TWm3P.png" alt="enter image description here" /></a></p>
<p>So, except for SA1N that only separates Foxp3 from CD8 T cells, I found that other tumor types can't separate Fopx3 from CD8 T cells. Do you think it is related to technical issues, such as values of dimensionality or PCA I put? Since I am not an expert about scRNA, I don't know the exact meaning of dims or PCA, though. Is there anyone who can tell me what's going on? And, do you know how to make the right clustering to get exact T cell subpopulations?</p>
<p>I do appreciate your valuable time and thank you for reading this long story.</p>

CC BY-SA 4.0
2020-05-29T02:16:48.393Z
It could have something to do with the parameter you chose for clustering. Usually, people try many sets of parameters. Also, it is somewhat known that CD4,CD8, and Treg can be really hard to be separated.
Phoenix Mu
CC BY-SA 4.0
2020-05-29T23:14:40.027Z
You can check where the CD4, CD8, and Foxp3 clusters are. The short explanation is that the genes that you are interested in, are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs
StupidWolf
CC BY-SA 4.0
2020-05-29T23:15:32.127Z
Then it's a matter of deciding the resolution to find the clusters.. You can try and play around with this.. I would only have time tomorrow to look at this data and write a solution
StupidWolf
CC BY-SA 4.0
2020-05-31T00:48:42.387Z
@StupidWolf, Hello again. I solved this issue by increasing the resolution up to 1.8 or even higher. I also happened to change the value of min.cell or min.features or others, such as nFeature_RNA, nFeature_RNA, but it appeared that these are "not that critical" enough to change the result. Your advice sounds really important, but I don't understand the part "are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs". Are PCs the PCA? Do you mind If I have a basic question about the terms?  What is the difference between PCA and dimensionality?
raiora
CC BY-SA 4.0
2020-05-31T08:52:54.320Z
You can check this for a visual guide to PCs, https://setosa.io/ev/principal-component-analysis/
StupidWolf
CC BY-SA 4.0
2020-05-31T08:54:25.957Z
hi @raiora, glad you manage to find the cluster ! sorry i haven't found time yet to look at it.. Yes the PCs are the results from PCA. You have a lot of genes, and you can do what is called dimension reduction as you have pointed to. PCA is a method of dimension reduction which finds the linear combination of genes that explains the most variance in your data. The data is projected onto so called principal components which is what I referred to as PCs. This is what you are looking at with the elbow plot, how much variance each PCs explain
StupidWolf
CC BY-SA 4.0
2020-05-31T08:58:33.290Z
So from the first step, ```RunPCA(tumors, features = VariableFeatures(object = tumors))```, you can check whether foxp3 etc is in VariableFeatures(object = tumors)), because this will decide whether the PCs capture the variation in t-reg cells. By default i think 2000 genes are used? Then you can look under your seurat object, tumours@reductions$pca@feature.loadings to see what is the loading for foxp3 etc
StupidWolf
<scrnaseq><seurat><clustering>
Clustering issue about overlapping CD4, CD8, and Foxp3 genes in mouse tumor models
# Clustering issue about overlapping CD4, CD8, and Foxp3 genes in mouse tumor models

(Asked by: username_0 on 2020-05-29)

bioinformatics scientists.
I hope your work goes smoothly well and be safe.

I'm stuck in trouble at the step of clustering.

Briefly, I downloaded fastq files from a published study (<https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861>). There are 24 samples that come from 6 different mouse tumors.
Each tumor consists of two biological replicates where one biological replicate includes two technical replicates. So, I ran the command "cellranger count", and generated 24 "filtered\_feature\_bc\_matrix" files. For example, for the MC38, bio1-tech1, bio1-tech2, bio2-tech1, bio2-tech2 samples exist. FYI, These samples contain tumor cells as well as tumor-infiltrating immune cells.

And then, I carried out the command "cellranger aggr" to combine all matrix files into one matrix file.
As suggested Seurat tutorial (<https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html>), Standard pre-processing workflow, linear dimensional reduction, determining dimensionality, clustering, and differential expression were done. Below are the commands I put.

> tumors <- NormalizeData(tumors)
>
> tumors <- FindVariableFeatures(tumors, selection.method = "vst", nfeatures = 2000)
>
> all.genes <- rownames(tumors)
>
> tumors <- ScaleData(tumors, features = all.genes)
>
> tumors <- RunPCA(tumors, features = VariableFeatures(object = tumors))
>
> tumors <- JackStraw(tumors, num.replicate = 100)
>
> tumors <- ScoreJackStraw(tumors, dims = 1:20)
>
> ElbowPlot(tumors)
>
> tumors <- FindNeighbors(tumors, dims = 1:10) # I chose dims from 5, 6, 7, 8, 9, 10 or 20 as a cutoff.
>
> tumors <- FindClusters(tumors, resolution = 1.2)
>
> tumors <- RunUMAP(tumors, dims = 1:10)
>
> VlnPlot(tumors, features = c("Cd79a", "Foxp3", "Cd3e", "Cd4", "Zap70", "Cd8a", "Cd14", "Il2ra"))

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

Interestingly, as you can see here, Foxp3, Cd4, Cd8a are found to be overlapped. This phenomenon was the same when I drew vlnplot using different dimensionality from dims=1:5 to 1:10 or 1:20.

Next, I performed this single-cell RNA sequencing analysis using each tumor; mc38, b16, sa1n, ll2, emt6, and ct26. Representatively, MC38, SA1N, EMT6, and CT26 are presented as below.

For mc38, CD4, CD8, and Foxp3 are overlapped.
[![enter image description here](https://i.sstatic.net/Uji2V.png)](https://i.sstatic.net/Uji2V.png)
For sa1n, CD4, CD8, and Foxp3 are separated.
[![enter image description here](https://i.sstatic.net/kd7cg.png)](https://i.sstatic.net/kd7cg.png)
For emt6, CD4, CD8, and Foxp3 are overlapped.
[![enter image description here](https://i.sstatic.net/N11NV.png)](https://i.sstatic.net/N11NV.png)
For ct26, CD4, CD8, and Foxp3 are overlapped.
[![enter image description here](https://i.sstatic.net/TWm3P.png)](https://i.sstatic.net/TWm3P.png)

So, except for SA1N that only separates Foxp3 from CD8 T cells, I found that other tumor types can't separate Fopx3 from CD8 T cells. Do you think it is related to technical issues, such as values of dimensionality or PCA I put? Since I am not an expert about scRNA, I don't know the exact meaning of dims or PCA, though. Is there anyone who can tell me what's going on? And, do you know how to make the right clustering to get exact T cell subpopulations?

I do appreciate your valuable time and thank you for reading this long story.

**Tags:** scrnaseq, seurat, clustering
> **Comments:**
> **username_1**: It could have something to do with the parameter you chose for clustering. Usually, people try many sets of parameters. Also, it is somewhat known that CD4,CD8, and Treg can be really hard to be separated.  
> **username_2**: You can check where the CD4, CD8, and Foxp3 clusters are. The short explanation is that the genes that you are interested in, are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs  
> **username_2**: Then it's a matter of deciding the resolution to find the clusters.. You can try and play around with this.. I would only have time tomorrow to look at this data and write a solution  
> **username_0**: @StupidWolf, Hello again. I solved this issue by increasing the resolution up to 1.8 or even higher. I also happened to change the value of min.cell or min.features or others, such as nFeature_RNA, nFeature_RNA, but it appeared that these are "not that critical" enough to change the result. Your advice sounds really important, but I don't understand the part "are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs". Are PCs the PCA? Do you mind If I have a basic question about the terms?  What is the difference between PCA and dimensionality?  
> **username_2**: You can check this for a visual guide to PCs, https://setosa.io/ev/principal-component-analysis/  
> **username_2**: hi @raiora, glad you manage to find the cluster ! sorry i haven't found time yet to look at it.. Yes the PCs are the results from PCA. You have a lot of genes, and you can do what is called dimension reduction as you have pointed to. PCA is a method of dimension reduction which finds the linear combination of genes that explains the most variance in your data. The data is projected onto so called principal components which is what I referred to as PCs. This is what you are looking at with the elbow plot, how much variance each PCs explain  
> **username_2**: So from the first step, ```RunPCA(tumors, features = VariableFeatures(object = tumors))```, you can check whether foxp3 etc is in VariableFeatures(object = tumors)), because this will decide whether the PCs capture the variation in t-reg cells. By default i think 2000 genes are used? Then you can look under your seurat object, tumours@reductions$pca@feature.loadings to see what is the loading for foxp3 etc  

---

## 0 answers

