2022-11-02 单细胞数据整合

自用学习

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

x=list.files()

dir = c('BC2/', "BC21/")
names(dir) = c('BC2',  'BC21')      

#读取数据
scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)
}

#标准化 寻找高变基因
for (i in 1:length(scRNAlist)) {
  scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]])
  scRNAlist[[i]] <- FindVariableFeatures(scRNAlist[[i]], selection.method = "vst",nfeatures = 3000)
}

#设置并行运算
library(future)
plan("multiprocess", workers =3)
options(future.globals.maxSize = 2000 * 1024^2)
#寻找锚定点
?FindIntegrationAnchors
scRNA.anchors <- FindIntegrationAnchors(object.list = scRNAlist,anchor.features = 2000)

scRNA1 <- IntegrateData(anchorset = scRNA.anchors)

整合后的assay有两个 一个是RNA 有29738个基因，一个是integrated 有2000个基因 默认的assay（activate.assay）是aintegrated

2000个锚定点 用于整合数据 而非画图 后期作图要将activate.assay改为RNA

#对数据进行scale
scRNA1=ScaleData(scRNA1)
scale后的数据存放在integrated里面

RNA里面的scale data为0 ，因此，在使用seurat锚定点方法时，要做两次scale data

scRNA1 <- RunPCA(scRNA1, npcs = 30, verbose = T)
# t-SNE and Clustering

scRNA1 <- FindNeighbors(scRNA1, reduction = "pca", dims = 1:20)
scRNA1 <- FindClusters(scRNA1, resolution = 0.8)
scRNA1 <- RunUMAP(scRNA1, reduction = "pca", dims = 1:20)
colnames(scRNA1@meta.data)
scRNA1 <- RunTSNE(scRNA1, dims = 1:20)
colnames(scRNA1@meta.data)
DimPlot(scRNA1, reduction = "umap", group.by = "orig.ident")
DimPlot(scRNA1, reduction = "umap", label = TRUE)

DefaultAssay(scRNA1) <- "RNA"
scRNA <- ScaleData(scRNA1)

Harmony

install.packages("devtools")
library(devtools)
install_github("immunogenomics/harmony")
library(harmony)
BiocManager::install("SingleCellExperiment")


scRNA.2=Read10X("BC2/")
scRNA.21=Read10X("BC21/")
scRNA.21 = CreateSeuratObject(scRNA.21 ,project="sample_21",min.cells = 3, min.features = 200)
scRNA.2= CreateSeuratObject(scRNA.2 ,project="sample_2",min.cells = 3, min.features = 200)

scRNA_harmony <- merge(scRNA.21, y=c(scRNA.2 ))
scRNA_harmony <- NormalizeData(scRNA_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)

system.time({scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")})
###问题 harmony之后的数据在哪里？


###一定要指定harmony
scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)

scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:16)
?DimPlot
plot1 =DimPlot(scRNA_harmony, reduction = "umap",label = T) 
plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident') 
#combinate
plotc <- plot1+plot2
plotc