文献复现及简介—胰腺癌的单细胞水平肿瘤异质性

生信菜鸟团

发布于 2023-09-09 17:12:05

6301

❝群里的一个小伙伴提出来想让我帮忙复现cell里面的部分单细胞数据的图，今天我先来简单的梳理一下文献还有一部分复现代码到降维分群部分。同时我们群里的小伙伴会跟大家分享一下这篇文献，今晚八点，感兴趣的可以关注一下。 ❞

文献简介

「Microenvironment drives cell state, plasticity, and drug response in pancreatic cancer.」「Cell. 2021 Dec 9;184(25):6119-6137.e26.doi: 10.1016/j.cell.2021.11.017.」

文章中的数据集

https://singlecell.broadinstitute.org/single_cell/study/SCP1644

「除了Biopsy活检样本和Organoid类器官样本，文章中还有对公共数据库挖掘的数据。」https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131886

文章样品制备流程

「转移患者的样本通过核心针活检收集并分离，之后将活检细胞分配用于scRNA-seq，并通过下游串行scRNA-seq取样开发与患者匹配的类器官。」

「提出复现的小伙伴比较关心活检部分的复现，活检数据是23042个细胞，其中恶性7740，非恶性细胞15320个。」

「推断的拷贝数变异(CNV)改变评分分离了每个活检中的假定癌和非癌细胞，并与参考靶向DNA-seq显示了高度一致性。」

「23个病人的CNV score图」

「细胞亚群注释图」

部分代码复现

「step1:导入数据」

rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(data.table)
library(stringr)
getwd()
setwd("../")
###### step1:导入数据 ######  
#数据是csv格式的就常规读取即可。
ct=fread("./SCP1664_raw/Biopsy_RawDGE_23042cells.csv",data.table = F)
ct[1:4,1:4]
rownames(ct)=ct[,1]
ct=ct[,-1]
sce =CreateSeuratObject(counts =  ct ,
                        min.cells = 5,
                        min.features = 300 )

head(rownames(sce@meta.data))
tail(rownames(sce@meta.data))
phe=str_split(rownames(sce@meta.data),'_',simplify = T)
head(phe)
tail(phe)
sce@meta.data$orig.ident=paste0('',phe[,2])
table(sce@meta.data$orig.ident)
sce.all=sce
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
tail(sce.all@meta.data, 10)
table(sce.all$orig.ident)

「step2: QC省略」

「step3: 降维分群」

dir.create("2-harmony")
getwd()
setwd("2-harmony")

sce=sce.all 
sce
sce <- NormalizeData(sce, 
                         normalization.method = "LogNormalize",
                         scale.factor = 1e4) 
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))

library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T ) 
sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
                     dims = 1:15) 
sce.all=sce
#设置不同的分辨率，观察分群效果
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
  sce.all=FindClusters(sce.all, #graph.name = "CCA_snn",
                       resolution = res, algorithm = 1)
}
colnames(sce.all@meta.data)
apply(sce.all@meta.data[,grep("RNA_snn",colnames(sce.all@meta.data))],2,table)
head(colnames(sce.all))
table(sce.all$orig.ident)

#接下来分析
sel.clust = "RNA_snn_res.0.8"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
saveRDS(sce.all, "sce.all_int.rds")
getwd()
setwd('../')

「step4: 检查分群」

###### 检查常见分群情况  ######
dir.create("3-cell")
setwd("3-cell")  
getwd()

#22971
DimPlot(sce.all, reduction = "umap", group.by = "seurat_clusters",label = T) 
DimPlot(sce.all, reduction = "umap", group.by = "RNA_snn_res.0.8",label = T) 
ggsave('umap_by_RNA_snn_res.0.8.pdf',width = 8,height = 6)

sce.all=RunTSNE(sce.all,  dims = 1:15, 
                reduction = "harmony")
p_tsne=DimPlot(sce.all, reduction = "tsne", group.by = "RNA_snn_res.0.8",label = T);p_tsne
ggsave('tsne_by_RNA_snn_res.0.8.pdf',width = 8,height = 6)