update PatternLine

NAMESPACE

@@ -37,6 +37,7 @@ export(GetTDEseqAssayData)
 export(IsMatrixEmpty)
 export(NormDataTDEseq)
 export(PatternHeatmap)
+export(PatternLine)
 export(SetAssayData)
 export(TDEseq)
 export(Version)

R/TDEseq.R

@@ -214,9 +214,14 @@ tdeseq.default <- function(object,
     }
     dfTDEseqResults$ChangePoint<-ChangePoint
 
+	fit.inc <- do.call(cbind,res.tdeseq[seq(2,length(res.tdeseq),9)])
+	fit.dec <- do.call(cbind,res.tdeseq[seq(3,length(res.tdeseq),9)])
+	fit.cov <- do.call(cbind,res.tdeseq[seq(4,length(res.tdeseq),9)])
+	fit.con <- do.call(cbind,res.tdeseq[seq(5,length(res.tdeseq),9)])
 
-
-	return(dfTDEseqResults)
+	ModelFits<-GetModelFits(dfTDEseqResults,list(fit.inc,fit.dec,fit.cov,fit.con),stage.id)
+	out<-list(dfTDEseqResults=dfTDEseqResults,ModelFits=ModelFits)
+	return(out)
 }## end function 
 
 
@@ -307,7 +312,9 @@ tdeseq.Assay <- function(object,
 
 
 	## store the scaled data in the slot
-	object <- SetAssayData(object = object, slot = 'tde', new.data = new.data)
+	object <- SetAssayData(object = object, slot = 'tde', new.data = new.data$dfTDEseqResults)
+	object <- SetAssayData(object = object, slot = 'fit', new.data = new.data$ModelFits)
+	return(object)
 	## store top number of features in tde slot
 }## end func
 

R/Visualization.R

@@ -146,24 +146,27 @@ return(p)
 
 #' PatternLine: Line plot to show the pattern specific temporal genes. Please first install ggplot2.
 #' 
-#' @param seuobj A seurat object with UMAP embeddings been calculated
-#' @param features Genes to be shown in feature plot
+#' @param tdeobj A tdeseq object
+#' @param features Genes to be shown in line plot
 #' @author Yue Fan, Shiquan Sun
-
+#' @export
 PatternLine<-function (obj, feature.show = NULL, cols = NULL) 
 {
     if(is.null(cols))
 	{
 	cols=c('#E50C7D',"#E34627","#A22066","#A474A4","#2D8573","#E1DE15","#C16728","#2578BE","#738DC8","#C0C0C0", '#7d8d8e','#2a24d0','#a2292e','#274382','#838d36')
 	}
     suppressPackageStartupMessages(library("ggplot2"))
-	feature.show=feature.show[feature.show%in%rownames(obj@ModelFits)]
+	dat <- GetTDEseqAssayData(obj,slot='fit')
+	feature.show=feature.show[feature.show%in%rownames(dat)]
 	if(length(feature.show)==0)
 	{
 	stop(paste0("Genes to be shown are not specified!!"))
 	}
-    dat = obj@ModelFits[feature.show, ]
-    time = sort(unique(obj@Metadata$Time_Origin))
+
+    dat = dat[feature.show, ]
+	time=as.numeric(colnames(dat))
+
 	N=length(feature.show)
 	if(N!=1)
 	{

R/tdeseq.object.R

@@ -74,7 +74,8 @@ Assay <- setClass(
 		tde = "data.frame",
 		svg = "ANY",
 		de.truth = "ANY",
-		misc = "ANY"))
+		misc = "ANY",
+		fit = "matrix"))
 
 
 ###########################################
@@ -948,7 +949,8 @@ GetTDEseqAssayData.Assay <- function(object,
 	'tde',
 	'dr',
 	'svg',
-	'de.truth'), ...) {
+	'de.truth',
+	'fit'), ...) {
 	slot <- slot[1]
 	slot <- match.arg(arg = slot)
 	return(slot(object = object, name = slot))
@@ -1008,7 +1010,8 @@ SetAssayData.Assay <- function(object,
 	'dr',
 	'tde',
 	'svg',
-	'de.truth'), ...) {
+	'de.truth',
+	'fit'), ...) {
 
 	slot <- slot[1]
 	slot <- match.arg(arg = slot)

R/util.R

@@ -191,6 +191,23 @@ combined_pvalue <- unlist(lapply(res, ComputeACAT))
 return(combined_pvalue)
 }
 
+GetModelFits<-function(res_dat,res,stage)
+{
+sig_gene_idx=which(res_dat$SignificantDE=='Yes')
+ModelFits=matrix(0,nrow=length(sig_gene_idx),ncol=length(unique(stage)))
+rownames(ModelFits)=res_dat$gene[sig_gene_idx]
+colnames(ModelFits)=sort(unique(stage))
+pattern=c('Growth','Recession','Trough','Peak')
+for(i in sig_gene_idx)
+{
+for(j in sort(unique(stage)))
+{
+ModelFits[res_dat$gene[i],j+1]=mean(res[[match(res_dat$pattern[i],pattern)]][which(stage==j),i])
+}
+}
+return(ModelFits)
+}
+
 
 ComputeACAT <- function(Pvals, Weights=NULL){
  #### check if there is NA