1.下載并整理數(shù)據(jù)
演示數(shù)據(jù)：GSE231920厅篓，3個(gè)treat+3個(gè)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 = "_")分成多個(gè)部分(比如GSM7306054_sample1_barcodes.tsv.gz被分成GSM7306054遥巴、sample1传睹、barcodes.tsv.gz三個(gè)部分)媒抠，我們需要的是第2個(gè)部分（i=2）,即sample1可霎。
1.1字符串處理的函數(shù)
stringr包里的函數(shù)：
str_sub()：用于從字符串中提取子字符串科乎。可以根據(jù)起始位置和結(jié)束位置來選擇部分字符串前酿。
text <- "Hello, world!"
# 提取從第7個(gè)字符到最后的子字符串
result <- str_sub(text, start = 7)
print(result)
# 輸出結(jié)果為 "world!"
str_remove()：用于從字符串中移除指定的子字符串患雏。可以根據(jù)正則表達(dá)式或者普通字符串來指定需要移除的內(nèi)容罢维。
text <- "Data analysis is fun."
# 移除包含 "is" 的子字符串
result <- str_remove(text, "is")
print(result)
# 輸出結(jié)果為 "Data analysis fun."
str_detect()：用于檢測(cè)字符串中是否包含指定的模式或者子字符串淹仑。返回邏輯值（TRUE或FALSE），表示是否檢測(cè)到匹配
texts <- c("apple", "banana", "orange", "grape")
# 檢測(cè)包含 "an" 的字符串
result <- str_detect(texts, "an")
print(result)
# 輸出結(jié)果為 TRUE FALSE TRUE FALSE
str_replace()：用于替換字符串中的匹配內(nèi)容肺孵≡冉瑁可以根據(jù)正則表達(dá)式或者普通字符串來替換匹配到的內(nèi)容。
text <- "R programming is great!"
# 替換 "great" 為 "awesome"
result <- str_replace(text, "great", "awesome")
print(result)
# 輸出結(jié)果為 "R programming is awesome!"
1.2 文件處理函數(shù)
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] "單細(xì)胞.Rproj" "單細(xì)胞學(xué)小組geo數(shù)據(jù)" "山脊圖.png"
#[25] "未命名.R" "氣泡圖.png" "注釋.png"
#[28] "注釋umap.png" "注釋氣泡.png" "sce.all.Rdata"
dir(pattern = ".R$") #列出工作目錄下以.R結(jié)尾的文件 $表示以……結(jié)尾
# [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" "單細(xì)胞.Rproj"
##[10] "未命名.R"
file.create("douhua.txt") #用代碼創(chuàng)建文件
## [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這個(gè)函數(shù)里
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [[1]]
## [1] 1
##
## [[2]]
## [1] 2
##
## [[3]]
## [1] 3
##
## [[4]]
## [1] 4
1.4 自定義函數(shù)
比如分情況討論裝包可以寫成一個(gè)函數(shù)：
my_install = function(pkg){
if (!require(pkg,character.only = T))install.packages(pkg,update = F,ask = F)
}
my_install("tidyr")
#就會(huì)把”tidyr”代入到大括號(hào)里的pkg位置
#安裝多個(gè)包
ps = c("tidyr","dplyr","stringr")
lapply(ps, my_install)
1.5 整理成Read10X要求的格式
每個(gè)樣本單獨(dú)一個(gè)文件夾初婆，文件夾名是樣本名蓬坡，每個(gè)文件夾里3個(gè)固定名稱的文件，不可以有前綴磅叛。
# 為每個(gè)樣本創(chuàng)建單獨(dú)的文件夾屑咳，都放在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"
#每個(gè)樣本的三個(gè)文件復(fù)制到單獨(dú)的文件夾
#自定義函數(shù)可不起名字，直接放在lapply第二個(gè)參數(shù)的位置上
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
#輸出只留了前三個(gè)做示范
#[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
#輸出只留了前三個(gè)做示范
#[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
修改文件名時(shí)，路徑時(shí)要從工作目錄之下開始寫
2.批量讀取
rm(list = ls())
library(Seurat)
rdaf = "sce.all.Rdata"
if(!file.exists(rdaf)){
f = dir("01_data/")
scelist = list() #創(chuàng)建空的列表敲董，下面的for循環(huán)每執(zhí)行一次紫皇，scelist里面就會(huì)多一個(gè)元素。
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]]))#輸出每個(gè)文件的基因數(shù)和細(xì)胞數(shù)
}
sce.all = merge(scelist[[1]],scelist[-1]) #合并多個(gè)對(duì)象
sce.all = JoinLayers(sce.all)
#merge后腋寨，每個(gè)樣本的表達(dá)矩陣是一個(gè)獨(dú)立的的layer聪铺，JoinLayers是合并為一個(gè)表達(dá)矩陣
set.seed(335)
sce.all = subset(sce.all,downsample=700)#每個(gè)樣本抽700個(gè)細(xì)胞
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) #每個(gè)樣本多少細(xì)胞
##
## sample1 sample2 sample3 sample4 sample5 sample6
## 700 700 700 700 700 700
sum(table(Idents(sce.all)))#細(xì)胞總數(shù)
## [1] 4200
3.質(zhì)控指標(biāo)
除了線粒體基因(MT-)，核糖體(RP[SL])和紅細(xì)胞基因(HB[^(P)])也是常見的過濾指標(biāo)萄窜，也可以畫畫看铃剔，如果有離群值，也可以過濾掉查刻。
#過濾掉線粒體基因键兜，核糖體和紅細(xì)胞基因
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
#根據(jù)小提琴圖指定指標(biāo)去掉離群值
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普气，它需要的計(jì)算資源少，且準(zhǔn)確程度高佃延，是最受歡迎的方法现诀。
**需要?jiǎng)h掉原來的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"就會(huì)按照樣本來分配顏色于未。檢查一下去除樣本間批次效應(yīng)的效果如何，在umap圖上面看陡鹃，糊在一起烘浦，沒有每個(gè)樣本自成一簇，就說明效果還可以萍鲸。如果效果不太好闷叉，可以比較一下其他的去除批次效應(yīng)的方法(一般沒啥必要)。
