Setting Up a Single Cell RNA Sequencing Pipeline That Actually Works
Single Cell Analysis has become the default approach for most labs studying heterogeneous tissues. The promise is simple: you can see variation between individual cells rather than getting averaged-out signals from a bulk homogenate. The reality is messier than the brochures suggest. I have run this workflow enough times to know where it breaks and how to fix it. Start with the raw data. You will receive FASTQ files from your sequencing facility, usually demultiplexed by sample index. Don't skip the QC check on those reads. I once processed a full batch only to discover 40% of the reads were adapter contamination because the facility had underclipped the reads. FastQC will tell you this immediately. Spending 10 minutes on that check saved me two days of wasted compute time. For alignment, Cell Ranger from 10x Genomics is the standard for Chromium data. It takes roughly 30 minutes to 2 hours depending on library complexity and your hardware. If you are working with Smart-seq2 or other full-length protocols, use Alevin or Salmon in single-cell mode instead. They are significantly faster and memory-efficient. Alevin can process a typical dataset in about 15 minutes on a standard 16-core machine with 32GB RAM, compared to hours with traditional aligners like STARsolo.
The Filtering Step Everyone Gets Wrong
Most tutorials suggest filtering based on mitochondrial percentage, gene count, and UMIs. The standard cutoffs are mitochondrial reads above 20% indicating dead cells, fewer than 200 genes suggesting empty droplets or ambient RNA, and more than 2500 genes potentially marking doublets. These are reasonable starting points, but they are not universal rules. Here is the problem I ran into last year. I was analyzing neuronal tissue, and neurons naturally have very high mitochondrial content. Applying a 20% mitochondrial cutoff removed nearly half my neuron population before any downstream analysis. The workaround was to set the mitochondrial threshold based on the bimodal distribution in my data rather than an arbitrary number. I used the ggplot2 library in R to visualize the distribution and identified a natural gap around 35% mitochondrial reads. Everything below that threshold was retained. This preserved the genuine neuronal signal while still removing the truly compromised cells. Another thing nobody warns you about is ambient RNA correction. When cells lyse during preparation, their RNA floats around and gets captured in other droplets. Libraries like SoupX or decontX can estimate and remove this contamination. I found that not applying ambient RNA correction to my immune cell datasets resulted in false differential expression calls, particularly for genes that are highly expressed in myeloid cells contaminating lymphocyte samples.
Clustering and Dimensionality Reduction
After filtering, you normalize the data and identify highly variable genes. Seurat's FindVariableFeatures function works well for this. The default selection of 2000 variable genes is a reasonable starting point. I usually check the dispersion plot to confirm that the selected genes show genuine biological variance rather than technical noise. For dimensionality reduction, PCA is standard. You want enough principal components to capture biological structure without including too much noise. The ElbowPlot in Seurat helps here, but I also run clustree to see how clusters change across different resolution parameters before committing to a final clustering. This prevents the common mistake of picking a resolution that either over-splits or under-splits your cell populations. UMAP and t-SNE are both viable for visualization. UMAP preserves more global structure, which matters when you are trying to understand relationships between cell states. t-SNE distorts distances more aggressively but can sometimes reveal finer substructure. I use both. UMAP for the big picture and t-SNE when I need to resolve closely related subtypes.
Get the Full Details
Annotation and Marker Identification
Cell type annotation is where many projects stall. Automated tools like SingleR can provide a first pass using reference datasets, but manual curation is essential. SingleR tends to misannotate rare cell types and cells in transition states. I cross-reference automated predictions with known marker genes from CellMarker and PanglaoDB, then verify by inspecting feature plots for each candidate marker. When identifying marker genes for a cluster, default parameters often miss important signals. I adjust the min.pct parameter to require expression in at least 10% of cells rather than the default 25%, and I set logfc.threshold to 0.25 instead of 0.25 or higher. This catches markers that are biologically meaningful but modestly expressed. The tradeoff is more false positives, so I always validate with independent methods before drawing conclusions.
Handling Batch Effects in Single Cell Analysis
Batch effect is probably the most underestimated problem in this field. If you process samples on different days, with different reagent lots, or from different donors, the technical variation can easily exceed the biological signal. Harmony and BBKNN are effective correction methods. Harmony integrates faster and works well for large datasets, typically processing 50,000 cells in under a minute. BBKNN is more memory-efficient for very large datasets but can over-correct when biological variation correlates with batch structure. I learned this the hard way when I integrated datasets from three different sequencing runs and noticed that cell type proportions shifted dramatically after integration. The correction had partially removed genuine biological differences between samples. The solution was to run integration diagnostics using kBET and silhouette scores before and after correction to quantify how much biological structure was being preserved.
Downstream Analysis and Validation
Differential expression between conditions should not rely solely on bulk-like methods applied to single cell data. Wilcoxon rank-sum tests implemented in Seurat are conservative but reliable. For more power, consider MAST or DESeq2 with pseudo-bulk aggregation. The pseudo-bulk approach, where you sum counts per cell type per sample before running differential expression, generally produces more reproducible results and fewer false positives. It also properly accounts for sample-level variability rather than treating each cell as an independent observation, which is a statistical error I see constantly in published work. Trajectory inference with Monocle3 or Slingshot can reveal developmental or transition pathways. These tools assume a branching or linear structure, which may not exist in your data. I always validate trajectory hypotheses with independent markers and, when possible, with spatial transcriptomics data to confirm that the inferred transitions correspond to actual spatial arrangements in the tissue. The field moves fast. New methods emerge regularly, and what was standard two years ago is often superseded. The core principles remain the same: check your data at every step, don't blindly apply default parameters, and validate your findings. Single Cell Analysis is powerful, but it demands careful attention to detail at every stage of the workflow.
