M3Drop_Vignette.R
In M3Drop: Michaelis-Menten Modelling of Dropouts in single-cell RNASeq

### R code from vignette source 'M3Drop_Vignette.Rnw'

###################################################
### code chunk number 1: M3Drop_Vignette.Rnw:50-52
###################################################
library(M3Drop)
library(M3DExampleData)


###################################################
### code chunk number 2: M3Drop_Vignette.Rnw:59-64
###################################################
counts <- Mmus_example_list$data
labels <- Mmus_example_list$labels
total_features <- colSums(counts >= 0)
counts <- counts[,total_features >= 2000]
labels <- labels[total_features >= 2000]


###################################################
### code chunk number 3: M3Drop_Vignette.Rnw:72-73
###################################################
norm <- M3DropConvertData(counts, is.counts=TRUE)


###################################################
### code chunk number 4: M3Drop_Vignette.Rnw:78-79
###################################################
norm <- M3DropConvertData(log2(norm+1), is.log=TRUE, pseudocount=1)


###################################################
### code chunk number 5: M3Drop_Vignette.Rnw:90-101
###################################################
K <- 49
S_sim <- 10^seq(from=-3, to=4, by=0.05)
MM <- 1-S_sim/(K+S_sim)
plot(S_sim, MM, type="l", lwd=3, xlab="Expression", ylab="Dropout Rate", 
	xlim=c(1,1000))
S1 <- 10; P1 <- 1-S1/(K+S1);
S2 <- 750; P2 <- 1-S2/(K+S2);
points(c(S1,S2), c(P1,P2), pch=16, col="grey85", cex=3)
lines(c(S1, S2), c(P1,P2), lwd=2.5, lty=2, col="grey35")
mix <- 0.5;
points(S1*mix+S2*(1-mix), P1*mix+P2*(1-mix), pch=16, col="grey35", cex=3)


###################################################
### code chunk number 6: M3Drop_Vignette.Rnw:113-114
###################################################
M3Drop_genes <- M3DropFeatureSelection(norm, mt_method="fdr", mt_threshold=0.01)


###################################################
### code chunk number 7: M3Drop_Vignette.Rnw:151-152
###################################################
count_mat <- NBumiConvertData(Mmus_example_list$data, is.counts=TRUE)


###################################################
### code chunk number 8: M3Drop_Vignette.Rnw:161-167
###################################################
DANB_fit <- NBumiFitModel(count_mat)

# Smoothed gene-specific variances
par(mfrow=c(1,2))
stats <- NBumiCheckFitFS(count_mat, DANB_fit)
print(c(stats$gene_error,stats$cell_error))


###################################################
### code chunk number 9: M3Drop_Vignette.Rnw:180-181
###################################################
NBDropFS <- NBumiFeatureSelectionCombinedDrop(DANB_fit, method="fdr", qval.thres=0.01, suppress.plot=FALSE)


###################################################
### code chunk number 10: M3Drop_Vignette.Rnw:202-203
###################################################
HVG <- BrenneckeGetVariableGenes(norm)


###################################################
### code chunk number 11: M3Drop_Vignette.Rnw:213-215
###################################################
heat_out <- M3DropExpressionHeatmap(M3Drop_genes$Gene, norm, 
		cell_labels = labels)


###################################################
### code chunk number 12: M3Drop_Vignette.Rnw:227-230
###################################################
cell_populations <- M3DropGetHeatmapClusters(heat_out, k=4, type="cell")
library("ROCR") 
marker_genes <- M3DropGetMarkers(norm, cell_populations)


###################################################
### code chunk number 13: M3Drop_Vignette.Rnw:242-244
###################################################
head(marker_genes[marker_genes$Group==4,],20) 
marker_genes[rownames(marker_genes)=="Cdx2",] 


###################################################
### code chunk number 14: M3Drop_Vignette.Rnw:259-261
###################################################
heat_out <- M3DropExpressionHeatmap(NBDropFS$Gene, norm, 
		cell_labels = labels)


###################################################
### code chunk number 15: M3Drop_Vignette.Rnw:278-280
###################################################
heat_out <- M3DropExpressionHeatmap(HVG, norm, 
		cell_labels = labels)