❝本来这周打算把上周欠下的NMF复现做一下,今天尝试了一下没有成功做出来,出了点小问题。 正好群里有小伙伴问能不能复现一下这篇文献,就来简单复现一下。 ❞
文献标题:
文献数据集:
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
getwd()
# setwd("../")
###### step1:导入数据 ######
library(stringr)
fs = list.files('./raw/',pattern = '^GSM')
#执行这一步需要解压tar -xvf
samples=str_split(fs,'_',simplify = T)[,1]
lapply(unique(samples),function(x){
y=fs[grepl(x,fs)]
folder=paste0("raw/", str_split(y[1],'_',simplify = T)[,1])
dir.create(folder,recursive = T)
#为每个样本创建子文件夹
file.rename(paste0("raw/",y[1]),file.path(folder,"barcodes.tsv.gz"))
#重命名文件,并移动到相应的子文件夹里
file.rename(paste0("raw/",y[2]),file.path(folder,"features.tsv.gz"))
file.rename(paste0("raw/",y[3]),file.path(folder,"matrix.mtx.gz"))
})
dir='./raw'
sceList = lapply(samples,function(pro){
print(pro)
sce =CreateSeuratObject(counts = Read10X(file.path(dir,pro)) ,
project = pro ,
min.cells = 5,
min.features = 300 )
return(sce)
})
sce.all=merge(x=sceList[[1]],
y=sceList[ -1 ],
add.cell.ids = gsub('^GSM[0-9]*_','',samples) )
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@meta.data$orig.ident)
#分组
sce=sce.all
sce$group<-ifelse(grepl("GSM5077732",sce$orig.ident),"Normal","Tumor")
table(sce$group)
dir.create("2-harmony")
getwd()
setwd("2-harmony")
# sce.all=readRDS("../1-QC/sce.all_qc.rds")
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)
head(colnames(sce.all))
table(sce.all$orig.ident)
table(sce.all$group)
#接下来分析,按照分辨率为0.5进行
sel.clust = "RNA_snn_res.0.5"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident)
saveRDS(sce.all, "sce.all_int.rds")
其中GSM5077732数据是normal,其余均为tumor
feats <- c("nFeature_RNA", "nCount_RNA","percent_mito", "percent_ribo", "percent_hb")
p1=VlnPlot(sce.all, features = feats, pt.size = 0.01, ncol = 2,group.by = "orig.ident") +
NoLegend()
p1
library(ggplot2)
ggsave(filename="Vlnplot.pdf",plot=p1,height = 8,width = 15)
getwd()
setwd('../')
###### step5:检查常见分群情况 ######
dir.create("3-cell")
setwd("3-cell")
getwd()
#sce.all=readRDS("./2-harmony/sce.all_int.rds")
DimPlot(sce.all, reduction = "umap", group.by = "seurat_clusters",label = T)
DimPlot(sce.all, reduction = "umap", group.by = "RNA_snn_res.0.5",label = T)
ggsave('umap_by_RNA_snn_res.0.5.pdf',width = 8,height = 6)
library(ggplot2)
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A','CD19', 'CD79A', 'MS4A1' ,
'IGHG1', 'MZB1', 'SDC1',
'CD68', 'CD163', 'CD14',
'TPSAB1' , 'TPSB2', # mast cells,
'MKI67','TOP2A','KLRC1',
'RCVRN','FPR1' , 'ITGAM' ,
'FGF7','MME', 'ACTA2',
'PECAM1', 'VWF',
'KLRB1','NCR1', # NK
'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1',
'MKI67' ,'TOP2A',"NKG7","S100A8","S100A9" )
library(stringr)
genes_to_check=str_to_upper(unique(genes_to_check))
genes_to_check
p_all_markers <- DotPlot(sce.all, features = genes_to_check,
assay='RNA' ) + coord_flip()
p_all_markers
p_all_markers
p_umap=DimPlot(sce.all, reduction = "umap",
group.by = "RNA_snn_res.0.5",label = T,label.box = T)
library(patchwork)
p_all_markers+p_umap
ggsave('markers_umap_3.pdf',width = 13,height = 6)
DimPlot(sce.all, reduction = "umap",
group.by = 'group',label = T)
ggsave('group_umap.pdf',width = 7,height = 6)
# 需要自行看图,定细胞亚群:
# 0,5:T
# 4:NK
# 1,12:B
# 6,16:cycling
# 7,10,11:macro
# 3:mono
# 2,8,9,13,15,17,19:epi
# 18:stromal
# 14:Endo
celltype=data.frame(ClusterID=0:19,
celltype= 0:19)
#定义细胞亚群
celltype[celltype$ClusterID %in% c( 1,12 ),2]='B'
celltype[celltype$ClusterID %in% c( 0,5),2]='T'
celltype[celltype$ClusterID %in% c( 4),2]='NK'
celltype[celltype$ClusterID %in% c(14),2]='Endo'
celltype[celltype$ClusterID %in% c( 6,16 ),2]='cycling'
celltype[celltype$ClusterID %in% c( 7,10,11),2]='macro'
celltype[celltype$ClusterID %in% c( 3),2]='mono'
celltype[celltype$ClusterID %in% c(2,8,9,13,15,17,19),2]='epi'
celltype[celltype$ClusterID %in% c(18),2]='stromal'
head(celltype)
celltype
table(celltype$celltype)
sce.all@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
sce.all@meta.data[which(sce.all@meta.data$RNA_snn_res.0.5 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(sce.all@meta.data$celltype)
th=theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
library(patchwork)
p_all_markers=DotPlot(sce.all, features = genes_to_check,
assay='RNA' ,group.by = 'celltype' ) + coord_flip()+th
p_umap=DimPlot(sce.all, reduction = "umap", group.by = "celltype",label = T,label.box = T)
p_all_markers+p_umap
ggsave('markers_umap_by_celltype_2.pdf',width = 13,height = 8)