Seurat单细胞分析常见代码-01
2021-11-19 本文已影响0人
whitebird
1.如何通过Seurat对象获取细胞对应的坐标?
UMAP_coord <- data_obj@reductions$umap@cell.embeddings
TSNE_coord <- data_obj@reductions$tsne@cell.embeddings
pbmc_umap_coord <- data_obj@reductions$umap
pbmc_umap_coord2 <- as.data.frame(pbmc_umap_coord)
2.如何通过Seurat对象获取原始表达数据?
DefaultAssay(data_obj) <- 'RNA'
countdata <- GetAssayData(data_obj, slot = "counts")
write.table(countdata,"./finial_result/raw_count_combine.txt")
# GetAssayData(data_obj, slot = "counts")等同于data_obj@assays$RNA@counts
GetAssayData(data_obj, slot = "counts")
data_obj@assays$RNA@counts
3.提取特定基因阳性表达的细胞？
# 单个基因
DefaultAssay(data_obj) <- "RNA"
gene_name <- "CD4"
names(which(data_obj@assays$RNA@counts[gene_name,]>0))
WhichCells(object = data_obj, expression = CD4 > 0, slot = "counts")
# 也可用FetchData()提取基因表达量
# for total co-positive cells:
Co_Positive <- function(object,gene1,gene2){
CellData <- FetchData(object, vars = c(gene1,gene2), slot = "data")
print(head(CellData))
Gene1_total_cells <- sum(CellData[,gene1]>0)
cat("Total cells positive for", gene1,Gene1_total_cells," ")
Gene2_total_cells <- sum(CellData[,gene2]>0)
cat("Total cells positive for", gene2,Gene2_total_cells)
Co_Postitive_cell_ID <- subset(CellData, CellData[,gene1]>0)
Co_Postitive_cell_ID <- subset(Co_Postitive_cell_ID, Co_Postitive_cell_ID[,gene2]>0)
print(head(Co_Postitive_cell_ID))
Co_Postitive_cell_Total <- nrow(Co_Postitive_cell_ID)
cat("Total cells positive for",gene1," & ",gene2," ", Co_Postitive_cell_Total)
}
# 从总样本总细胞中提取出CD4 CXCR5 Foxp3三阳性的细胞
WhichCells(data_obj, slot = 'counts', expression = CD4 > 0 | CXCR5 > 0 | FOXP3 > 0)
How to subset using AND, working on raw counts as above.
WhichCells(data_obj, slot = 'counts', expression = CD4 > 0 & CXCR5 > 0 & FOXP3 > 0)
# 提取细胞
select_cells <- WhichCells(data_obj, slot = 'counts', expression = CD4 > 0 & CXCR5 > 0 & FOXP3 > 0)
select_obj <- subset(data_obj, cells = select_cells)
4.修改tsne/umap图的配色
.cluster_cols <- c(
"#DC050C", "#FB8072", "#1965B0", "#7BAFDE", "#882E72",
"#B17BA6", "#FF7F00", "#FDB462", "#E7298A", "#E78AC3",
"#33A02C", "#B2DF8A", "#55A1B1", "#8DD3C7", "#A6761D",
"#E6AB02", "#7570B3", "#BEAED4", "#666666", "#999999",
"#aa8282", "#d4b7b7", "#8600bf", "#ba5ce3", "#808000",
"#aeae5c", "#1e90ff", "#00bfff", "#56ff0d", "#ffff00")
p1 <- DimPlot(data_obj, reduction = "umap",label=T,group.by = "cluster",
cols = colorRampPalette(.cluster_cols,space="rgb")(length(unique(ColonData$cluster))))
save_plot("DimPlot_UMAP_clustering.png", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL, dpi=600)
save_plot("DimPlot_UMAP_clustering.pdf", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL)
5.查看tsne/umap图的默认配色
# 如celltype最后一个为Unkown，用灰色表示
library(scales)
my_color_palette <- c(hue_pal()(3), "gray")
Idents(data_obj) <- "celltype"
p1 <- DimPlot(data_obj, reduction = "tsne", label = TRUE, label.size = 3, pt.size=0.1, raster=FALSE) + scale_color_manual(values = my_color_palette)+labs(title="")
6.加载R包
方式1：
使用suppressMessages运行的时候不显示提示信息
suppressMessages({
library(Seurat)
library(cowplot)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(RColorBrewer)
library(reshape2)
})
方式2：
有的包在加载的时候会显示很多信息，可以使用 suppressPackageStartupMessages 函数来抑制这些信息的输出。
invisible(lapply(c("R.utils","Seurat","dplyr","kableExtra","ggplot2","scater",
"scran","BiocSingular","Matrix","cowplot"), function(x) {
suppressPackageStartupMessages(library(x,character.only = T))
}))
方式3：
easypackages能轻松加载多个R包
# Load libraries
library(easypackages)
packages <- c('clusterProfiler', 'org.Mm.eg.db', 'org.Hs.eg.db', 'GSEABase', 'biomaRt', 'enrichplot','plyr', 'dplyr', 'cowplot', 'ggplot2', 'patchwork')
libraries(packages)
7.给cluster定义细胞类型
data_obj$celltype <- "NA"
data_obj@meta.data[data_obj@meta.data$cluster %in% c("B0", "B2", "B6"),]$celltype <- "Naïve B"
data_obj@meta.data[data_obj@meta.data$cluster %in% c("B1","B3","B4","B5"),]$celltype <- "Memory B"
data_obj@meta.data[data_obj@meta.data$cluster %in% c("B7"),]$celltype <- "Plasma B"
data_obj@meta.data[data_obj@meta.data$cluster %in% c("B8"),]$celltype <- "Unknown"
8.小鼠基因转成首字母大写
library(Hmisc)
gene_list$Gene <- capitalize(tolower(gene_list$Gene))
9.统计不同resolution下的cluster数目分布
cluster_list <- sapply(grep("RNA_snn_res",sort(colnames(data_obj@meta.data)),value = TRUE),function(x) length(unique(data_obj@meta.data[,x])))
cluster_list <- as.data.frame(cluster_list) %>% tibble::rownames_to_column("resolution")
colnames(cluster_list) <- c("resolution", "cluster_nums")
write.csv(cluster_list, file = file.path(sample_analysis_path, "cluster_resolution_range_table.csv"), row.names = F)
10.可视化不同resolution下的umap图
DiffResolution_plot <- lapply(seq(0.05,0.3,by=0.05),function(i){
p <- DimPlot(data_obj, reduction = "umap",label=T,group.by = paste("RNA_snn_res.",i,sep=""),cols = .cluster_cols)+
labs(title=paste("Resolution = ",i,sep=""))+NoLegend()
return(p)
})
p_res <- cowplot::plot_grid(plotlist = DiffResolution_plot,ncol=3)
save_plot("Cluster_with_Different_resolution_UMAP.png", p_res, base_height = 12, base_aspect_ratio = NULL, base_width = 20, dpi=600)
save_plot("Cluster_with_Different_resolution_UMAP.pdf", p_res, base_height = 12, base_aspect_ratio = NULL, base_width = 20)
11.指定resolution为当前cluster
Idents(data_obj) <- selected_resolution
data_obj$cluster <- factor(Idents(data_obj), levels = sort(as.numeric(levels(Idents(data_obj)))))
12.统计每个样本的celltype占比
total_populations <- data_obj@meta.data %>% group_by(celltype) %>% summarize (total.pop = n())
# gets the proportion of cells for each cell type within a group by dividing by the total
count_populations <- data_obj@meta.data %>% group_by_at(vars(orig.ident, "celltype")) %>% summarize (n = n())
count_populations <- left_join(count_populations, total_populations, by = "celltype")
count_populations <- count_populations %>% mutate (proportion = n/total.pop)
count_populations <- count_populations %>% arrange(celltype, orig.ident)
colnames(count_populations) <- c("sample", "celltype", "count", "total.pop", "proportion")
write.csv(count_populations, "celltype_proportion_table.csv")
