1.下载并整理数据
演示数据：GSE231920，3个treat+3个control
untar("GSE231920_RAW.tar",exdir = "GSE231920_RAW")#解包
library(stringr)
fs = paste0("GSE231920_RAW/",dir("GSE231920_RAW/"))
samples = dir("GSE231920_RAW/") %>% str_split_i(pattern = "_",i = 2) %>% unique();samples
## 从文件名中提取到样本名，用str_split_i把文件名称按照"_"作为分隔符(pattern = "_")分成多个部分(比如GSM7306054_sample1_barcodes.tsv.gz被分成GSM7306054、sample1、barcodes.tsv.gz三个部分)，我们需要的是第2个部分（i=2）,即sample1。
1.1字符串处理的函数
stringr包里的函数：
str_sub()：用于从字符串中提取子字符串。可以根据起始位置和结束位置来选择部分字符串。
text <- "Hello, world!"
# 提取从第7个字符到最后的子字符串
result <- str_sub(text, start = 7)
print(result)
# 输出结果为 "world!"
str_remove()：用于从字符串中移除指定的子字符串。可以根据正则表达式或者普通字符串来指定需要移除的内容。
text <- "Data analysis is fun."
# 移除包含 "is" 的子字符串
result <- str_remove(text, "is")
print(result)
# 输出结果为 "Data analysis fun."
str_detect()：用于检测字符串中是否包含指定的模式或者子字符串。返回逻辑值（TRUE或FALSE），表示是否检测到匹配
texts <- c("apple", "banana", "orange", "grape")
# 检测包含 "an" 的字符串
result <- str_detect(texts, "an")
print(result)
# 输出结果为 TRUE FALSE TRUE FALSE
str_replace()：用于替换字符串中的匹配内容。可以根据正则表达式或者普通字符串来替换匹配到的内容。
text <- "R programming is great!"
# 替换 "great" 为 "awesome"
result <- str_replace(text, "great", "awesome")
print(result)
# 输出结果为 "R programming is awesome!"
1.2 文件处理函数
dir() # 列出工作目录下的文件
# [1] "2024-06-18-151153.jpg" "2024-06-18-151153.png" "anno.txt"
# [4] "featureplot.png" "GSE231920_RAW" "GSE231920_RAW.tar"
# [7] "heatmap.png" "input" "markers.Rdata"
#[10] "my_markers.txt" "obj.Rdata" "plot_zoom.png"
#[13] "plot_zoom2.png" "plot_zoom3.png" "plot_zoom4.png"
#[16] "ref_BlueprintEncode.RData" "ref_Human_all.RData" "sce.all.Rdata"
#[19] "sce.Rdata" "single_ref" "violin.png"
#[22] "单细胞.Rproj" "单细胞学小组geo数据" "山脊图.png"
#[25] "未命名.R" "气泡图.png" "注释.png"
#[28] "注释umap.png" "注释气泡.png" "sce.all.Rdata"
dir(pattern = ".R$") #列出工作目录下以.R结尾的文件 $表示以……结尾
# [1] "未命名.R"
dir(pattern = ".R") #列出工作目录下包含.R的文件
##[1] "GSE231920_RAW" "GSE231920_RAW.tar" "markers.Rdata"
##[4] "obj.Rdata" "ref_BlueprintEncode.RData" "ref_Human_all.RData"
## [7] "sce.all.Rdata" "sce.Rdata" "单细胞.Rproj"
##[10] "未命名.R"
file.create("douhua.txt") #用代码创建文件
## [1] TRUE
file.exists("douhua.txt") #某文件在工作目录下是否存在
## [1] TRUE
file.remove("douhua.txt") #用代码删除文件
## [1] TRUE
file.exists("douhua.txt") #删掉了就不存在
## [1] FALSE
## 可以批量的新建和删除
f = paste0("douhua",1:10,".txt")
file.create(f)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
file.remove(f)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
1.3 lapply
lapply(1:4, print) #就是把1、2、3、4依次代入print这个函数里
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [[1]]
## [1] 1
##
## [[2]]
## [1] 2
##
## [[3]]
## [1] 3
##
## [[4]]
## [1] 4
1.4 自定义函数
比如分情况讨论装包可以写成一个函数：
my_install = function(pkg){
if (!require(pkg,character.only = T))install.packages(pkg,update = F,ask = F)
}
my_install("tidyr")
#就会把”tidyr”代入到大括号里的pkg位置
#安装多个包
ps = c("tidyr","dplyr","stringr")
lapply(ps, my_install)
1.5 整理成Read10X要求的格式
每个样本单独一个文件夹，文件夹名是样本名，每个文件夹里3个固定名称的文件，不可以有前缀。
# 为每个样本创建单独的文件夹，都放在01_data下面
ctr = function(s){
ns = paste0("01_data/",s)
if(!file.exists(ns))dir.create(ns,recursive = T)
}
lapply(samples, ctr)
list.files(path = "01_data")
## [1] "sample1" "sample2" "sample3" "sample4" "sample5" "sample6"
#每个样本的三个文件复制到单独的文件夹
#自定义函数可不起名字，直接放在lapply第二个参数的位置上
lapply(fs, function(s){
#s = fs[1]
for(i in 1:length(samples)){
#i = 1
if(str_detect(s,samples[[i]])){
file.copy(s,paste0("01_data/",samples[[i]]))
}
}
})
image.png
#所有文件改名，去掉前缀。
on = paste0("01_data/",dir("01_data/",recursive = T));on
#输出只留了前三个做示范
#[1] "01_data/sample1/GSM7306054_sample1_barcodes.tsv.gz"
# [2] "01_data/sample1/GSM7306054_sample1_features.tsv.gz"
# [3] "01_data/sample1/GSM7306054_sample1_matrix.mtx.gz"
nn = str_remove(on,"GSM\\d+_sample\\d_");nn
#输出只留了前三个做示范
#[1] "01_data/sample1/barcodes.tsv.gz" "01_data/sample1/features.tsv.gz"
#[3] "01_data/sample1/matrix.mtx.gz"
file.rename(on,nn)
image.png
修改文件名时，路径时要从工作目录之下开始写
2.批量读取
rm(list = ls())
library(Seurat)
rdaf = "sce.all.Rdata"
if(!file.exists(rdaf)){
f = dir("01_data/")
scelist = list() #创建空的列表，下面的for循环每执行一次，scelist里面就会多一个元素。
for(i in 1:length(f)){
pda <- Read10X(paste0("01_data/",f[[i]]))
scelist[[i]] <- CreateSeuratObject(counts = pda,
project = f[[i]],
min.cells = 3,
min.features = 200)
print(dim(scelist[[i]]))#输出每个文件的基因数和细胞数
}
sce.all = merge(scelist[[1]],scelist[-1]) #合并多个对象
sce.all = JoinLayers(sce.all)
#merge后，每个样本的表达矩阵是一个独立的的layer，JoinLayers是合并为一个表达矩阵
set.seed(335)
sce.all = subset(sce.all,downsample=700)#每个样本抽700个细胞
save(sce.all,file = rdaf)
}
load(rdaf)
head(sce.all@meta.data)
## orig.ident nCount_RNA nFeature_RNA
## AAAGAACAGTCTGTAC-1_1 sample1 3884 1377
## AAAGGGCTCTCGCGTT-1_1 sample1 1615 721
## AAAGGTACAATTGGTC-1_1 sample1 3806 1215
## AAAGTGAGTATCGAAA-1_1 sample1 5554 1456
## AAAGTGATCTCAACCC-1_1 sample1 3111 1328
## AACAAAGGTAAGGTCG-1_1 sample1 3658 1243
table(sce.all$orig.ident) #每个样本多少细胞
##
## sample1 sample2 sample3 sample4 sample5 sample6
## 700 700 700 700 700 700
sum(table(Idents(sce.all)))#细胞总数
## [1] 4200
3.质控指标
除了线粒体基因(MT-)，核糖体(RP[SL])和红细胞基因(HB[^(P)])也是常见的过滤指标，也可以画画看，如果有离群值，也可以过滤掉。
#过滤掉线粒体基因，核糖体和红细胞基因
sce.all[["percent.mt"]] <- PercentageFeatureSet(sce.all, pattern = "^MT-")
sce.all[["percent.rp"]] <- PercentageFeatureSet(sce.all, pattern = "^RP[SL]")
sce.all[["percent.hb"]] <- PercentageFeatureSet(sce.all, pattern = "^HB[^(P)]")
head(sce.all@meta.data, 3)
## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rp percent.hb
## AAAGAACAGTCTGTAC-1_1 sample1 3884 1377 4.763131 23.09475 0
## AAAGGGCTCTCGCGTT-1_1 sample1 1615 721 0.619195 35.85139 0
## AAAGGTACAATTGGTC-1_1 sample1 3806 1215 7.777194 37.12559 0
#过滤掉离群值，画小提琴
VlnPlot(sce.all,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt",
"percent.rp",
"percent.hb"),
ncol = 3,pt.size = 0.5, group.by = "orig.ident")
image.png
#根据小提琴图指定指标去掉离群值
sce.all = subset(sce.all,percent.mt < 20&
nCount_RNA < 40000 &
nFeature_RNA < 6000)
table(sce.all@meta.data$orig.ident)
##
## sample1 sample2 sample3 sample4 sample5 sample6
## 687 622 683 686 678 675
4.整合降维聚类分群
多样本的整合，使用harmony，它需要的计算资源少，且准确程度高，是最受欢迎的方法。
**需要删掉原来的obj.Rdata！！！！！
f = "obj.Rdata"
library(harmony)
if(!file.exists(f)){
sce.all = sce.all %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData(features = rownames(.)) %>%
RunPCA(pc.genes = VariableFeatures(.)) %>%
RunHarmony("orig.ident") %>%
FindNeighbors(dims = 1:15, reduction = "harmony") %>%
FindClusters(resolution = 0.5) %>%
RunUMAP(dims = 1:15,reduction = "harmony") %>%
RunTSNE(dims = 1:15,reduction = "harmony")
save(sce.all,file = f)
}
load(f)
ElbowPlot(sce.all)
image.png
p1 = DimPlot(sce.all, reduction = "umap",label = T,pt.size = 0.5)+ NoLegend();p1
image.png
DimPlot(sce.all, reduction = "umap",pt.size = 0.5,group.by = "orig.ident")
image.png
group.by = "orig.ident"就会按照样本来分配颜色。检查一下去除样本间批次效应的效果如何，在umap图上面看，糊在一起，没有每个样本自成一簇，就说明效果还可以。如果效果不太好，可以比较一下其他的去除批次效应的方法(一般没啥必要)。
