This page generates analysis code in your browser. It does not upload data or execute the workflow on the Dormanomics server. Review the generated parameters, software references and input requirements, then run and validate the code in your own computing environment.
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)