单细胞数据复现-肺癌文章代码复现2

原创

小胡子刺猬的生信学习123

发布于 2022-05-07 12:16:21

2.1K60

代码可运行

文章被收录于专栏：文献分享及代码学习文献分享及代码学习

运行总次数：0

代码可运行

昨天是先对数据初步的质量进行过滤，今天是对过滤后的数据进行标准化和后面开始进行分群。其实还是比较流程的，但是还是有些东西是很细节的，是一些学习的片段的推文比不了的，主要是因为相对的是这是文章已经发表的文章，思路相对比较成熟。

数据标准化

一开始并不明白sct矫正和normalization标准化有什么区别，查了一些博主的文章，明白了以下几点：

SCT矫正选用的基因数目更多（3000个），大于normalization的2000个，相比而言，范围更多，同时有可能多选用的这1000个基因对下游分析会造成很大的影响。
SCT矫正不适宜做后面的DE（差异基因分析），主要是因为SCT矫正后面得到的是rnacount的残差，不是基本的count值，目前也有博主对两种方法的差异分析做了test，发现sct矩阵获得的结果相对较少；但是目前今年1月份seurat也推出了sct矩阵也可以做差异分析，首先需要在前面进行PrepSCTFindMarkers，但是可以发现目前的文章还有很多教程还是针对RNA矩阵做后面的findallmarkers，所以一般的建议是前面做normaliztion和SCT，后面做分析的时候，对其进行不同的矩阵的调取。
SCT矫正的时候是一行命令就做了nor里面的三个命令计算，但是相对于也比较费时间。

在作者的代码中，发现他把regress中加入rnacount和核糖体的基因来进行矫正，其实这部分是可以调整的，比较想去除线粒体基因的影响，可以放入线粒体的，这一步就是相对比较自己的研究内容来看。

# Data normalization；conserve.memory:节约资源

seu_obj <- SCTransform(seu_obj, verbose = T, vars.to.regress = c("nCount_RNA", "pMT"), conserve.memory = T)

saveRDS(seu_obj, file = "seurat_objects/all_SCTransform.RDS")

分群结果可视化

##pca降维，选取前面的50个维度对后面进行分析
seu_obj <- RunPCA(seu_obj)
ElbowPlot(seu_obj, ndims = 50)
##穷极循环，去看哪个dims比较适合后面的分群，一般是看到这些图的结果没有特别大的变化，表明这个dims可取
for (i in c(15, 20)) {
  umaptest <- RunUMAP(seu_obj, dims = 1:i, verbose = F)
  print(DimPlot(umaptest, reduction = "umap", group.by = "orig.ident") + labs(title = paste0(i, " dimensions")))
  print(FeaturePlot(umaptest, features = c("EPCAM", "PTPRC"), sort.cell = T))
  print(FeaturePlot(umaptest, features = c("MARCO", "KIT"), sort.cell = T))
  print(FeaturePlot(umaptest, features = c("FOXJ1", "AGER"), sort.cell = T))
  print(FeaturePlot(umaptest, features = c("JCHAIN", "VWF"), sort.cell = T))
  remove(umaptest)
}
#选取上面获得的适宜的dims做后面的umap的降维
seu_obj <- RunUMAP(seu_obj, dims = 1:15, verbose = T)
seu_obj <- FindNeighbors(seu_obj, dims = 1:15)
#穷极循环，分辨率test，一般是分辨率越高，分群越多，这里都没有对循环的结果进行保存，可以参照对marker基因进行可视化的保存图方式进行更改
for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  seu_obj <- FindClusters(seu_obj, resolution = i)
  print(DimPlot(seu_obj, reduction = "umap") + labs(title = paste0("resolution: ", i)))
}
##对分群的结果进行可视化
for (i in c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")) {
  print(FeaturePlot(seu_obj, features = i, coord.fixed = T, sort.cell = T))
}
##选用两种group的展示方式
DimPlot(seu_obj, group.by = "orig.ident")
DimPlot(seu_obj, group.by = "SCT_snn_res.0.5", label = T)

对部分marker基因进行可视化

挑选了一些比较的重要的marker基因进行可视化，这里发现分辨率是有选取不一样的分组的，跟刚刚的0.5是有差别的，怀疑是为了让分群的结果更少一些，有利于后期的绘图。

# Main cell type annotation,viridis(10):在颜色面板中选取10个颜色

mainmarkers <- c("PECAM1", "VWF", "ACTA2", "JCHAIN", "MS4A1", "PTPRC", "CD68", "KIT", "EPCAM", "CDH1", "KRT7", "KRT19")

for (i in seq_along(mainmarkers)) {
  FeaturePlot(seu_obj, features = mainmarkers[i], coord.fixed = T, order = T, cols = viridis(10))
  #ggsave2(paste0("FeaturePlot_mainmarkers_", mainmarkers[i], ".png"), path = "output/annotation", width = 10, height = 10, units = "cm")
}

DotPlot(seu_obj, features = mainmarkers, group.by = "SCT_snn_res.0.2") + 
  coord_flip() + 
  scale_color_viridis()
#ggsave2("DotPlot_mainmarkers.png", path = "output/annotation", width = 30, height = 8, units = "cm")

##label：标签
DimPlot(seu_obj, group.by = "SCT_snn_res.0.2", label = T, label.size = 5)
#ggsave2("DimPlot_all_clusters.png", path = "output/annotation", width = 20, height = 20, units = "cm")

人工分群注释后，加入注释结果

这里也是跟上面一样，选用的是0.2分辨率的分群结果，然后把注释的结果读到mata.data中。

##使用Idents（）函数可查看不同细胞的分群
Idents(seu_obj) <- seu_obj$SCT_snn_res.0.2
annotation_curated_main <- read_excel("../data/curated_annotation/curated_annotation_main.xlsx")
new_ids_main <- annotation_curated_main$main_cell_type
##levels是因子水平向量；RenameIdents ()函数 : 细胞簇注释名更改
names(new_ids_main) <- levels(seu_obj)
seu_obj <- RenameIdents(seu_obj, new_ids_main)
seu_obj@meta.data$main_cell_type <- Idents(seu_obj)

其他分组信息加入mata.data中

在一开始的时候，为了将每个样本的数据读进去，加载了一个excel表格，但是有其他的信息没用，所以这里是为了加载其他的信息进去，比如样本的来源，样本的取样时间。样本的分组等等，都是比较有利于后面的图形展示的分组用的，所以我觉得这一步很必要，可以用于后续不同的数据的矩阵图形展示。

# Add metadata
## fetchdata:抓取数据
metatable <- read_excel("../data/metadata/patients_metadata.xlsx")
metadata <- FetchData(seu_obj, "orig.ident")
metadata$cell_id <- rownames(metadata)
metadata$sample_id <- metadata$orig.ident
##left_join：左连接，保留X中得所有观测，这样可以保留原表数据
metadata <- left_join(x = metadata, y = metatable, by = "sample_id")
rownames(metadata) <- metadata$cell_id

##AddMetaData：添加元数据列；AddMetaData()中传递给元数据参数的内容必须是具有与数据@meta.data中的行名匹配的行名的数据框
seu_obj <- AddMetaData(seu_obj, metadata = metadata)

细胞周期蛋白计算

在文章中，作者是在后面对细胞周期蛋白进行了讨论的，但是一般都是在分析前面的时候，开始进行细胞周期蛋白的计算，也可以在sct矫正的时候加入细胞周期蛋白，去去除细胞周期蛋白的影响，这里主要的是考虑我们做的项目的物种来源，以及这部分的内容会不会会对后面的分析产生影响，其实还是需要不断的去尝试以及选择最适的分析方法。

# Cell cycle scoring

### add cell cycle, cc.genes loaded with Seurat

s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

score_cc <- function(seu_obj) {
  seu_obj <- CellCycleScoring(seu_obj, s.genes, g2m.genes)
  seu_obj@meta.data$CC.Diff <- seu_obj@meta.data$S.Score - seu_obj@meta.data$G2M.Score
  return(seu_obj)
}

seu_obj <- score_cc(seu_obj)

FeatureScatter(seu_obj, "G2M.Score", "S.Score", group.by = "Phase", pt.size = .1) +
  coord_fixed(ratio = 1)

下面是我在另一个博主那里摘取的代码。

if(T){
  g2m_genes <- cc.genes$g2m.genes
  g2m_genes <- CaseMatch(search=g2m_genes, match=rownames(seu_obj))
  s_genes <- cc.genes$s.genes    
  s_genes <- CaseMatch(search=s_genes, match=rownames(seu_obj))
  seu_obj <- CellCycleScoring(seu_obj, g2m.features=g2m_genes, s.features=s_genes)
  tmp <- RunPCA(seu_obj, features = c(g2m_genes, s_genes), verbose = F)
  p <- DimPlot(tmp, reduction = "pca", group.by = "orig.ident")
  ggsave("CellCycle_pca.png", p, width = 8, height = 6)
  rm(tmp)
}

亚群细分

作者将肺癌样本主要分为了以下的三个分群，因此根据注释结果，将每个亚群进行了提取，随后了对数据进行了ScaleData

处理，使得数据处于高斯分布状态。将主要提取的这三类亚群进行数据保存，后面的分析将单独对这几个亚群进行单独分析。

# Subset, rescale and save RDS files
###subset and rescale
saveRDS(seu_obj, file = "seurat_objects/all.RDS")
Idents(seu_obj) <- seu_obj@meta.data$main_cell_type
epi <- subset(seu_obj, idents = "Epithelial")
imm <- subset(seu_obj, idents = "Immune")
str <- subset(seu_obj, idents = "Stromal")
epi <- ScaleData(epi)
imm <- ScaleData(imm)
str <- ScaleData(str)
epi
imm
str
saveRDS(epi, file = "seurat_objects/epi.RDS")
saveRDS(imm, file = "seurat_objects/imm.RDS")
saveRDS(str, file = "seurat_objects/str.RDS")

fig1图片的复现

根据样本来源、病人来源等对ump的结果进行可视化，这里也是相对于比较流程。

# Plots for figure 1
###r plots for figure 1
DimPlot(seu_obj, group.by = "tissue_type", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1B.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(seu_obj, group.by = "patient_id", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1C.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(seu_obj, group.by = "main_cell_type", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1D_umap.png", path = "../results", width = 15, height = 15, units = "cm")

#mutate()总会添加新的列到已有列的末尾；rev ()函数用于返回数据对象的反向版本
cell_types <- FetchData(seu_obj, vars = c("sample_id", "main_cell_type", "tissue_type")) %>% 
  mutate(main_cell_type = factor(main_cell_type, levels = c("Stromal", "Immune", "Epithelial"))) %>% 
  mutate(sample_id = factor(sample_id, levels = rev(c("p018t", "p019t", "p023t", "p024t", "p027t", "p028t", "p030t", "p031t", "p032t", "p033t", "p034t", "p018n", "p019n", "p027n", "p028n", "p029n", "p030n", "p031n", "p032n", "p033n", "p034n"))))

##mapping：ggplot2绘图，代表映射；aes()函数将长格式的数据集映射到不同的图层
##coord_flip 横纵坐标位置转换
ggplot(data = cell_types) + 
  geom_bar(mapping = aes(x = sample_id, fill = main_cell_type, ), position = "fill", width = 0.75) +
  scale_fill_manual(values = use_colors) +
  coord_flip()
ggsave2("Fig1D_barplot.pdf", path = "../results", width = 15, height = 30, units = "cm")
##去除seu_obj这个环境变量
remove(seu_obj)