做海量单细胞测序，大概大部分的人都会从seraut开始，功能有很多，主要是获得matrix后的数据可视化，在分析的过程中，细胞类型注释是比较l考验操作人员的问题，具体的方法有很多，我在文中做了简单的介绍，其中对于操作者而言最简单的方式大概使用软件进行自动化注释了，所以也对singleR进行了一个测试，供大家参考~~
1 首先先加载以下我们需要的包
suppressMessages({
library(reshape2)
library(Seurat)
library(corrplot)
library(dplyr)
library(grid)
library(cowplot)
library(grid)
library(tidyverse)})
2 加载完包之后读入我们需要的数据
data1 <- Read10X('./Marrow_V2_1_matrix_10X') #读取10X数据
data2 <- Read10X('./Marrow_V2_2_matrix_10X')
#分析公共数据的时候遇到了都是csv的文件，其中barcode还两列，尝试了各种方式，最后用的
#x = Seurat::ReadMtx(mtx = "MTX文件", cells = "barcode文件",features = "gene文件",mtx.transpose = TRUE,feature.sep = ",",cell.sep = ",")#还有参数可以指定用哪一些来做barcode和gene.不过这个函数好像只有Seurat4可以用
3 读入数据后将数据整合成可用于seraut分析的格式,并且将这2个数据进行合并
PRO1 <- CreateSeuratObject(counts = data1,project = "Bong_marrow1", min.cells = 5) #将数据整理成10X可分析的对象
PRO2 <- CreateSeuratObject(counts = data2,project = "Bong_marrow2", min.cells = 5)
PRO<- merge(x=PRO1,y=PRO2) #将2个样本合并
4 统计数据的线粒体的含量，如果有需要，也可以统计数据的核糖体基因数据
PRO[["percent.mt"]] <- PercentageFeatureSet(PRO, pattern = "^MT-") #统计样本中的线粒体含量比例，如果线粒体统计出来都是0，可以试试把MT改成mt，可能是大小写的原因
#PRO[["percent.ribo"]] <- PercentageFeatureSet(PRO, pattern = "^Rp[sl]") #看核糖体RNA的含量
如果想可视化结果可以
VlnPlot(PRO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),group.by='orig.ident', pt.size = 0,ncol = 3)#这一步也可以用group.by的参数设置展示分组，一般会选择按照样本展示,如果想要有点，可以把pt.size = 0设置成=1.
image.png
5 接下来我们对数据进行以下处理，做均一化和数据缩放，降维分析等，可以看一下umap的降维图
PRO <- NormalizeData(object = PRO) #数据均一化，由于在测序过程中存在测序偏好性的问题，可能每个细胞的测序并不完全一致，所以Seraut对其进行均一化，默认的方法是log1p，即每个基因测得的数目除以该细胞的总的reads数目然后乘以10000，然后加1后取log值，通过这种方法，测序深度不同的细胞也可以放在一起分析，当然了，测序深度的差异最好不要太大。
PRO <- ScaleData(object = PRO)#可以理解为数据缩放，很多时候当数据表达量过大的时候比如1和10000，在图上就很难展现了，所以这个时候将数据缩放，就可以把绘图的差距拉小，也方便后边的降维分析，一般数据缩放采用z-score的方法
PRO <- FindVariableFeatures(object = PRO) #选择高可变基因#genes.use<- head(HVFInfo(object = PRO),2000)
PRO <- RunPCA(object=PRO,features = VariableFeatures(object = PRO))
PRO <- FindNeighbors(PRO, reduction = "pca", dims = 1:20) #首先计算每个细胞的KNN，也就是计算每个细胞之间的相互距离，依据细胞之间邻居的overlap来构建snn graph。
PRO <- FindClusters(PRO,resolution = 0.8, algorithm = 1) #根据SNN的方法来分群，具体的算法看不明白，就先了解是干嘛的吧
PRO <- RunTSNE(object=PRO,dims.use=1:20,do.fast=TRUE,check_duplicates = FALSE)#用tsne的方法进行降维
PRO <- RunUMAP(PRO, reduction = "pca", dims = 1:20)#用umap的方法进行降维
DimPlot(PRO, reduction = "umap", label = TRUE)
如果想保存图片的话
#P<-DimPlot(object = PRO, reduction = "umap",label = TRUE,repel = TRUE)
#print(P)
#ggplot2::ggsave(filename = "umap.pdf")) # 图片尺寸可以调整width = 13, height = 10
image
6 细胞类型注释
在获得这样一个分群图后，那我们可能关注，每个cluster是什么细胞类型，注释的原理是看每个细胞的marker基因的表达的情况，常用的方法有3种，一种是直接看每个cluster中高表达哪些基因，这些基因属于那种细胞类型高表达的基因。 第二种方法是用我们已知的marker基因看fetureplot 第三种就是自动化注释，现在有很多的自动化软件也很成熟，大家也可以试试，在这里主要介绍手动注释的2种方法，说实话，没有数据库的帮助，手动注释还是很难和消耗时间的。
第一种方法
subc<-levels(x=PRO@active.ident) #找到每个cluster的marker，然后将这些marker写到excel表里边，后期可以根据每个cluster里的的marker进行注释
for (l in subc){
cluster.markers <- FindMarkers(object=PRO,ident.1=l,min.pct=0.1,logfc.threshold = 0.25)
write.table(cluster.markers,file=paste(l,'_diffgenes.xls',sep=''),sep='\t',quote=F,row.names=T)}
dev.off()
#cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)当然也可以找针对某个cluster的marker
#luster1.markers <- FindMarkers(pbmc, ident.1 = 1,ident.2=2, min.pct = 0.25) #可以看cluster1和cluster2的差异基因，如果想看指定基因的差异，可以用features 这个参数指定基因。
第二种方法，以Ms4a1和Cd79a是B细胞marker，Cd3d,Trac是T细胞marker为例
FeaturePlot(PRO, features = c("Ms4a1", "Cd3d", ))
第三种方法指的是自动化注释，常用的自动化注释软件有SingleR，Cellassign，Celaref，scTPA等等，后期有时间也会做个注释测评。
image
7 细胞类型重命名
一般第一种方式和第二种方式都是需要一些文献积累的的，或者也可以直接用数据库，注释后我们则可以将其进行标记
#Idents(PRO)<-'RNA_snn_res.0.8'
cluster <- c("Neutrophils", "Neutrophils","Monocytes","Neutrophils","Monocytes","GMP","GMP","pDCs","ProEryth","T cells","HSCs","GMP","cDCs","B cells","Plasma Cells","Neutrophils","ProEryth","Basophils","Unknown","Monocytes","Neutrophils","Pre_B","Unknown","Unknown","pDCs","Erythroblasts","Monocytes","B cells","Unknown","pDCs","Unknown")
names(cluster) <- levels(PRO)
PRO <- RenameIdents(PRO, cluster)
PRO[["cluster"]]<- Idents(object = PRO)#将注释信息加入到meta.data中
PRO1<-subset(PRO,idents = 'Unknown', invert = TRUE) #有一些cluster比较难注释，为了不妨碍后边的分析，我们这里先把它过滤掉。
DimPlot(PRO, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()#这一步也可以通过group.by看不同分组的umap图。
image
8细胞亚型细分
很多时候我们会对我们感兴趣的细胞亚型进行细分，细分后还需要映射回原来的cluster，在这里也记录一下方法
PRO1<-subset(PRO,idents = 'T cells')
PRO1 <- FindNeighbors(PRO1, reduction = "pca", dims = 1:20) #首先计算每个细胞的KNN，也就是计算每个细胞之间的相互距离，依据细胞之间邻居的overlap来构建snn graph。
PRO1 <- FindClusters(PRO1,resolution = 0.3, algorithm = 1) #根据SNN的方法来分群，具体的算法看不明白，就先了解是干嘛的吧
PRO1 <- RunTSNE(object=PRO1,dims.use=1:20,do.fast=TRUE,check_duplicates = FALSE)#用tsne的方法进行降维
PRO1 <- RunUMAP(PRO1, reduction = "pca", dims = 1:20)#用umap的方法进行降维
DimPlot(PRO1) #可以看一下拉出来的分群
image.png
对细分的群进行重命名
cluster <- c("T1", "T2")
names(cluster) <- levels(PRO1)
PRO1 <- RenameIdents(PRO1, cluster)
PRO1[["cluster1"]]<- Idents(object = PRO1)#将注释信息加入到meta.data中
DimPlot(PRO1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
image.png
将重新注释的cluster映射回去
meta <- PRO@meta.data
meta$cell_name <- rownames(meta)
meta$cluster <- as.character(PRO@active.ident)
ident <- meta$cluster
names(ident) <- rownames(meta)
meta1 <- PRO1@meta.data
ident1 <- as.character(PRO1@active.ident)
names(ident1) <- names(PRO1@active.ident)
for(i in names(ident1)){
if( i %in% meta$cell_name){
ident[i] <- ident1[i]}
}
meta$cluster <- ident
PRO@meta.data <- meta
Idents(PRO)<-'cluster'
DimPlot(PRO,label=TRUE)
image.png
如图，T细胞被分为了T1，T2 2个cluster
9 除了分群以外，大家可能还关注每种细胞所占的比例
cellnum <- table(PRO@active.ident, PRO@meta.data[,"orig.ident"]) #统计每个样本种每种细胞亚型的细胞数目
freq<- prop.table(x=table(PRO@active.ident ,PRO@meta.data[,"orig.ident"]),margin=2) #统计每个每个样本中每个细胞亚型的比例
head（freq）#每种细胞类型占的比例
image
测试2 SingleR
1 首先加载需要用到的包
suppressMessages({
library(Seurat)
library(tidyverse)
library(SingleR)
library(celldex)
})
2 下载需要用的数据库
#获取需要的数据库
Mouse.se=MouseRNAseqData() #在这里我们使用小鼠的数据库，由于下载比较慢，所以也可以第一次加载后保存下来
#hpca.se=HumanPrimaryCellAtlasData() ##第一次载入会下载数据集，可能会慢一些，后面在用时就不用下载了
#Blue.se=BlueprintEncodeData()
#Immune.se=DatabaseImmuneCellExpressionData()
#Nover.se=NovershternHematopoieticData()
#MonacoIm.se=MonacoImmuneData()
#ImmGen.se=ImmGenData() #(鼠)
#Mouse.se=MouseRNAseqData() #(鼠)
saveRDS(Mouse.se,'Mouse.se.rds')
3 用singleR进行注释
PRO_for_SingleR <- GetAssayData(PRO, slot="data") ##获取标准化矩阵
RO.mouse <- SingleR(test = PRO_for_SingleR, ref = Mouse.se, labels = Mouse.se$label.main) # 使用小鼠参考数据集，main大类注释，也可使用fine小类注释，不过准确性不好确定
PRO@meta.data$labels <-PRO.mouse$labels #加标签
p<-DimPlot(PRO, group.by = c("cluster", "labels"),reduction = "umap",label=TRUE) #画图
ggplot2::ggsave(p,filename = "singleRumap.pdf",width = 15, height = 10)
image.png
第一张图是人工注释的，第二张图是singleR注释的，Emmmmm，自动化还是有点问题的，在大小群的定义上可能并不如人那么智能，不过像这种有发育关系的组织类型，可能也的确对自动化注释不太友好，当然也可以用别的数据库试试，结果仅供参考吧
关于初步的分析大部分内容是这些，如果大家对某个基因在不同的cluster里边的表达感兴趣可以用小提琴图VlnPlot，多个基因在不同cluster里边的表达可以用DotPlot和DoHeatmap来画，有什么小的需求都可以通过（？函数名称，比如？DotPlot）的方式看看参数，其实我们大部分的需求，写代码的的人都想到了，我们可以直接通过函数的参数来解决，有问题评论区见~~
网友评论