单细胞Day7-多样本数据-1

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

https://satijalab.org/seurat/articles/seurat5_integration