Seurat (R) tutorial¶

How to run Ouroboros on a Seurat (R) object¶

Ouroboros is a python package, so running it on a Seurat object requires saving your raw counts as a CSV file and running the command-line implementation of Ouroboros on the saved CSV. The ouroboros output can then be loaded into back into R.

Dependencies:

Seurat
SeuratObject
Matrix
Zellconverter
SingleCellExperiment

In [ ]:

Copied!





suppressPackageStartupMessages({
  library(Seurat)
  library(SeuratObject)
  library(Matrix)
  library(zellkonverter)
  library(SingleCellExperiment)
  library(SeuratData)
})
suppressPackageStartupMessages({
  library(Seurat)
  library(SeuratObject)
  library(Matrix)
  library(zellkonverter)
  library(SingleCellExperiment)
  library(SeuratData)
})

We'll demonstrate Ouroboros using pbmc3k, a standard dataset of 2,700 PBMCs (~13,700 genes) that ships as raw counts. We subset it to 100 cells to remove computational time and written object size.

In [10]:

Copied!

InstallData("pbmc3k")
data("pbmc3k")

pbmc3k   # ~2700 cells x ~13714 genes, raw counts in the RNA assay
InstallData("pbmc3k")
data("pbmc3k")

pbmc3k   # ~2700 cells x ~13714 genes, raw counts in the RNA assay

Installing package into ‘/projects/pangen/analysis/hmac/applications/miniconda3/lib/R/library’
(as ‘lib’ is unspecified)

In [12]:

Copied!

set.seed(330)
cells <- colnames(GetAssayData(pbmc3k, assay = "RNA", layer = "counts"))
pbmc100 <- pbmc3k[, sample(cells, 100)]
set.seed(330)
cells <- colnames(GetAssayData(pbmc3k, assay = "RNA", layer = "counts"))
pbmc100 <- pbmc3k[, sample(cells, 100)]

Next write the raw counts as a Scanpy h5ad object to run ouroboros on. p.s. this tends to take a long time (~ 15 minutes)

In [17]:

Copied!





# Seurat -> SingleCellExperiment
sce <- as.SingleCellExperiment(pbmc100)

# confirm a raw-counts assay is present
assayNames(sce)          # expect "counts" (and usually "logcounts")

# write .h5ad with RAW counts in .X
writeH5AD(sce, "R_example_output/pbmc_small.h5ad", X_name = "counts")
# Seurat -> SingleCellExperiment
sce <- as.SingleCellExperiment(pbmc100)

# confirm a raw-counts assay is present
assayNames(sce)          # expect "counts" (and usually "logcounts")

# write .h5ad with RAW counts in .X
writeH5AD(sce, "R_example_output/pbmc_small.h5ad", X_name = "counts")

'counts'
'logcounts'

There is still some formatting of that h5ad that needs to be done so Ouroboros will recognize it. Run this in a terminal (replace the h5ad path with your own):

conda activate ouroboros_env

# go to your output dir
cd R_example_output

python - <<'PY'
import h5py, numpy as np, anndata as ad
path = "pbmc_small.h5ad"
vlen = h5py.special_dtype(vlen=str)
with h5py.File(path, "r+") as f:
    if "layers/raw_counts" not in f:
        f.copy("X", "layers/raw_counts")

    # strip modern 'dict' tags (already done, but harmless to repeat)
    def fix_dict(name, obj):
        if isinstance(obj, h5py.Group):
            et = obj.attrs.get("encoding-type")
            et = et.decode() if isinstance(et, (bytes, bytearray)) else et
            if et == "dict":
                del obj.attrs["encoding-type"]; obj.attrs.pop("encoding-version", None)
    f.visititems(fix_dict)

    # flatten obs/var to index-only: drop categorical/group columns the old reader can't parse
    for key in ("obs", "var"):
        if key not in f: continue
        g = f[key]
        idx = g.attrs.get("_index", "_index")
        idx = idx.decode() if isinstance(idx, (bytes, bytearray)) else idx
        for col in list(g.keys()):
            if col != idx:
                del g[col]
        g.attrs.create("column-order", np.array([], dtype=vlen))
PY

Double check it worked with this:

python - <<'PY'
import anndata as ad
a = ad.read_h5ad("pbmc_small.h5ad")
print("OK:", a.shape, "| genes:", list(a.var_names[:3]), "| raw_counts:", "raw_counts" in a.layers)
PY

Then run Ouroboros on test object on command line

Call:

ouroboros \
    --data pbmc_small.h5ad \
    --data_type h5ad \
    --outdir .

Read Ouroboros output into Seurat object¶

In [ ]:

Copied!





emb <- read.csv(
  "R_example_output/ouroboros_embeddings_pseudotimes.csv",
  row.names = 1,          # first column is the cell_id index pandas wrote
  check.names = FALSE     # keep column names like dim1/dim3 exactly as-is
)

# sanity check: do the IDs line up with the Seurat object?
head(rownames(emb))
head(colnames(pbmc100))
mean(rownames(emb) %in% colnames(pbmc100))   # want this to be 1 

# add all Ouroboros columns to the cell metadata
pbmc100 <- AddMetaData(pbmc100, metadata = emb)

head(pbmc100@meta.data)
emb <- read.csv(
  "R_example_output/ouroboros_embeddings_pseudotimes.csv",
  row.names = 1,          # first column is the cell_id index pandas wrote
  check.names = FALSE     # keep column names like dim1/dim3 exactly as-is
)

# sanity check: do the IDs line up with the Seurat object?
head(rownames(emb))
head(colnames(pbmc100))
mean(rownames(emb) %in% colnames(pbmc100))   # want this to be 1 

# add all Ouroboros columns to the cell metadata
pbmc100 <- AddMetaData(pbmc100, metadata = emb)

head(pbmc100@meta.data)

'AAACCGTGCTTCCG'
'AAAGAGACGCGAGA'
'AAGATTACCCGTTC'
'AAGCAAGAGGTGTT'
'AAGTTCCTCATTCT'
'AATGTTGATCTACT'

'AAACCGTGCTTCCG'
'AAAGAGACGCGAGA'
'AAGATTACCCGTTC'
'AAGCAAGAGGTGTT'
'AAGTTCCTCATTCT'
'AATGTTGATCTACT'

1

A data.frame: 6 × 14
	orig.ident	nCount_RNA	nFeature_RNA	seurat_annotations	RNA_snn_res.0.5	seurat_clusters	dim1	dim2	dim3	KNN_phase	cell_cycle_pseudotime	south	dormancy_depth	G0_classification
	<fct>	<dbl>	<int>	<fct>	<fct>	<fct>	<dbl>	<dbl>	<dbl>	<chr>	<dbl>	<chr>	<dbl>	<chr>
AAACCGTGCTTCCG	pbmc3k	2639	960	CD14+ Mono	1	1	-0.2796167	0.9579706	-0.06408316	G0	NA	True	-0.5078842	quiescence
AAAGAGACGCGAGA	pbmc3k	3033	1058	CD14+ Mono	1	1	-0.7858221	0.5954063	0.16725660	G0	NA	True	-0.8482282	senescence
AAGATTACCCGTTC	pbmc3k	2762	928	CD14+ Mono	1	1	-0.6899688	0.5372884	0.48504025	G0	NA	True	-0.8712022	senescence
AAGCAAGAGGTGTT	pbmc3k	1684	857	NK	0	0	-0.5400499	-0.5703587	0.61889990	G2M	NA	True	-0.4522001	quiescence
AAGTTCCTCATTCT	pbmc3k	2334	862	FCGR3A+ Mono	1	1	-0.0417582	0.4772747	0.87776140	G0	NA	True	-0.4849494	quiescence
AATGTTGATCTACT	pbmc3k	2672	960	Memory CD4 T	0	0	-0.6850837	0.7279406	-0.02761939	G0	NA	True	-0.7152675	senescence

In [ ]:

Copied!





pbmc100 <- NormalizeData(pbmc100)
pbmc100 <- FindVariableFeatures(pbmc100)
pbmc100 <- ScaleData(pbmc100)
pbmc100 <- RunPCA(pbmc100, npcs = 30)
pbmc100 <- RunUMAP(pbmc100, dims = 1:20, n.neighbors = 15)

DimPlot(pbmc100, reduction = "umap")
pbmc100 <- NormalizeData(pbmc100)
pbmc100 <- FindVariableFeatures(pbmc100)
pbmc100 <- ScaleData(pbmc100)
pbmc100 <- RunPCA(pbmc100, npcs = 30)
pbmc100 <- RunUMAP(pbmc100, dims = 1:20, n.neighbors = 15)

DimPlot(pbmc100, reduction = "umap")

Warning message in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, :
“pseudoinverse used at -2.019”
Warning message in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, :
“neighborhood radius 0.32006”
Warning message in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, :
“reciprocal condition number  0”
Warning message in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, :
“There are other near singularities as well. 0.031008”
Centering and scaling data matrix

PC_ 1 
Positive:  CST3, AIF1, LST1, FCER1G, TYROBP, CFD, FTL, NPC2, FTH1, S100A11 
	   LGALS1, FCN1, LYZ, CTSS, PSAP, S100A9, SERPINA1, CD68, RP11-290F20.3, SAT1 
	   BRI3, TYMP, SPI1, MNDA, COTL1, IFITM3, IFI30, MS4A7, FGR, OAZ1 
Negative:  MALAT1, RPS6, RPL13, LTB, ISG20, CD3E, CD7, CXCR4, RPLP0, CTSW 
	   LYAR, FKBP11, CCL5, CCR7, NKG7, GZMA, UBXN4, NOSIP, RARRES3, TRABD2A 
	   TRAT1, DENND2D, RPL10, S1PR4, HSPE1, BBC3, RPS14, PBXIP1, AQP3, TCF7 
PC_ 2 
Positive:  FCGR3A, UBLCP1, CST7, CETP, B2M, HAT1, EMR2, FGFBP2, GZMA, ABI3 
	   C2orf42, ZNF593, HSPA6, HIST3H2A, HS6ST1, ZBTB17, C1orf51, EXOC6B, PARP6, PPP1R12B 
	   NLRC4, CNNM3, RHOC, NKG7, GZMH, LILRA2, WRAP73, CTSW, GZMB, IFIT1 
Negative:  S100A8, LGALS2, CD14, GPX1, MS4A6A, CSF3R, S100A9, LYZ, FOLR3, S100A12 
	   RPL13, BNIP3L, FCN1, NCF1, CD79A, FCGR1A, HLA-DRB1, FES, PLBD1, GRN 
	   RNASE6, HLA-DRA, RBP7, AP1S2, IER3, BANK1, HLA-DQB1, SRGAP2, MS4A1, FPR1 
PC_ 3 
Positive:  GZMA, FGFBP2, NKG7, GZMB, FCRL6, CST7, S100A9, CTSW, PRF1, GNLY 
	   S100A8, XCL2, SPON2, GZMH, MATK, BPGM, S1PR5, TYROBP, CCL4, MS4A6A 
	   CCL5, APMAP, PDIA4, PTGDR, CD14, AGTRAP, LGALS2, B2M, S100A4, KLRD1 
Negative:  HLA-DQA1, HLA-DQB1, CD79A, EAF2, LINC00926, HLA-DQA2, MS4A1, CD79B, TCL1A, HLA-DPA1 
	   IGLL5, POU2F2, CD74, FCER2, CD180, CD72, HLA-DPB1, CD1C, HLA-DRA, HLA-DRB1 
	   CCDC50, TAMM41, HLA-DOB, P2RX5, TSPAN13, BLNK, ZNF766, PPP1R14A, PKIG, NOTCH2 
PC_ 4 
Positive:  FGFBP2, FCRL6, MATK, GZMB, MARCKSL1, GZMA, BPGM, CST7, HLA-DQA1, GNLY 
	   NKG7, HLA-DQB1, SPON2, HLA-DQA2, GZMH, TTC14, S1PR5, CD79B, EAF2, CD79A 
	   PRF1, LAT2, LINC00926, IGLL5, HLA-DPA1, HLA-DOB, HLA-DPB1, CD74, DOCK5, TBC1D17 
Negative:  TCF7, SAMD9, COX4I1, MTIF3, HSPB1, RGCC, JUNB, HPCAL1, KLHL22, MAGED2 
	   CD3E, TRABD2A, SPATS2L, SSBP1, VIM, TRAT1, SMAP1, TNFRSF4, HIST3H2A, PARP6 
	   HS6ST1, CNNM3, EXOC6B, PPP1R12B, C1orf51, ZBTB17, ZNF593, C2orf42, NLRC4, HSPA6 
PC_ 5 
Positive:  ZNF503, DHRS9, ZNF652, CLCN6, PHGDH, INTS3, MGAT5, AC079305.10, FCER1A, AGBL5 
	   CCDC88A, ENHO, TMCC2, MTIF2, CD1C, MGST2, GPX3, DHRS1, LOH12CR2, ROGDI 
	   CCT5, CYSTM1, TFPT, CCDC50, PHACTR1, LARP1B, HDAC2, SRM, COA1, NRGN 
Negative:  MS4A1, CD79B, P2RY10, TCL1A, SWAP70, PKIG, CD79A, TCOF1, EMR2, POLD4 
	   CD72, SCPEP1, HLA-DOB, KIAA0125, NAT9, SIGLEC10, DUSP6, PILRA, LINC00926, MT-ND2 
	   SEC62, UBE2E1, FCER2, AAMP, APH1B, CTD-2006K23.1, TSPAN13, PDLIM1, FCGR3A, LINC01013 

Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
18:25:41 UMAP embedding parameters a = 0.9922 b = 1.112

18:25:41 Read 100 rows and found 20 numeric columns

18:25:41 Using Annoy for neighbor search, n_neighbors = 15

18:25:41 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

18:25:41 Writing NN index file to temp file /tmp/RtmppLebTQ/file5728e428c0b74

18:25:41 Searching Annoy index using 1 thread, search_k = 1500

18:25:41 Annoy recall = 100%

18:25:43 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 15

18:25:45 Initializing from normalized Laplacian + noise (using RSpectra)

18:25:45 Commencing optimization for 500 epochs, with 1866 positive edges

18:25:45 Using rng type: pcg

18:25:47 Optimization finished

No description has been provided for this image

In [19]:

Copied!

# Can colour UMAPs by Ouroboros variables
FeaturePlot(pbmc100, features = "cell_cycle_pseudotime")
# Can colour UMAPs by Ouroboros variables
FeaturePlot(pbmc100, features = "cell_cycle_pseudotime")

Seurat (R) tutorial¶

How to run Ouroboros on a Seurat (R) object¶

Read Ouroboros output into Seurat object¶

Ouroboros

Navigation

Table Of Contents