生新技能树单细胞GBM数据分析(SignleR以及Seurat 联合分析及细胞簇注释

学习是一种态度

图片来自网络

关于单细胞测序分析，本文主要参考生新技能树团队的帖子和代码，有部分内容属于自己的理解，在此非常感谢生新技能树团队无私的奉献。当然本帖子也参考了大量其他的贴子，参考内容和链接将会在相应地方放出，以便读者学习。
单细胞分析技术目前尚多，SignleR和Seruat是主流的分析，之前我们就Seruat 官网上的流程进行了初步的分析，这次的内容主要是结合SingleR 和Seruat 进行单细胞的注释。分析过程主要包括

1. 数据的读取和整理

数据的下载和读取
rm(list=ls())##清空环境变量
### 1、下载、探索、整理数据----
## 1.1 下载、探索数据
#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465
sessionInfo()
# 读取文件耗时比较长，请耐心等待
a <- read.table("/home/data/vip38/Data_resource/Single_Cell_Data_Test/GSE84465_GBM_All_data.csv.gz")##直接读取下载的数据压缩包
View(a)  ###得到的a是一个包含有23,465个features 和 3589个细胞的矩阵

#行名为symbol ID
#列名为sample，看上去像是两个元素的组合。
summary(a[,1:4]) ##对a 的数学相关特征(最大，最小，中位数，上下四分位数，平均值进行探索)
boxplot(a[,1:4])
head(rownames(a))
tail(rownames(a),10)

# 可以看到原文的counts矩阵来源于htseq这个计数软件，所以有一些不是基因的行需要剔除：
#  "no_feature"           "ambiguous"            "too_low_aQual"        "not_aligned"          "alignment_not_unique"
tail(a[,1:4],10)

a1=a[1:(nrow(a)-5),]


#原始counts数据

####样本相关的信息
#3,589 cells of 4 human primary GBM samples, accession number GSE84465
#2,343 cells from tumor cores and 1,246 cells from peripheral regions
b <- read.table("/home/data/vip38/Data_resource/Single_Cell_Data_Test/SraRunTable.txt",
                sep = ",", header = T)
b[1:4,1:4]
table(b$Patient_ID) # 4 human primary GBM samples
table(b$TISSUE) # tumor cores and peripheral regions
table(b$TISSUE,b$Patient_ID)###查看每个样本的肿瘤核心和外周组织的细胞数目


## 1.2 整理数据 
# tumor and peripheral 分组信息
head(colnames(a))###对细胞的命名进行探索，发现，其由plate的ID 和Well 的ID 组合而成
head(b$plate_id)
head(b$Well)
View(b)
#a矩阵列名（sample）并非为GSM编号，而主要是由相应的plate_id与Well组合而成

b.group <- b[,c("plate_id","Well","TISSUE","Patient_ID")]
paste0("X",b.group$plate_id[1],".",b.group$Well[1])
b.group$sample <- paste0("X",b.group$plate_id,".",b.group$Well)
head(b.group)
identical(colnames(a),b.group$sample)###确认a的行名和b的列明是一致的

# 筛选tumor cell
index <- which(b.group$TISSUE=="Tumor")
length(index)
group <- b.group[index,] #筛选的是行
head(group)
dim(group)###得到2343个肿瘤核心的细胞，5列

a.filt <- a[,index] #筛选的是列，筛选的是所有的肿瘤核心的细胞
dim(a.filt)
identical(colnames(a.filt),group$sample)

sessionInfo()

2. 通过Count 矩阵和meta 数据进行Seruat 对象构建，并进行数据的清洗和质量控制对数据的初步探索 主要基于以下一些指标

a.每个细胞中的检测的基因数目过滤细胞(每个细胞中检测到的基因不少于50个)
b基于每个基因在多少个细胞中检测到，过滤基因(每个基因至少在3个细胞里检测到)
c.每个细胞中的线粒体基因的占比，线粒体基因的占比越大，表明测序的时候，死细胞越多，通常线粒体基因的占比控制在小于5%或者10%以下，主要根据数据整体的质量进行评估
d.每个样本中的外源spike in的占比即ERCC 的比例(本次分析的过程中，为了和复现的原文保持一致的结果，ERCC 的比例控制在了小于40%

#############################################################
### 2、构建seurat对象，质控绘图----
# 2.1 构建seurat对象，质控
#In total, 2,343 cells from tumor cores were included in this analysis.
#quality control standards: 
#1) genes detected in < 3 cells were excluded; 筛选基因
#2) cells with < 50 total detected genes were excluded; 筛选细胞 
#3) cells with ≥ 5% of mitochondria-expressed genes were excluded. 筛选细胞
sessionInfo()
library("Seurat")

?CreateSeuratObject###可以通过查帮助文档进行创建seurat对象
sce.meta <- data.frame(Patient_ID=b.group$Patient_ID,
                       row.names = b.group$sample)
head(sce.meta)
table(sce.meta$Patient_ID)

# 这个函数 CreateSeuratObject 有多种多样的执行方式
scRNA = CreateSeuratObject(counts=a.filt,
                           meta.data = sce.meta,
                           min.cells = 3, 
                           min.features = 50)##细胞和基因的质控过滤，每个基因至少在3个细胞里表达，每个细胞至少检测得到50个features

#counts:a matrix-like object with unnormalized data with cells as columns and features as rows 
#meta.data:Additional cell-level metadata to add to the Seurat object
#min.cells: features detected in at least this many cells. 
#min.features:cells where at least this many features are detected.
head(scRNA@meta.data)
dim(scRNA@meta.data)
#nCount_RNA：the number of cell total counts
#nFeature_RNA：the number of cell's detected gene
summary(scRNA@meta.data)
scRNA@assays$RNA@counts[1:4,1:4]
# 可以看到，之前的counts矩阵存储格式发生了变化：4 x 4 sparse Matrix of class "dgCMatrix"

dim(scRNA)
#  20047  2342 仅过滤掉一个细胞
#########################################线粒体基因过滤
#接下来根据线粒体基因表达筛选低质量细胞
#Calculate the proportion of transcripts mapping to mitochondrial genes
table(grepl("^MT-",rownames(scRNA)))
#FALSE 

#20050 没有线粒体体基因
scRNA[[&#