How to make PHATE_HARMONY reduction ?

[agctu-zhang]

Hello. I read your article and would like to learn how to get phate reduction using harmony data, however, I have some confusion while reading the source code. It is as follows:

# Run PHATE on harmony values ####

tmp.phate <- phate(
  t(GetAssayData(myo.seurat,slot="data")),
  gamma=1,
  n.jobs = 25
) 
colnames(tmp.phate$embedding) <- paste0("phate_", 1:2)
tmp.phate$params$data <- NULL

myo.seurat[["phate"]] <- CreateDimReducObject(
  embeddings = tmp.phate$embedding,
  key = "phate_", 
  assay = 'RNA',
  misc = tmp.phate$params
)

Here, perform phate used expression matrix, not used harmony embeddings?

#     binning on phate_harmony values ####

nbins = 25
breaks = seq(
  min(myo.seurat@reductions$phate_harmony@cell.embeddings[,1]),
  max(myo.seurat@reductions$phate_harmony@cell.embeddings[,1]),
  diff(range(myo.seurat@reductions$phate_harmony@cell.embeddings[,1]))/nbins
)
binnames = as.factor(1:nbins)

myo.seurat$phate1.bins <- cut(
  myo.seurat@reductions$phate_harmony@cell.embeddings[,1]+10^-100,
  breaks = breaks,
  labels=binnames
)
myo.seurat$phate1.bins[is.na(myo.seurat$phate1.bins)] <- 1

I see "pahte_harmony" appearing here for the first time, but I can't find how to add "phate_harmony", is from "temp.phate"?

I would be very grateful if you could answer my questions!

[mckellardw]

My apologies! That snippet of code must have been accidentally deleted. We repeated the PHATE+Harmony analysis on other cell types here (https://github.com/mckellardw/scMuscle/blob/main/R_scripts/5S_subcompartment_analysis.R) though. The code you want can be found starting in line 335 (this is the code used for the final manuscript, sorry for the confusion!)

myo.seurat@reductions$harmony@cell.embeddings <- 
  myo.seurat@reductions$harmony@cell.embeddings[,myo.seurat@reductions$harmony@stdev %>% order(decreasing = T)]

# subset out cells that don't line up
remove.cells <- c(
  "oprescu_D2_ACTTTCACATAGGAGC", "oprescu_D2_CCTGCATAGCCTGGAA", "oprescu_D2_CTTCTCTCACATGAAA", 
  "oprescu_D2_GCACGTGGTACTTCCC", "oprescu_D2_TCTCACGGTGGACCTC", "oprescu_D2_TGGTAGTAGTGAACAT", 
  "oprescu_D2_TTACGTTAGCGTGAGT"
)
myo.seurat <- subset(
  myo.seurat,
  cells = Cells(myo.seurat)[!Cells(myo.seurat)%in%remove.cells]
)

# select high stddev dimensions

# n.pcs <- npcs(myo.seurat, reduction = "harmony")
# n.pcs = 32

harmony.dims <- myo.seurat@reductions$harmony@stdev %>% order(decreasing = T)
harmony.dims <- harmony.dims[1:35]

# run PHATE
tmp.phate <- phate(
  myo.seurat@reductions$harmony@cell.embeddings[,1:39], #cells x reduc dims
  gamma=1,
  n.jobs = 50
) 
colnames(tmp.phate$embedding) <- paste0("phate_harmony_", 1:2)
tmp.phate$params$data <- NULL

tmp.phate$embedding[,1] <- tmp.phate$embedding[,1]*-1 #flip x-axis

myo.seurat[["phate_harmony"]] <- CreateDimReducObject(
  embeddings = tmp.phate$embedding,
  key = "phate_harmony_", 
  assay = 'RNA',
  misc = tmp.phate$params
)

myo.seurat@reductions$phate_harmony@stdev <- 
  apply(myo.seurat@reductions$phate_harmony@cell.embeddings, 2, sd) #find std dev for phate vals

DimPlot(
  myo.seurat,
  cells=sample(colnames(myo.seurat)),
  reduction='phate_harmony',
  group.by = 'phate1.bins',
  pt.size = 0.5,
  # cols=unique(celltype.colors),
  label=F
)

[agctu-zhang]

Thank you for your reply, it totally solved my problem!