## download an object of 9 PBMC samples sample.file.url = "https://genome.med.nyu.edu/results/external/iCellR/example2/my.obj.Robj" # download the file download.file(url = sample.file.url, destfile = "my.obj.Robj", method = "auto") ### load iCellR and the object library(iCellR) load("my.obj.Robj") #### QC my.obj <- qc.stats(my.obj, s.phase.genes = s.phase, g2m.phase.genes = g2m.phase, which.data = "raw.data") summary(my.obj@stats) png('plot1_QC_stats.png',width = 15, height = 6, units = 'in', res = 300) stats.plot(my.obj, plot.type = "all.in.one", out.name = "UMI-plot", interactive = FALSE, cell.color = "slategray3", cell.size = 1, cell.transparency = 0.5, box.color = "red", box.line.col = "green") dev.off() ### For no filter my.obj@main.data <- my.obj@raw.data # or (optinal example) #my.obj <- cell.filter(my.obj, # min.mito = 0, # max.mito = 0.10, # min.genes = 200, # max.genes = 8000, # min.umis = 0, # max.umis = Inf) ### cell cycle my.obj <- cc(my.obj, s.genes = s.phase, g2m.genes = g2m.phase) png("cellCycle.png") pie(table(my.obj@stats$Phase)) dev.off() ### run PCA on top 2000 genes my.obj <- run.pca(my.obj, top.rank = 2000) ### find best genes for second round PCA or batch alignment my.obj <- find.dim.genes(my.obj, dims = 1:30,top.pos = 20, top.neg = 20) length(my.obj@gene.model) ########### Batch alignment (CPCA method) my.obj <- iba(my.obj,dims = 1:30, k = 10, ba.method = "CPCA", method = "gene.model", gene.list = my.obj@gene.model) ### impute data my.obj <- run.impute(my.obj,dims = 1:10,data.type = "pca", nn = 10) ### run tSNE and UMAP my.obj <- run.pc.tsne(my.obj, dims = 1:10) my.obj <- run.umap(my.obj, dims = 1:10) ### run KNetL my.obj <- run.knetl(my.obj, dims = 1:20, k = 400, dim.redux = "umap") ### cluster based on KNetL coordinates # The object is already clustered but here is an example: # my.obj <- iclust(my.obj, k = 300, data.type = "knetl", dim.redux = "umap") ### save object save(my.obj, file = "my.obj.Robj") ### plot 1 A= cluster.plot(my.obj,plot.type = "pca",interactive = F,cell.size = 0.5,cell.transparency = 1, anno.clust=F) B= cluster.plot(my.obj,plot.type = "umap",interactive = F,cell.size = 0.5,cell.transparency = 1,anno.clust=F) C= cluster.plot(my.obj,plot.type = "tsne",interactive = F,cell.size = 0.5,cell.transparency = 1,anno.clust=F) D= cluster.plot(my.obj,plot.type = "knetl",interactive = F,cell.size = 0.5,cell.transparency = 1,anno.clust=F) library(gridExtra) png('AllClusters.png', width = 12, height = 10, units = 'in', res = 300) grid.arrange(A,B,C,D) dev.off() ### plot 2 png('AllCondsClusts.png', width = 15, height = 15, units = 'in', res = 300) cluster.plot(my.obj, cell.size = 0.5, plot.type = "knetl", cell.color = "black", back.col = "white", cell.transparency = 1, clust.dim = 2, interactive = F,cond.facet = T) dev.off() ### plot 3 genelist = c("LYZ","MS4A1","GNLY","FCGR3A","NKG7","CD14","S100A9","CD3E","CD8A","CD4","CD19","KLRB1","LTB","IL7R","GZMH","CD68","CCR7","CD68","CD69","CXCR4","IFITM3","IL32","JCHAIN","VCAN","PPBP") rm(list = ls(pattern="PL_")) for(i in genelist){ MyPlot <- gene.plot(my.obj, gene = i, interactive = F, cell.size = 0.1, plot.data.type = "knetl", data.type = "main", scaleValue = T, min.scale = -2.5,max.scale = 2.0, cell.transparency = 1) NameCol=paste("PL",i,sep="_") eval(call("<-", as.name(NameCol), MyPlot)) } library(cowplot) filenames <- ls(pattern="PL_") B <- cluster.plot(my.obj,plot.type = "knetl",interactive = F,cell.size = 0.1,cell.transparency = 1,anno.clust=T) filenames <- c("B",filenames) library(cowplot) png('genes_KNetL-plot.png',width = 18, height = 15, units = 'in', res = 300) plot_grid(plotlist=mget(filenames)) dev.off() ### plots for porportions pdf('clust_cond_freq_info_pie.pdf',width = 10, height = 10) clust.cond.info(my.obj, plot.type = "pie", normalize.ncell = TRUE, my.out.put = "plot", normalize.by = "percentage") dev.off() pdf('clust_cond_freq_info_bar.pdf',width = 10, height = 10) clust.cond.info(my.obj, plot.type = "bar", normalize.ncell = TRUE,my.out.put = "plot", normalize.by = "percentage") dev.off() my.obj <- clust.cond.info(my.obj, plot.type = "bar", normalize.ncell = TRUE) data <- (my.obj@my.freq) write.table(data, "myFrequencies.tsv", sep = "\t", row.names =F) pdf('clust_cond_freq_info_pie.cond.pdf',width = 10, height = 10) clust.cond.info(my.obj, plot.type = "pie.cond", normalize.ncell = TRUE, my.out.put = "plot", normalize.by = "percentage") dev.off() pdf('clust_cond_freq_info_bar.cond.pdf',width = 10, height = 10) clust.cond.info(my.obj, plot.type = "bar.cond", normalize.ncell = TRUE,my.out.put = "plot", normalize.by = "percentage") dev.off() ###### Cluster QC png('cluster_mito_ratio.png',width = 6, height = 6, units = 'in', res = 300) clust.stats.plot(my.obj, plot.type = "box.mito", interactive = F) dev.off() png('cluster_gene_cov.png',width = 6, height = 6, units = 'in', res = 300) clust.stats.plot(my.obj, plot.type = "box.gene", interactive = F) dev.off() ### cell cycle per cluster png("cluster_cc_pie.png",width = 6, height = 6, units = 'in', res = 300) clust.stats.plot(my.obj, plot.type = "pie.cc", interactive = F) dev.off() png("cluster_cc_bar.png",width = 6, height = 6, units = 'in', res = 300) clust.stats.plot(my.obj, plot.type = "bar.cc", interactive = F) dev.off() ### more plots library(gridExtra) A= cluster.plot(my.obj,plot.type = "pca",col.by = "conditions",interactive = F,cell.size = 0.5) B= cluster.plot(my.obj,plot.type = "umap",col.by = "conditions",interactive = F,cell.size = 0.5) C= cluster.plot(my.obj,plot.type = "tsne",col.by = "conditions",interactive = F,cell.size = 0.5) D= cluster.plot(my.obj,plot.type = "knetl",col.by = "conditions",interactive = F,cell.size = 0.5) png('AllConds.png', width = 8, height = 8, units = 'in', res = 300) grid.arrange(A,B,C,D) dev.off() cluster.plot(my.obj, cell.size = 2, plot.type = "knetl", cell.color = "black", back.col = "white", col.by = "clusters", cell.transparency = 1, clust.dim = 2, interactive = T, density = F, out.name = "2d_knetl_clusters")