2
$\begingroup$

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

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 For sa1n, CD4, CD8, and Foxp3 are separated. enter image description here For emt6, CD4, CD8, and Foxp3 are overlapped. enter image description here For ct26, CD4, CD8, and Foxp3 are overlapped. enter image description here

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.

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ 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. $\endgroup$
    – Phoenix Mu
    May 29 '20 at 2:16
  • $\begingroup$ 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 $\endgroup$
    – StupidWolf
    May 29 '20 at 23:14
  • 1
    $\begingroup$ You can check this for a visual guide to PCs, setosa.io/ev/principal-component-analysis $\endgroup$
    – StupidWolf
    May 31 '20 at 8:52
  • 1
    $\begingroup$ 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 $\endgroup$
    – StupidWolf
    May 31 '20 at 8:54
  • 1
    $\begingroup$ 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 $\endgroup$
    – StupidWolf
    May 31 '20 at 8:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.