Using plot1cell to draw better s

2023-07-05  本文已影响0人  MYS_bio_man

1. Preparations before using plot1cell [or packages installing]

# turbo/fasten your R while connecting to network
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")))


dev.packages <- c("chris-mcginnis-ucsf/DoubletFinder","Novartis/hdf5r","mojaveazure/loomR")
devtools::install_github(dev.packages)

bioc.packages <- c("biomaRt","GenomeInfoDb","EnsDb.Hsapiens.v86","GEOquery","simplifyEnrichment","ComplexHeatmap")
BiocManager::install(bioc.packages)

devtools::install_github("TheHumphreysLab/plot1cell")
## or the development version, devtools::install_github("HaojiaWu/plot1cell")

2. Load packages and example data

library(plot1cell)
iri.integrated <- Install.example() # example data was a seurat object
# too slow
# searching "GSE139107" in geo and download yourself
# then reading data
# the folowed scripts are packaged in function "Install.example"
getGEOSuppFiles(GEO = "GSE139107") # not run/passed while download yourself 
setwd("GSE139107/") # not run/passed while download yourself 
all_files <- list.files(pattern = "dge")
all_ct <- list()
for (i in 1:length(all_files)) {
    ct_data <- read.delim(all_files[i])
    ct_data <- Matrix(as.matrix(ct_data), sparse = T)
    all_ct[[i]] <- ct_data
    print(all_files[i])
}
all.count <- RowMergeSparseMatrices(all_ct[[1]], all_ct[-1])
meta.data <- read.delim("GSE139107_MouseIRI.metadata.txt.gz")
all.count <- all.count[, rownames(meta.data)]
iri <- CreateSeuratObject(counts = all.count, min.cells = 0,
                          min.features = 0, meta.data = meta.data)
iri.list <- SplitObject(iri, split.by = "orig.ident")
iri.list <- lapply(X = iri.list, FUN = function(x) {
    x <- NormalizeData(x, verbose = FALSE)
    x <- FindVariableFeatures(x, verbose = FALSE)
})
features <- SelectIntegrationFeatures(object.list = iri.list)
iri.list <- lapply(X = iri.list, FUN = function(x) {
    x <- ScaleData(x, features = features, verbose = FALSE)
    x <- RunPCA(x, features = features, verbose = FALSE)
})
anchors <- FindIntegrationAnchors(object.list = iri.list,
                                  reference = c(1, 2), reduction = "rpca", dims = 1:50) ## Fast integration using reciprocal PCA (RPCA)
iri.integrated <- IntegrateData(anchorset = anchors, dims = 1:50)
iri.integrated <- ScaleData(iri.integrated, verbose = FALSE)
iri.integrated <- RunPCA(iri.integrated, verbose = FALSE)
iri.integrated <- RunUMAP(iri.integrated, dims = 1:25, min.dist = 0.2)
iri.integrated <- SetIdent(iri.integrated, value = "celltype")
levels(iri.integrated) <- c("PTS1", "PTS2", "PTS3", "NewPT1",
                            "NewPT2", "DTL-ATL", "MTAL", "CTAL1", "CTAL2", "MD",
                            "DCT", "DCT-CNT", "CNT", "PC1", "PC2", "ICA", "ICB",
                            "Uro", "Pod", "PEC", "EC1", "EC2", "Fib", "Per", "Mø",
                            "Tcell")
levels(iri.integrated@meta.data$Group) <- c("Control", "4hours",
                                            "12hours", "2days", "14days", "6weeks")
DefaultAssay(iri.integrated) <- "RNA"

colnames(iri.integrated@meta.data) # check colnames
head(iri.integrated) # check data

3. Using plot_circlize to draw a circle plot

circ_data <- prepare_circlize_data(iri.integrated, scale = 0.8 ) # scaled for plot,[to make it bigger or smaller], has no relation to your data
set.seed(666666) # sed 666666 seems wonderfull
# 设置细胞分群信息的颜色
cluster_colors<-rand_color(length(levels(iri.integrated)))
group_colors<-rand_color(length(names(table(iri.integrated$Group))))
rep_colors<-rand_color(length(names(table(iri.integrated$orig.ident))))

### plot and save figures
# 绘制细胞分群圈图
pdf('circlize_plot.pdf', width = 6, height = 6)
plot_circlize(circ_data,do.label = T, pt.size = 0.1, col.use = cluster_colors ,bg.color = 'white', kde2d.n = 200, repel = T, label.cex = 1)
# 添加细胞群注释信息
add_track(circ_data, group = "Group", colors = group_colors, track_num = 2) ## can change it to one of the columns in the meta data of your seurat object
add_track(circ_data, group = "orig.ident",colors = rep_colors, track_num = 3) ## can change it to one of the columns in the meta data of your seurat object
dev.off()
image.png

4. Using complex_dotplot_single to draw a dot plots for marker genes

png(filename =  'dotplot_single.png', width = 4, height = 6,units = 'in', res = 100) # pdf or png, two formats you shold learn
complex_dotplot_single(seu_obj = iri.integrated, feature = "Havcr1",groups = "Group")
dev.off()
image.png

设置groups和splitby参数对多个分组信息进行分割绘图

iri.integrated@meta.data$Phase<-plyr::mapvalues(iri.integrated@meta.data$Group, from = levels(iri.integrated@meta.data$Group), to = c("Healthy",rep("Injury",3), rep("Recovery",2)))
iri.integrated@meta.data$Phase<-as.character(iri.integrated@meta.data$Phase)
png(filename =  'dotplot_single_split.png', width = 4, height = 6,units = 'in', res = 100)
complex_dotplot_single(iri.integrated, feature = "Havcr1",groups = "Group",splitby = "Phase")
dev.off()
image.png

To visualize the same gene on multiple group factors, simply add more group factor IDs to the groups argument.

png(filename =  'dotplot_more_groups.png', width = 8, height = 6,units = 'in', res = 100)
complex_dotplot_single(seu_obj = iri.integrated, feature = "Havcr1",groups= c("Group","Replicates"))
dev.off()
image.png

Each group factor can be further splitted by its own factor if the splitby argument is provided. Note that in this case, the order of the group factors needs to match the order of splitby factors.

iri.integrated@meta.data$ReplicateID<-plyr::mapvalues(iri.integrated@meta.data$Replicates, from = names(table((iri.integrated@meta.data$Replicates))), to = c(rep("Rep1",3),rep("Rep2",3), rep("Rep3",1)))
iri.integrated@meta.data$ReplicateID<-as.character(iri.integrated@meta.data$ReplicateID)

png(filename =  'dotplot_more_groups_split.png', width = 9, height = 6,units = 'in', res = 200)
complex_dotplot_single(seu_obj = iri.integrated, feature = "Havcr1",groups= c("Group","Replicates"), splitby = c("Phase","ReplicateID"))
dev.off()
### In this example, "Phase" is a splitby factor for "Group" and "ReplicateID" is a splitby factor for "Replicates".
image.png

complex_dotplot_multiple can be used too

To visualize multiple genes in dotplot format, complex_dotplot_multiple should be used.

png(filename =  'dotplot_multiple.png', width = 10, height = 4,units = 'in', res = 300)
complex_dotplot_multiple(seu_obj = iri.integrated, features = c("Slc34a1","Slc7a13","Havcr1","Krt20","Vcam1"),group = "Group", celltypes = c("PTS1" ,   "PTS2"  ,  "PTS3"  ,  "NewPT1" , "NewPT2"))
dev.off()
image.png

5. Violin plot to show gene expression across groups

png(filename =  'vlnplot_single.png', width = 4, height = 6,units = 'in', res = 100)
complex_vlnplot_single(iri.integrated, feature = "Havcr1", groups = "Group",celltypes   = c("PTS1" ,   "PTS2"  ,  "PTS3"  ,  "NewPT1" , "NewPT2"))
dev.off()
image.png

Similar to complex_dotplot_single, the complex_vlnplot_single function also allows splitting the group factor by another factor with the argument splitby. [i've to say, this is the best fuc for me]

png(filename =  'vlnplot_single_split.png', width = 4, height = 6,units = 'in', res = 100)
complex_vlnplot_single(iri.integrated, feature = "Havcr1", groups = "Group",celltypes   = c("PTS1" ,   "PTS2"  ,  "PTS3"  ,  "NewPT1" , "NewPT2"), splitby = "Phase")
dev.off()
image.png

One gene/multiple group factors violin plot[one gene per multi groups and for more than one cell types]:

png(filename =  'vlnplot_multiple.png', width = 6, height = 6,units = 'in', res = 100)
complex_vlnplot_single(iri.integrated, feature = "Havcr1", groups = c("Group","Replicates"),celltypes   = c("PTS1" ,   "PTS2"  ,  "PTS3"  ,  "NewPT1" , "NewPT2"), font.size = 10)
dev.off()
image.png

If we add another gene "PRMT5"

png(filename =  'vlnplot_multiple.png', width = 6, height = 6,units = 'in', res = 100)
complex_vlnplot_single(iri.integrated, feature = c("Havcr1","PRMT5"), groups = c("Group","Replicates"),celltypes   = c("PTS1" ,   "PTS2"  ,  "PTS3"  ,  "NewPT1" , "NewPT2"), font.size = 10)
dev.off()

image.png

[to be continue...]

上一篇下一篇

猜你喜欢

热点阅读