5
$\begingroup$

I am fairly new to scRNAseq analysis and keep running into the same problem in the two datasets I am currently working with. We work with kidney inflammation and have two new mouse models, which we are now analysing in scRNAseq. I followed the usual Seurat workflow and ended up getting some nice UMAPs. With FindAllMarkers I was able to identify the clusters pretty well.

But once I do differential expression analysis (either with FindMarkers or pseudobulking) I got weird results.

Background There is a highly specialized epithelial celltype in the kidney (the podocyte) that specifically expresses a host of markers (like e.g. Nphs1 and Nphs2). They get downregulated when the kidney gets damaged (as it does in our disease models).

Problem

  1. I do not expect to see those markers to be present or regulated in any other celltype. But when I do the DE, those podocyte-specific marker genes pop up in almost all other clusters and are described to be highly regulated (downregulated, as I would expect in the pdoocytes).
  • So in my eyes, this points to "doublets of podocytes" and other celltypes being involved here. Maybe some damaged podocytes got into the same droplet as other celltypes and thus distort the analysis.
  1. I then used somewhat stringent filtering via scDblFinder* to remove a lot of doublets (~ 10 % of my dataset).
  • However, I still get the same problem.
  1. Then I analysed each batch (we have 3 in total) individually, but in each one this problem comes up.

Overall, I see 6.4% of doublets (that is 8,710 doublets in a total of 135,844 cells). I have two podocyte clusters with 1.7% doublets (that is 462 doublets in a total of 26,693 podocytes).

Some help / advice would be highly appreciated as I am becoming a little frustrated with this issue.

Thanks a lot!

*, Bioconductor package

$\endgroup$
0

3 Answers 3

2
$\begingroup$

An explanation in 3 parts

Final Edit: It turns out there was contamination in the data which you can read about in the comments, without being able to see it myself I could not know this. If you see this in an experiment and you're clustering from totally raw unfiltered data this maybe the case for you. Nevertheless as you subcluster within cell types of particular interest to get more refined marker genes the advice I provide below may be of use.

EDIT: The note above about looking at the violin plots may be useful, but again what your downstream analysis or experiment should impact your decision making.

TLDR: The presence of marker genes in other cell types is not surprising nor necessarily a problem. A marker gene that seurat finds are highly expressed in a cell type relative to the average in all cells not a gene that is uniquely expressed, if you are looking for uniquely expressed genes you will probably not find them.

1) A note on seurat

So most genes are expressed in most cells there is nothing wrong with your data or your analysis. Gene expression is very messy and the term "marker gene" is a bit of a misnomer. Note that other analysis pipelines like Monocle have different methods for determining what constitutes a marker gene.

Please consider this meta-analysis of gabanergic cerebral cortex neurons (i.e. inhibitory cortical neurons, interneurons etc etc).
https://ieeexplore.ieee.org/abstract/document/9669888

The fact is that most cell populations express a gradient of genes and that each gene likely has multiple roles in each cell. If your goal is to simply publish on expressed marker genes taking the top highest uniquely expressed genes which achieve significance in each cluster will be sufficient, just be sure to report ALL your settings and provide the meta data when you publish. If you are going on to do MERFISH or stain with the results you may want to be more selective.

FindaAllMarkers() uses a simple DEG analysis where it takes the average expression of each feature (the genes you selected, which I hope you did because it's not wroth looking at every gene for clustering just the top most variable) in the cluster against the average expression for these genes in all other clusters. Note that this means you are doing A TON of statistical tests so they use the bonferroni correction which would normally be considered to conservative but in this case is probably still too lenient. You should only be reporting on the very top 5 maybe 10 expressed genes. If your lucky one of these will be noticeably more prominent than the others.

2) Alternatives to seurat Monocle package has a different definition which arguably is better: http://cole-trapnell-lab.github.io/monocle-release/docs/#classifying-and-counting-cells though this package is less widely used. The monocle score accounts for how widely expressed across all clusters a marker is and the more unique to a cluster a gene is the greater its score in addition to its actual expression strength.
There is also a package called pandora, but I know next to nothing about it. Read up on them and decide which has features you like best.

Lastly I think part of your problem is just wanting biology to be a cleaner science than it is. The reality is that things we call cell types are simply transient states a population of cells passes through during development and various environmental perturbations, the reality is each cell is a dynamic system and you are sampling a portion of these graded changes. Some states are more stable than others and we have taken to calling these cell types but remember you can always subcluster within a cluster and pull out more clusters which can themselves be subclustered and so on and so forth. It calls into question the definition of a cell type even more and to be honest at this point finding cell types is somewhat of an art, the science is yet to be settled.

I encourage you to do more reading on gene expression within single cells and the definition of a cell type (or if such a thing even exists). Here is some reading to get you started.

Cembrowski, Mark S., and Vilas Menon. "Continuous variation within cell types of the nervous system." Trends in Neurosciences 41, no. 6 (2018): 337-348. https://www.sciencedirect.com/science/article/pii/S0166223618300596

Muñoz-Manchado, A. B., Gonzales, C. B., Zeisel, A., Munguba, H., Bekkouche, B., Skene, N. G., ... & Hjerling-Leffler, J. (2018). Diversity of interneurons in the dorsal striatum revealed by single-cell RNA sequencing and PatchSeq. Cell reports, 24(8), 2179-2190. https://www.sciencedirect.com/science/article/pii/S2211124718311562

If you really want to lean into this you should start digging into the concept of meta-cells, there is a package called Metacell2 which may be of use, be sure to read the documentation and methods section for Metacell (the first version) as Metacell2 won't make much sense without it:

Harris, Benjamin D., Megan Crow, Stephan Fischer, and Jesse Gillis. "Single-cell co-expression analysis reveals that transcriptional modules are shared across cell types in the brain." Cell Systems 12, no. 7 (2021): 748-756. -this is not metacell but is a meta analysis of cell markers using metacells with some pretty surprising and counter intuitive results about marker genes

Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.

Ensure you read all documentation carefully and understand the stats they are using. The methods use are a specific definition. Also a real marker gene will be robust to multiple methods of identifying them, if they disappear based on a change in normalization or in other packages you may have your answer there. Bear in mind if you publish your results should be reproducible and people will need to know exactly how you got your results so share your code and metadata please.

3) What do you want to get from this analysis?

As I mentioned earlier if your goal is to get published, simply finding some markers in seurat may suffice. If your goal is to do further experiments based on these results consider your hypothesis carefully. If you want more certainty pull some similar kidney scRNA seq data from GEO and compare your results to their or integrate your data with theirs seurat has a tutorial you can follow along with to do this https://satijalab.org/seurat/articles/integration_introduction.html#performing-integration-on-datasets-normalized-with-sctransform-1 Note you can also use this pipeline to do batch correction, that may help you. The sct transform give pretty good results.

Sorry nature isn't like a biology textbook =(, I'm not saying that to be facetious, I just want to remind you that living things are very complex. Lean into that complexity and make some interesting discoveries or just accept it and dig out some results you can publish. Talk with your supervisor and explain the issue. Hope this helps!

4) I lied sorry this is a response in 4 parts, make sure you report exactly what you did and share your metadata and pipeline

Did I mention your results should be reproducible? How can one do this? Here are some tips. Make sure your keeping your analysis on github, so you can keep track of all the changes it goes through and revert to previous versions if you change things. Keep all your code, its worth having a file that is full of old code you threw away cause it was messy or you didn't end up using it. Not broken code but stuff you don't think is useful. Report your final results, share your data and make sure a monkey could follow you instructions to reproduce your results as well as being able to double check they get the same thing. Use a conda environment to ensure software updates will not destroy your code for future lab members or people who read your future paper. Save all file paths as variables at the top of your code in your scripts so people can easily change the file structure without having to go through your code and struggle finding which read.csv() functions they need to edit.

$\endgroup$
3
  • 1
    $\begingroup$ Thanks so much for your detailed answer Angus, this was really very helpful for me as a beginner. As mentioned in my answer below, we discussed the dataset with more experienced bioinformaticians and it seems like there was a big amount of ambient RNA contamination in the dataset. Correction for this with the DecontX package cleaned things up a good bit (even though of course there is still some celltype-specific transcripts floating around in the other celltypes which is expected to a certain degree, as you mentioned). $\endgroup$
    – Smeerlap
    Dec 16, 2022 at 7:34
  • 1
    $\begingroup$ Glad you caught that confound. I work in brain tissue, which is much messier often people are doing lots of sub clustering, i.e. they do the initial round of clustering but then select only the neurons and redo the process. Then they might pull out only the excitatory or inhibitory cells (the two major types of neurons) and subcluster within just those cells again. As you imagine it gets messy. I can see how if at the initial round of clustering you are seeing this extreme level of non-specificity that would suggest contamination. Glad you caught it. $\endgroup$ Dec 16, 2022 at 9:12
  • $\begingroup$ Hi @Smeerlap could you please very kindly, upvote and/or mark this answer as "accepted". It is a very cool answer as you've said. Seurat vs. Monocle was very helpful for me and first time I've seen a comparison on this site. $\endgroup$
    – M__
    Dec 16, 2022 at 12:41
2
$\begingroup$

The violin plots show only little expression in non-podocyte clusters, but the podocyte marker transcripts are only the most obvious manifestation of this problem. Also, small but significant amounts of endothelial and other celltype transcripts show up in most of the other clusters. Identification of the celltype works fine but these contaminating counts throw off most of the downstream analyses.

By now, we did manage to find a solution for the problem, though: It turns out the problem was due to high amounts of ambient RNA contamination. I tried three different packages to remove the contaminating counts (DIEM, DecontX and SoupX) and found that DecontX (and SoupX) worked best for my dataset.

$\endgroup$
1
$\begingroup$

You should make violin plots of the genes of interest (x-axis is your different clusters (i.e. cell-types) and y-axis is expression). Maybe they are differentially expressed but at the same time not expressed very highly.

Also why does their expression matter in other cell types? Do you have strong evidence to expect that these genes would only be expressed in podocytes?

$\endgroup$
1
  • $\begingroup$ Hi! Thanks for your help. It seems like the problem was due to ambient RNA contamination (see my post below). $\endgroup$
    – Smeerlap
    Dec 16, 2022 at 7:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.