单细胞实战之细胞周期矫正，双胞体和RNA污染去除——入门到进阶(初级篇2）

load('../phe.Rdata')
phe$group = sce$group = ifelse(sce$orig.ident %in% c("GSM5688706_WGC","GSM5688707_JCA","GSM5688708_LS-CRC3"),"left","right")
head(phe) 
gplots::balloonplot(
  table(phe$celltype,phe$orig.ident)
)

x = phe$orig.ident
y = phe$celltype
tbl =  table(x,y)
df = dcast(as.data.frame(tbl),x~y)
ptb = round(100*tbl/rowSums(df[,-1]),2)
df = as.data.frame(ptb)
colnames(df)=c('orig.ident','celltype','Percentage')
head(df)
gpinfo=unique(phe[,c('orig.ident','group')]);gpinfo
df=merge(df,gpinfo,by='orig.ident')
head(df)

p <- ggplot(df,aes(x=group,y=Percentage))+
  geom_boxplot(outlier.alpha=0, aes(fill=celltype))+facet_grid(~celltype,scales = "free")+
  geom_jitter(aes(x = group, y = Percentage,color=celltype))+ 
  scale_color_npg() +
  scale_fill_npg() +
  theme_bw()+
  theme(axis.text.x = element_text(face = "bold",angle = 45, hjust=1, vjust=0.5,size = 14), 
        legend.position= "none",
        strip.background = element_rect(color="black",fill = "steelblue3", linetype = "solid"),
        strip.text = element_text(face = "bold", color = "black",hjust = 0.5, size = 12,),
        plot.title=element_text(face="bold.italic",size="20", color="brown",hjust = 0.5),
        axis.text.y = element_text(face = "bold",size = 14) ,
        axis.title = element_text(size = 14),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) 
p   
pro='Percentage-of-each--celltype-in-orig.ident-by-group'
ggsave(paste0(pro,"_boxplot.pdf"),width = 12,height = 6)

rm(list=ls())
source('scRNA_scripts/lib.R')
library(DoubletFinder)
library(BiocParallel)
getwd()
dir.create('4-Corrective-data/')
setwd('4-Corrective-data/') 

register(MulticoreParam(workers = 4, progressbar = TRUE))
sce.all <- readRDS("../2-harmony/sce.all_int.rds")
sp='human' 
load('../phe.Rdata') 
identical(rownames(phe) , colnames(sce.all))
sce.all=sce.all[,colnames(sce.all) %in% rownames(phe) ]
sce.all@meta.data = phe[match(colnames(sce.all),rownames(phe) ),]
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
DimPlot(sce.all) 
colnames(sce.all@meta.data)

scRNA <- sce.all

# Seurat包自带了一组细胞周期标记基因，这些基因来自 Tirosh 等人在 2015 年发表的研究（通常是指用于单细胞 RNA 测序数据中的细胞周期分析）。
# 这些基因被用来识别细胞所处的细胞周期阶段，如G1、S、G2/M 期。
s.genes <- intersect(cc.genes$s.genes,rownames(scRNA))
g2m.genes <- intersect(cc.genes$g2m.genes,rownames(scRNA))
scRNA <- CellCycleScoring(scRNA,
                        s.features = s.genes,
                        g2m.features = g2m.genes,
                        set.ident = T)
colnames(scRNA@meta.data)
p1=VlnPlot(scRNA,features =  "S.Score" ,pt.size = 0,group.by = 'celltype') +NoLegend()
p2=VlnPlot(scRNA,features =  "G2M.Score" ,pt.size = 0,group.by = 'celltype')+NoLegend()
p1+p2
table(scRNA$Phase,scRNA$celltype)
  #     B cells endothelial cells epithelial/cancer cells fibroblasts mast cells myeloid cells plasma proliferative cells T/NK cells VSMCs
  # G1     2027               437                    1478         955        468          2763   1832                  68       5855   365
  # G2M     857                56                     637          35         78           336    164                 396       2557    14
  # S       920                56                     306          50         50           357    100                 339       2033    32

# features = rownames(scRNA)的情况需要前后一致
scRNA_new <- ScaleData(scRNA,
                       vars.to.regress = c("S.Score", "G2M.Score"),
                       features = rownames(scRNA))

# 查看不进行细胞周期和细胞之后处理之后的数据情况
# 分别查看celltype和细胞周期基因标记情况
p1=DimPlot(scRNA_new,reduction="pca",group.by='celltype') + DimPlot(scRNA,reduction="pca",group.by='celltype');p1
p2=DimPlot(scRNA_new,reduction="pca") + DimPlot(scRNA,reduction="pca") +NoLegend();p2
p1/p2

# 在这里我们重新对scRNA和scRNA_new数据进行PCA以及后续流程并进行结果对比
seurat <- RunPCA(scRNA,
                 features = VariableFeatures(object = scRNA),
                 verbose = FALSE)
seurat <- RunHarmony(seurat,"orig.ident")
seurat <- FindNeighbors(seurat, dims = 1:15, 
                     reduction = "harmony", 
                     verbose = FALSE) 
seurat <- FindClusters(seurat,resolution = 0.8, 
                    verbose = FALSE)
seurat <- RunUMAP(seurat,dims=1:15,reduction = "harmony",
                     umap.method = "uwot", metric = "cosine")
seurat <- DimPlot(seurat,reduction = 'umap',label = T);seurat

# scRNA_new
seuratObj <- RunPCA(scRNA_new,
                    features = VariableFeatures(object = scRNA_new),
                    verbose = FALSE)
seuratObj <- RunHarmony(seuratObj,"orig.ident")
seuratObj <- FindNeighbors(seuratObj, dims = 1:15, 
                     reduction = "harmony", 
                     verbose = FALSE) 
seuratObj <- FindClusters(seuratObj,resolution = 0.8, 
                    verbose = FALSE)
seuratObj <- RunUMAP(seuratObj,dims=1:15,reduction = "harmony",
                     umap.method = "uwot", metric = "cosine")
seuratObj <- DimPlot(seuratObj,reduction = 'umap',label = T);seuratObj

rm(list=ls())
source('scRNA_scripts/lib.R')
library(DoubletFinder)
library(BiocParallel)
getwd()
dir.create('4-Corrective-data/')
setwd('4-Corrective-data/') 

register(MulticoreParam(workers = 4, progressbar = TRUE))
sce.all <- readRDS("../2-harmony/sce.all_int.rds")
sp='human' 
load('../phe.Rdata') 
identical(rownames(phe) , colnames(sce.all))
sce.all=sce.all[,colnames(sce.all) %in% rownames(phe) ]
sce.all@meta.data = phe[match(colnames(sce.all),rownames(phe) ),]
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
DimPlot(sce.all) 
colnames(sce.all@meta.data)

# 定义主要参数和函数
run_doublet_finder <- function(sce_object, doublet_rate, pcs = 1:15, sct = FALSE) {
  # 参数扫描，选择最优 pK 值
  sweep.res <- paramSweep(sce_object, PCs = pcs, sct = sct)
  sweep.stats <- summarizeSweep(sweep.res, GT = FALSE)
  bcmvn <- find.pK(sweep.stats)
  pK_bcmvn <- bcmvn$pK[which.max(bcmvn$BCmetric)] %>% as.numeric()
  
  # 计算 homotypic doublets 和调整的 doublet 数目
  homotypic.prop <- modelHomotypic(sce_object$seurat_clusters)
  nExp_poi <- round(doublet_rate * ncol(sce_object))
  nExp_poi.adj <- round(nExp_poi * (1 - homotypic.prop))
  
  # 鉴定 doublets
  sce_object <- doubletFinder(sce_object, 
                              PCs = pcs, 
                              pN = 0.25, 
                              pK = pK_bcmvn, 
                              nExp = nExp_poi.adj, 
                              reuse.pANN = FALSE, 
                              sct = sct)
  return(sce_object)
}

# 运行 DoubletFinder 分析
sce_list <- SplitObject(scRNA, split.by = "orig.ident")
# check每个样本的细胞数
sce_list
# doublet_rate中的值需要根据
sce_list[["GSM5688706_WGC"]] <- run_doublet_finder(sce_list[["GSM5688706_WGC"]], 
                                                   doublet_rate = 0.018)

# 选择一个可视化一下
DimPlot(sce_list[["GSM5688706_WGC"]], 
        reduction = "umap", 
        group.by = "DF.classifications_0.25_24_59")

# 遍历每个子对象，直接修改列名
sce_list <- lapply(sce_list, function(sce) {
  # 修改 pANN 列名（匹配以 "pANN_0.25" 开头的列）
  pANN_col <- grep("^pANN_0.25", colnames(sce@meta.data), value = TRUE)
  if (length(pANN_col) > 0) {
    colnames(sce@meta.data)[colnames(sce@meta.data) %in% pANN_col] <- "pANN_0.25"
  }
  
  # 修改 DF.classifications 列名（匹配以 "DF.classifications_0.25" 开头的列）
  DF_col <- grep("^DF\\.classifications_0.25", colnames(sce@meta.data), value = TRUE)
  if (length(DF_col) > 0) {
    colnames(sce@meta.data)[colnames(sce@meta.data) %in% DF_col] <- "DF.classifications"
  }
  
  return(sce)
})

sce.all <- merge(sce_list[[1]], y= sce_list[ -1 ]) #add.cell.ids =  samples) 
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all <- JoinLayers(sce.all)
dim(sce.all[["RNA"]]$counts )
colnames(sce.all@meta.data)

sce.all <- NormalizeData(sce.all,normalization.method = "LogNormalize",
                         scale.factor = 1e4) 
sce.all <- FindVariableFeatures(sce.all)
sce.all <- ScaleData(sce.all)
sce.all <- RunPCA(sce.all, features = VariableFeatures(object = sce.all))
sce.all <- RunHarmony(sce.all, "orig.ident")
sce.all <- RunUMAP(sce.all,  dims = 1:15,reduction = "harmony")
DimPlot(sce.all,group.by = "DF.classifications")

# 如果需要过滤双细胞的话
sce.all.filtered <- subset(sce.all, subset = DF.classifications == "Singlet")
DimPlot(sce.all.filtered,group.by = "DF.classifications")

set.seed(123)
# 得到表达矩阵
counts <- GetAssayData(object = sce.all.filtered, slot = "counts")
decontX_res <- decontX(counts)
sce.all.filtered$contamination <- decontX_res$contamination
sce.all.filtered$contamination

# 可视化
library(ggplot2)
FeaturePlot(sce.all.filtered, 
            features = 'contamination', 
            raster=FALSE       # 细胞过多时候需要加这个参数
           ) + 
  scale_color_viridis_c()+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())+
  xlab('UMAP_1')+
  ylab('UMAP_2')

# contamination值在0-1之间,官方没有给出具体的参考值
# 有老师认为0.2可能是一个比较合适的界值
sce_filt <- sce.all.filtered[,sce.all.filtered$contamination < 0.2]
DimPlot(sce_filt,label = T)+NoLegend()
sce_filt
# An object of class Seurat 
# 21308 features across 24080 samples within 1 assay 
# Active assay: RNA (21308 features, 2000 variable features)
#  3 layers present: scale.data, data, counts
#  3 dimensional reductions calculated: pca, harmony, umap

qsave(sce_filt,"sce.all.qs")
setwd("..")