By Nils Trost

This calculates the cross correlation for all genes with a gene of interest.

# y is a cell x gene matrix, as can be found in the @data or @counts slots of a Seurat assay object. (dgCMatrix)
y <- srt@assays$RNA@data
 
# Get the expression of the gene of interest (POU5F1 in this case) in the cells ordered by the pseudo time lambda.
ts.GoI <- y["POU5F1", order(lambda, na.last = NA)]
 
# Calculate the cross correlation for all genes with the gene of interest.
cc.GoI <- apply(y, 1, function(z) {
  # order by pseudo time lambda
  ts.z <- z[order(lambda, na.last = NA)]
  # calculate the cross corellation along the entire length of the pseudotime (lag.max)
  ct <- ccf(ts.GoI, ts.z, lag.max = y@Dim[2], plot = FALSE)
  ac <- as.vector(ct$acf)
  names(ac) <- as.vector(ct$lag)
  return(ac)
})
 
# Filter out NAs and NaNs
cc.GoI.cln <- cc.GoI[, !colSums(is.na(cc.GoI)) > 0]
 
# Assign a simple scoring get to genes that peak at a similar timepoint as the gene of interest by giving
# more weight to cross correlation close to 0 lag and calculating the mean of the cross correlation.
cc.score <- apply(cc.GoI.cln, 2, function(x) {
  pos <- as.numeric(names(x))
  weight <- dnorm(pos, mean = 0, sd = 50) * 1
  return(mean(x * weight))
})
 
# order and get the top 100 genes.
ordered.cc <- cc.score[order(cc.score, decreasing = T)]
top.cc.genes.GoI <- names(ordered.cc)[2:100]