library(Seurat)
library(DoubletFinder)
# Step 1.1: Load 10x data
data <- Read10X(data.dir="{tenx_dir}")
so <- CreateSeuratObject(counts=data, project="Dormancy", min.cells=3)
# Step 1.2: Calculate QC metrics
so[["percent.mt"]] <- PercentageFeatureSet(so, pattern="^mt-")
so[["percent.rb"]] <- PercentageFeatureSet(so, pattern="^Rp[sl]")
# Step 1.3: Filter cells
so <- subset(so, subset=nCount_RNA > {min_umi} & nFeature_RNA > {min_genes} & percent.mt < {max_mt})
# Step 1.4: Normalize and find variable features
so <- NormalizeData(so)
so <- FindVariableFeatures(so, selection.method="vst", nfeatures=2000)
# Step 1.5: Doublet removal with DoubletFinder
so <- ScaleData(so)
so <- RunPCA(so, features=VariableFeatures(object=so))
so <- RunUMAP(so, dims=1:20)
pK_data <- paramSweep(so, PCs=1:20, sct=FALSE)
pK_df <- summarizeSweep(pK_data, GT=FALSE)
bc_pK <- find.pK(pK_df)
optimal_pK <- as.numeric(as.character(bc_pK$pK[which.max(bc_pK$BCmetric)]))
so <- doubletFinder(so, PCs=1:20, pN=0.25, pK=optimal_pK, nExp=round({doublet_rate}*ncol(so)), reuse.pANN=FALSE)
df_col <- grep("DF.classifications", colnames(so@meta.data), value=TRUE)[1]
so <- so[, so@meta.data[[df_col]] == "Singlet"]
# Save
saveRDS(so, "01_seurat_qc.rds")
print(paste("Cells after QC:", ncol(so)))library(Seurat)
library(SingleR)
library(celldex)
# Step 2.1: Load data
so <- readRDS("01_seurat_qc.rds")
# Step 2.2: PCA and UMAP
so <- ScaleData(so)
so <- RunPCA(so, features=VariableFeatures(object=so), npcs=50)
so <- FindNeighbors(so, dims=1:{pc_dims})
so <- FindClusters(so, resolution={resolution})
so <- RunUMAP(so, dims=1:{pc_dims})
# Step 2.3: SingleR annotation
ref <- MouseRNAseqData()
pred <- SingleR(test=GetAssayData(so, layer="data"), ref=ref, labels=ref$label.fine)
so$cell_type <- pred$labels
# Step 2.4: Visualize
DimPlot(so, reduction="umap", label=TRUE, repel=TRUE) + NoLegend()
ggsave("02_umap_clusters.pdf", width=8, height=6)
DimPlot(so, reduction="umap", group.by="cell_type", label=TRUE, repel=TRUE) + NoLegend()
ggsave("02_umap_celltypes.pdf", width=10, height=6)
# Step 2.5: Find markers
markers <- FindAllMarkers(so, only.pos=TRUE, min.pct=0.25, logfc.threshold=0.25)
top_markers <- markers %>% group_by(cluster) %>% top_n(n=10, wt=avg_log2FC)
write.csv(top_markers, "02_cluster_markers.csv")
saveRDS(so, "02_seurat_clustered.rds")library(Seurat)
library(tidyverse)
so <- readRDS("02_seurat_clustered.rds")
# Step 3.1: DEG per cell type
Idents(so) <- "cell_type"
all_types <- unique(so$cell_type)
deg_list <- list()
for(ct in all_types) {
cells <- subset(so, idents=ct)
if(length(unique(cells${group_col})) >= 2) {
deg <- FindMarkers(cells, ident.1="{ident1}", ident.2="{ident2}",
group.by="{group_col}", logfc.threshold=0.25)
deg$cell_type <- ct
deg$gene <- rownames(deg)
deg_list[[ct]] <- deg
}
}
all_deg <- bind_rows(deg_list)
write.csv(all_deg, "03_celltype_DEG.csv")
# Step 3.2: Metabolic scoring
oxphos_genes <- c("Cox4i1","Cox5a","Cox5b","Cox6c","Atp5f1a","Atp5f1b","Ndufa1","Sdhb","Uqcrq")
glycolysis_genes <- c("Hk2","Pkm","Pfkp","Eno1","Gapdh","Ldha","Pgk1")
so <- AddModuleScore(so, features=list(oxphos_genes), name="OXPHOS_score")
so <- AddModuleScore(so, features=list(glycolysis_genes), name="Glycolysis_score")
# Step 3.3: Visualization
VlnPlot(so, features=c("OXPHOS_score1","Glycolysis_score1"), pt.size=0, ncol=2)
ggsave("03_metabolic_scores.pdf", width=12, height=5)
saveRDS(so, "03_seurat_scored.rds")library(monocle3)
library(Seurat)
# Step 4.1: Convert Seurat to Monocle
cds <- as.cell_data_set(so)
# Step 4.2: Preprocess and reduce
cds <- cluster_cells(cds, reduction_method="UMAP")
cds <- learn_graph(cds)
# Step 4.3: Order cells
root_cells <- colnames(cds)[cds@colData$cell_type == "{root_cell}"]
cds <- order_cells(cds, root_cells=root_cells)
# Step 4.4: Plot pseudotime
plot_cells(cds, color_cells_by="pseudotime", label_cell_groups=FALSE,
label_leaves=FALSE, label_branch_points=FALSE, graph_label_size=1.5)
ggsave("04_pseudotime.pdf", width=8, height=6)
# Step 4.5: Differential expression along trajectory
deg_graph <- graph_test(cds, neighbor_graph="principal_graph", cores=4)
top_degs <- deg_graph %>% arrange(q_value) %>% head(100)
write.csv(top_degs, "04_trajectory_DEG.csv")
# Step 4.6: Gene dynamics
genes_to_plot <- c("Ucp1","Pparg","Cidea","Prdm16","Cebpb")
plot_genes_in_pseudotime(cds[rowData(cds)$gene_short_name %in% genes_to_plot,],
color_cells_by="condition", min_expr=0.1)
ggsave("04_gene_pseudotime.pdf", width=10, height=8)