实战探究五个参数对UMAP图可视化的影响

rm(list=ls())
options(stringsAsFactors = F) 
source('../scRNA_scripts/lib.R')
getwd()

###### step1: 导入数据 ######   
# 付费环节 800 元人民币
# 参考：https://mp.weixin.qq.com/s/tw7lygmGDAbpzMTx57VvFw
dir='../inputs/'
samples=list.files( dir )
samples = samples[str_detect(samples,"singlecell_count_matrix.txt")]
sceList = lapply(samples,function(pro){ 
  # pro=samples[7] 
  print(pro)  
  ct <- data.table::fread( file.path(dir,pro),
                           data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  #ct=t(ct) 
  sce =CreateSeuratObject(counts =  ct ,
                          project =  gsub('_singlecell_count_matrix.txt.gz','',pro)  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
}) 
do.call(rbind,lapply(sceList, dim))
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids = samples  ) 
names(sce.all@assays$RNA@layers)
sce.all[["RNA"]]$counts 
# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all <- JoinLayers(sce.all)
dim(sce.all[["RNA"]]$counts )

as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident) 

library(stringr)
sce.all$orig.ident=str_split(colnames(sce.all),'[-_]',simplify = T)[,1]
table(sce.all$orig.ident) 
sce.all

# 如果为了控制代码复杂度和行数 
# 可以省略了质量控制环节
###### step2: QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')
sp='human'
getwd()

sce.all.filt <- NormalizeData(sce.all.filt, 
                           normalization.method = "LogNormalize",
                           scale.factor = 1e4) 

## 评估参数的影响
variable_nums <- c(2000,3000,4000)
npcs_nums <- c(30,50,70,90,110)
neighbors_nums <- c(20,30,40,50)
min_dist_nums <- c(0.001,0.01,0.1,0.3,0.5) 
# https://zhuanlan.zhihu.com/p/352461768
# RunUMAP: n.neighbors default is 30[5-50]; min.dist default is 0.3[0.001-0.5]
if(F){
variable_nums=2000
npcs_nums=50
neighbors_nums=30
min_dist_nums=0.001
}

for (i in variable_nums){
  for (j in npcs_nums){
    tmp_sce <- FindVariableFeatures(sce.all.filt, selection.method = "vst",
                                    nfeatures = i)
    tmp_sce <- ScaleData(tmp_sce)
    tmp_sce <- RunPCA(tmp_sce, features = VariableFeatures(object = tmp_sce), 
                      verbose = FALSE, npcs= j)
    tmp_sce <- RunHarmony(tmp_sce, "orig.ident")
    tmp_sce <- FindNeighbors(tmp_sce, dims = 1:10, verbose = FALSE,reduction = "harmony")
    tmp_sce <- FindClusters(tmp_sce, resolution = 0.5, verbose = FALSE)
    table(tmp_sce$seurat_clusters)
    
    for (n in neighbors_nums){
      for (d in min_dist_nums){
        tmp_sce <- RunUMAP(tmp_sce, dims = 1:15, umap.method = "uwot", metric = "cosine", reduction = "harmony", n.neighbors=n, min.dist=d)
        
        pdf(paste0("pca_object_v",i,".pca_",j,"neighbors_",n,"dist_",d,".pdf"))
        p1 <- DimPlot(tmp_sce,label = T,reduction = "umap",raster=F)
        print(p1)
        genes_to_check = c("CD3E","CD3D","MS4A1","CD79A")
        T_genes = c("CD3E","CD3D")
        B_genes = c("MS4A1","CD79A")
        p2 <- DotPlot(tmp_sce , features = genes_to_check )   +
        coord_flip()  +
        theme(axis.text.x=element_text(angle=45,hjust = 1))
        print(p2)
        p3 <- FeaturePlot(object = tmp_sce, features = T_genes,raster=F) 
        print(p3)
        p4 <- FeaturePlot(object = tmp_sce, features = B_genes,raster=F) 
        print(p4)        
        dev.off()
      }
    }
  }
}