R/FlagCluster.R

In Oscope: Oscope - A statistical pipeline for identifying oscillatory genes in unsynchronized single cell RNA-seq

#' @title Flag gene clusters with small within-cluster phase differences and/or small within-cluster sine scores
#' @usage FlagCluster(SineRes, KMRes, Data, qt, thre=pi/4, qtincluster=.5, qtinpermu=.9 ,Seed=1)
#' @param SineRes output of OscopeSine() function
#' @param KMRes output of KMRes() function
#' @param Data a gene-by-sample (isoform-by-sample) matrix indicating the rescaled expression of genes/isoforms.
#' all values should be between [-1, 1].
#' @param qt,thre Define a gene pair's linear score as 
#' min(eta, pi-eta), in which eta is defined as phase shift mod pi.
#' A cluster will be flagged if the 
#' qt th quantile of within-cluster linear score is less than thre. 
#' @param qtincluster,qtinpermu
#' To define clusters with small within-cluster sine scores, for each cluster we generate
#' permuted data of these genes (different cell permutation for each gene). 
#' We calculate the within-cluster sine scores within
#' the cluster of permuted genes, then infer whether the sine scores
#' in the cluster of interest are greater than those generated by the permuted genes.
#' A cluster will be flagged if its qtincluster th quantile in the original data is less than 
#' its qtinpermu th quantile in permuted data.
#' @param Seed seed
#' @return Output: RemoveID: a vector of cluster numbers that are flagged;
#' SineCompreList: sine score and permuted sine score for each cluster;  LinearList: linear score of each cluster
#' @examples aa <- sin(seq(0,1,.1))
#' bb <- sin(seq(0.5,1.5,.1))
#' cc <- sin(seq(0.9,1.9,.1))
#' tmp <- matrix(sin(rnorm(330)),ncol=11)
#' rownames(tmp) <- paste0("tmp",1:30)
#' Dat <- rbind(aa, bb, cc, tmp)
#' res1 <- OscopeSine(Dat)
#' res2 <- OscopeKM(res1, quan=.8, maxK=5)
#' res <- FlagCluster(res1, res2, Dat)
#' @author Ning Leng

FlagCluster <- function(SineRes, KMRes, Data, qt=.9,thre=pi/4, qtincluster=.5, qtinpermu=.9, Seed=1){


set.seed(Seed)
expect_is(SineRes, "list")
expect_is(KMRes, "list")
Len <- length(KMRes)
ShiftList <- sapply(1:Len,function(i)SineRes$ShiftMat[KMRes[[i]],KMRes[[i]]],simplify=FALSE)
for(i in 1:Len)diag(ShiftList[[i]]) <- NA
ShiftList_new <- sapply(1:Len,function(i){
mm <- ShiftList[[i]] %% pi
M0 <- abs(mm-0)
M90 <- abs(mm-pi)
out <- ifelse(M0<M90,M0,M90)
},
simplify=FALSE)

message("========= by sine score =========")

P <- rep(NA,Len)
PerList <- vector("list",Len)
for(i in 1:Len){
	per <- PermuCut(Data[KMRes[[i]],], length(KMRes[[i]]))
#for(i in 1:Len){	
	Incluster0 <- SineRes$SimiMat[KMRes[[i]],KMRes[[i]]]
	Incluster <- Incluster0[upper.tri(Incluster0)]
	Inpermu <- per[[1]][upper.tri(per[[1]])]
	message("cluster ",i)
	print(quantile(Incluster))
	message("cluster ",i," permuted data ")
	print(quantile(Inpermu))
	P[i] <- ifelse(quantile(Incluster,qtincluster)>quantile(Inpermu,qtinpermu), FALSE, TRUE)
	PerList[[i]] <- cbind(Incluster, Inpermu)
	#P[i] <- suppressWarnings(t.test(Inpermu, Incluster, alternative = "less")$p.value)
  #P[i] <- suppressWarnings(ks.test(Inpermu, Incluster)$p.value)
}
Pmark <- which(P)
#Pmark <- which(P>pvalthre)
message("\n flagged cluster:", Pmark)

message("========= by phase shift =========")

V <- rep(NA, Len)
for(i in 1:Len){
message("cluster ",i)
tmp <- as.vector(ShiftList_new[[i]])
tmpnoNA <- tmp[which(!is.na(tmp))]
print(summary(tmpnoNA))
V[i] <- quantile(tmpnoNA, qt)
}
expect_is(V,c("numeric","integer"))
Vmark <- which(V<thre)
message("flagged cluster:", Vmark)

Out <- list(FlagID=unique(c(Vmark,Pmark)), SineCompareList=PerList,LinearList=ShiftList_new,
						FlagID_bysine = Pmark, FlagID_byphase=Vmark)
}