小洁老师的小班课——单细胞专题Day1

不爱吃豆沙包

已于 2024-07-04 09:37:04 修改

阅读量417

点赞数 4

分类专栏：丘山的读书笔记文章标签：学习笔记 r语言

于 2024-07-04 09:36:13 首次发布

本文链接：https://blog.csdn.net/m0_71952366/article/details/140170620

版权

丘山的读书笔记专栏收录该内容

17 篇文章 7 订阅

订阅专栏

单细胞专题-Day1

第一件事当然是装包啦~

需要安装的R包：

cran包：'tidyverse','msigdbr', 'patchwork','SeuratObject','Seurat'

Bioconductor包： 'sva', 'monocle', 'GOplot','GSVA','plotmo','regplot','scRNAseq','BiocStyle','celldex','SingleR','BiocParallel'

装包超实用小函数：require 对于加载成功的包，返回逻辑值T，未装的包返回F，用于装包时代码的分情况讨论。

单细胞数据库

GEO
Single Cell Portal
Human Cell Atlas
Single Cell Expression Atlas
UCSC

不同类型文件的读取方法：https://mp.weixin.qq.com/s/W7szy-Kg6G1N1ENHNRjGiw 需要时查询即可。

多样本数据整理

跟着小洁老师，不用再手动整理文件夹！代码操作yyds！！

解压：

untar("GSE231920_RAW.tar",exdir = "GSE231920_RAW")#解包
library(stringr)
fs = paste0("GSE231920_RAW/",dir("GSE231920_RAW/"))
fs#查看每个样本的数据文件

从文件名中提取样本名

samples = dir("GSE231920_RAW/") %>% str_split_i(pattern = "_",i = 2) %>% unique();samples

不同样本需要修改代码，该例子为分割GSM7306054_sample1_barcodes.tsv.gz这一文件名

dir()列出工作目录下文件

为每个样本创建单独文件夹

##定义一个函数，创建01_data及下面的样本文件夹，s为样本名字
ctr = function(s){
  ns = paste0("01_data/",s)
  if(!file.exists(ns))dir.create(ns,recursive = T)
}
##使用lapply函数批量创建文件夹,ctr为刚创建的函数，samples为样本名的向量
lapply(samples,ctr)

每个样本的三个文件分别复制到文件夹

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]]))
    }
  }
})

更改所有文件名，使其符合读取要求

on = paste0("01_data/",dir("01_data/",recursive = T));on

nn = str_remove(on,"GSM\d+_sample\d_");nn

file.rename(on,nn)

d是正则表达式中的任意数字，可以借鉴AI的写法

批量读取

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是合并为一个表达矩阵
  save(sce.all,file = rdaf)   ##储存文件
}
load(rdaf)
head(sce.all@meta.data)

table(sce.all$orig.ident) #每个样本多少细胞

sum(table(Idents(sce.all)))#细胞总数

数据质控

去除离群值/双细胞等异常指标的细胞

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)

##画出小提琴图
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")

#根据小提琴图指定指标去掉离群值
sce.all = subset(sce.all,percent.mt < 20&
                   nCount_RNA < 40000 &
                   nFeature_RNA < 6000)
table(sce.all@meta.data$orig.ident)

PS：MAD法去除离群值

降维聚类分群

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") %>%        ##harmony去批次方法
    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)    #肘形图用来选择主成分数量

p1 =  DimPlot(sce.all, reduction = "umap",label = T,pt.size = 0.5)+ NoLegend();p1     ##查看umap图
DimPlot(sce.all, reduction = "umap",pt.size = 0.5,group.by = "orig.ident")     ##按样本着色

降维聚类分群中，FindClusteres函数的resolution可以调整，数值越大，分群越细，一般在0.3-1.2之间。

细胞类型注释

##手动注释
#FindAllmarkers函数
library(dplyr)
f = "markers.Rdata"
if(!file.exists(f)){
  markers <- FindAllMarkers(seu.obj, only.pos = TRUE,min.pct = 0.25)
  save(markers,file = f)
}
load(f)

a = read.table("markers_total.txt",sep = ",")
gt = split(a[,2],a[,1])

DotPlot(sce.all, features = gt,cols = "RdYlBu") +
  RotatedAxis()+p1

unique(a$V1)
writeLines(paste0(0:20,","))
celltype = read.table("anno.txt",sep = ",") 
celltype

levels(Idents(sce.all))
levels(sce.all)

new.cluster.ids <- celltype$V2
names(new.cluster.ids) <- levels(sce.all)
sce <- RenameIdents(sce.all, new.cluster.ids)
levels(Idents(sce))

sce <- subset(sce,idents = c('NULL'),invert=T)
levels(Idents(sce))

save(sce,file = "sce.Rdata")
p2 <- DimPlot(sce, reduction = "umap", 
              label = TRUE, pt.size = 0.5) + NoLegend()
p1+p2


##自动注释
library(celldex)
library(SingleR)
ls("package:celldex")

## [1] "BlueprintEncodeData"              "DatabaseImmuneCellExpressionData"
## [3] "HumanPrimaryCellAtlasData"        "ImmGenData"                      
## [5] "MonacoImmuneData"                 "MouseRNAseqData"                 
## [7] "NovershternHematopoieticData"

f = "ref_BlueprintEncode.RData"
#参考数据集，路径按需修改
if(!file.exists(f)){
  ref <- celldex::BlueprintEncodeData()
  save(ref,file = f)
}
ref <- get(load(f))
library(BiocParallel)
scRNA = sce.all
test = scRNA@assays$RNA$data
pred.scRNA <- SingleR(test = test, 
                      ref = ref,
                      labels = ref$label.main, 
                      clusters = scRNA@active.ident)
pred.scRNA$pruned.labels


new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
scRNA <- RenameIdents(scRNA,new.cluster.ids)
save(scRNA,file = "scRNA.Rdata")
p2 <- DimPlot(scRNA, reduction = "umap",label = T,pt.size = 0.5) + NoLegend()
p1+p2