library(HSMMSingleCell)
data("HSMM_expr_matrix")
data("HSMM_gene_annotation")
data("HSMM_sample_sheet")

suppressMessages(library(monocle))
suppressMessages(library(stringr))
suppressMessages(library(plyr))
suppressMessages(library(netbiov))

# read in data
R<-read.table("data/mouse_embryo/rpkms_Deng2014_preinplantation.txt")
# read metadata with sample names as rownames
M<-read.table("data/mouse_embryo/Metadata_mouse_embryo.csv",header = T, sep=",",row.names=1)

# prepare pheno and feature data
num_cells <- apply(R,1,function(x) sum(x>1))

genes <- data.frame(gene_short_name = rownames(R),num_cells_expressed=num_cells)
rownames(genes)<-rownames(R)
pd <- new("AnnotatedDataFrame", data = M)
fd <- new("AnnotatedDataFrame", data = genes)

# Create monocle data set - in this case we use tobit distribution for the expression since we have rpkms and no spike-ins.
cds <- newCellDataSet(as.matrix(R),
                       phenoData = pd,
                       featureData = fd,
                       lowerDetectionLimit=1,
                       expressionFamily=tobit(Lower=1))  

# find the expressed genes
expressed_genes <- row.names(subset(fData(cds), num_cells_expressed >= 10))

# estimate size factors 
cds <- estimateSizeFactors(cds)

savefile <- "data/mouse_embryo/monocle/monocle_de_genes.Rdata"
if (file.exists(savefile)){
  load(savefile)
}else{
  diff_test_res <- differentialGeneTest(cds[expressed_genes,],
                                      fullModelFormulaStr="~Stage")
  save(diff_test_res,file=savefile)
}
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
length(ordering_genes)

#Now define the ordering genes
cds <- setOrderingFilter(cds, ordering_genes)

cds <- reduceDimension(cds, max_components=2)
#Now that the space is reduced, it???s time to order the cells using the orderCells function as shown below.
cds <- orderCells(cds)
#Once the cells are ordered, we can visualize the trajectory in the reduced dimensional space.
plot_cell_trajectory(cds, color_by="Stage")

# plot by "State" - Monocle???s term for the segment of the tree.
plot_cell_trajectory(cds, color_by="State")

# we clearly see that State1 should be the root
# so now we can reorder the cells  again
cds <- orderCells(cds, root_state=1)

# plot with coloring by pseudotime
plot_cell_trajectory(cds, color_by="Pseudotime")

savefile <- "data/mouse_embryo/monocle/pseudotime_genes.Rdata"
if (file.exists(savefile)){
  load(savefile)
}else{
  diff_test_res2 <- differentialGeneTest(cds[expressed_genes,], fullModelFormulaStr="~sm.ns(Pseudotime)")
  save(diff_test_res2,file=savefile)
}

head(diff_test_res2[order(diff_test_res$qval),])

# take top 8 genes
plotgenes <- rownames(diff_test_res2[order(diff_test_res$qval)[1:8],])
plot_genes_in_pseudotime(cds[plotgenes,], color_by="Stage",ncol=2)

sig_gene_names <- row.names(subset(diff_test_res, qval < 1e-5))
length(sig_gene_names)

# With a strict cutoff we still have quite many significant genes, hard to produce a heatmap with all of them.
# sample at random 50 genes and plot heatmap
sel.genes <- sample(sig_gene_names,50)
plot_pseudotime_heatmap(cds[sel.genes,],
                        num_clusters = 6,
                        cores = 1,
                        show_rownames = T)

savefile <- "data/mouse_embryo/monocle/monocle_beam_res.Rdata"
if (file.exists(savefile)){
  load(savefile)
}else {
  BEAM_res <- BEAM(cds, branch_point=1, cores = 1)
  save(BEAM_res,file=savefile)
}


BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
sig_gene_names <- row.names(subset(BEAM_res, qval < 1e-5))
length(sig_gene_names)

# sample at random 50 signifiant genes and plot heatmap
sel.genes <- sample(sig_gene_names,50)


#plot_genes_branched_heatmap(cds[sel.genes,], branch_point = 1,
#                            num_clusters = 4,
#                            cores = 1,
#                            use_gene_short_name = T,
#                            show_rownames = T)

# we can plot some genes that separates along the branch points
plot_genes_branched_pseudotime(cds[sel.genes[1:12],],branch_point=1,ncol=3)

# Estimate RNA counts
rpc_matrix <- relative2abs(cds)
# Now, make a new CellDataSet using the RNA counts
cds2 <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
                       phenoData = pd,
                       featureData = fd,
                       lowerDetectionLimit=0.5,
                       expressionFamily=negbinomial.size())

# find the expressed genes
expressed_genes <- row.names(subset(fData(cds2), num_cells_expressed >= 10))

# estimate size factors 
cds2 <- estimateSizeFactors(cds2)
cds2 <- estimateDispersions(cds2)

savefile <- "data/mouse_embryo/monocle/monocle_de_genes2.Rdata"
if (file.exists(savefile)){
  load(savefile)
}else{
  diff_test_res <- differentialGeneTest(cds2[expressed_genes,],
                                      fullModelFormulaStr="~Stage")
  save(diff_test_res,file=savefile)
}
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
length(ordering_genes)

#Now define the ordering genes
cds2 <- setOrderingFilter(cds2, ordering_genes)

cds2 <- reduceDimension(cds2, max_components=2)
cds2 <- orderCells(cds2)
plot_cell_trajectory(cds2, color_by="Stage")

disp_table <- dispersionTable(cds2)
ordering_genes <- subset(disp_table,
                         mean_expression >= 0.5 &
                         dispersion_empirical >= 1 * dispersion_fit)$gene_id

cds3 <- setOrderingFilter(cds2,ordering_genes)
cds3 <- reduceDimension(cds3, max_components=2)
cds3 <- orderCells(cds3)
plot_cell_trajectory(cds3, color_by="Stage")

sessionInfo()