將創(chuàng)建對(duì)象的過(guò)程記錄

setwd("~/pyc/07-cellranger/01-Fudan_EEC_Normal/")

# library(scater)
library(dplyr)
library(Seurat)
library(patchwork)
library(SingleCellExperiment)

library(ComplexHeatmap)
library(ConsensusClusterPlus)
library(msigdbr)
library(fgsea)
library(dplyr)
# library(tibble)
library(DoubletFinder)
library(Signac)
library(ggplot2)
library(stringr)
library(SingleR)
library(psych)
library(clustree)
library(patchwork)
library(ggplot2)
library(reshape2)
library(RColorBrewer)
library(paletteer)
library(harmony)
library(RColorBrewer) 
library(viridis)
library(wesanderson)


##### 01-1-設(shè)置文件讀取目錄  #####
## 批量讀取數(shù)據(jù)
### 設(shè)置數(shù)據(jù)路徑與樣本名稱(chēng)
pwd <- getwd()
assays <- dir("./")
# (dir <- paste0(pwd, "/", assays,"/filtered_feature_bc_matrix/"))
(dir <- paste0(pwd, "/", assays,"/","outs/filtered_feature_bc_matrix/"))
# 按文件順序給樣本命名寞焙，名稱(chēng)不要以數(shù)字開(kāi)頭隧膏，中間不能有空格 
samples_name = c('T_223', 'T_224', 'N_227', 'T_229', 'N_230', 'N_231',
                 
                 'T_866', 'T_869', 'N_870', 'N_871')

##### 01-2-創(chuàng)建單細(xì)胞對(duì)象  #####
scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  #不設(shè)置min.cells過(guò)濾基因會(huì)導(dǎo)致CellCycleScoring報(bào)錯(cuò)：
  #Insufficient data values to produce 24 bins.  
  scRNAlist[[i]] <- CreateSeuratObject(counts, project=samples_name[i],
                                       min.cells=3, min.features = 200)
  #給細(xì)胞barcode加個(gè)前綴，防止合并后barcode重名
  scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])   
  #計(jì)算線(xiàn)粒體基因比例
  if(T){    
    scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-") 
  }
  #計(jì)算核糖體基因比例
  if(T){
    scRNAlist[[i]][["percent.rb"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^RP[SL]")
  }
  #計(jì)算紅細(xì)胞基因比例
  if(T){
    HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
    HB.genes <- CaseMatch(HB.genes, rownames(scRNAlist[[i]]))
    scRNAlist[[i]][["percent.HB"]]<-PercentageFeatureSet(scRNAlist[[i]], features=HB.genes) 
  }
}

### 給列表命名并保存數(shù)據(jù)
dir.create("../02_Integrate")
setwd("../02_Integrate/")

names(scRNAlist) <- samples_name
#system.time(save(scRNAlist, file = "Integrate/scRNAlist0.Rdata")) 
system.time(saveRDS(scRNAlist, file = "scRNAlist0.rds")) 
# "2023-11-08 15:07:35 CST" 在11月8日磁滚，到這一步。




#### 現(xiàn)在是2024-4-7，已經(jīng)過(guò)去好久了届吁！ ####
pwd <- getwd()
assays <- dir("./")
# (dir <- paste0(pwd, "/", assays,"/filtered_feature_bc_matrix/"))
(dir <- paste0(pwd, "/", assays,"/","filtered_feature_bc_matrix/"))
(dir <- dir[1:5])
# 按文件順序給樣本命名，名稱(chēng)不要以數(shù)字開(kāi)頭绿鸣，中間不能有空格 
samples_name = c('A_225', 'A_226', 'A_228', 'A_867', 'A_868')


scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  #不設(shè)置min.cells過(guò)濾基因會(huì)導(dǎo)致CellCycleScoring報(bào)錯(cuò)：
  #Insufficient data values to produce 24 bins.  
  scRNAlist[[i]] <- CreateSeuratObject(counts, project=samples_name[i],
                                       min.cells=3, min.features = 200)
  #給細(xì)胞barcode加個(gè)前綴疚沐，防止合并后barcode重名
  scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])   
  #計(jì)算線(xiàn)粒體基因比例
  if(T){    
    scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-") 
  }
  #計(jì)算核糖體基因比例
  if(T){
    scRNAlist[[i]][["percent.rb"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^RP[SL]")
  }
  #計(jì)算紅細(xì)胞基因比例
  if(T){
    HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
    HB.genes <- CaseMatch(HB.genes, rownames(scRNAlist[[i]]))
    scRNAlist[[i]][["percent.HB"]]<-PercentageFeatureSet(scRNAlist[[i]], features=HB.genes) 
  }
}
names(scRNAlist) <- samples_name

scRNAlist1 <- c(scRNAlist0, scRNAlist)

# 最后merge
##### 01-3-Merge #####
scRNA <- merge(scRNAlist1[[1]], scRNAlist1[2:length(scRNAlist1)])
scRNA 
# An object of class Seurat 
# 43055 features across 142654 samples within 1 assay 
# Active assay: RNA (43055 features, 0 variable features)


table(scRNA$orig.ident)
# A_225 A_226 A_228 A_867 A_868 N_227 N_230 N_231 N_870 N_871 T_223 T_224 T_229 T_866 T_869 
# 7570 11054  6112 10677  8117  8126  8001 11993  7518  8281  9914  4876 12560 13611 14244 
saveRDS(scRNA,file = 'scRNA_orig.rds')

# file = "/home/newdisk/pyc/07-cellranger/01-Fudan_EEC_Normal/scRNA_orig.rds"

rm(scRNAlist1)


# 拆分seurat對(duì)象
ifnb.list <- SplitObject(scRNA, split.by = "orig.ident")

最后的結(jié)果

比較有用的代碼就是

Seurat V4版本拆分對(duì)象

參考

Seurat 4.0 | 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)整合(scRNA-seq integration)
多個(gè)10x單細(xì)胞對(duì)象的合并和批次校正--seurat錨點(diǎn)整合+Harmony