### 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)