降维聚类，类型鉴定，差异分析
近年来，单细胞测序的成本逐渐降低，使得大众研究从组织水平转移到单个细胞水平的分析成为可能。本文主要是简单介绍一下SingleCellRNAseq（ScRNA）分析的基本流程，适用于大多数单细胞数据的前期分析，旨在介绍分析的流程和用到的方法，具体细节还需要感兴趣的你们深入挖掘。
主要分析流程
- 原始数据处理
- 细胞降维与聚类
- 细胞类群鉴定
- 基因差异表达分析
- 基因功能注释，细胞通讯，通路差异分析等
一，数据处理 -- Cellranger
-
在本篇文章中，我们都是基于10X来源的单细胞数据
Cellranger输入的每个样本包含3个文件：L1, R1,R2
rawdata.sample1.png
mkdir cellranger_output
cd cellranger_output
cellranger count \
--id= sample1 \
--transcriptome= "refdata-cellranger-GRCh38-3.0.0/" \
--fastqs= "../rawdata/sample1/" \
--r2-length = 98
-
tips: --id 是cellranger 输出文件，切记不可以与rawdata在同一文件夹下同名，否则会出现环境错误等奇怪的报错
Cellranger输出的每个样本也包含3个文件：barcode, features(genes),matrix
cellranger.out.png
很多时候，我们都可以从GEO网站上直接下载到上述3个文件的。
二，细胞降维与聚类 -- Seurat
2.1 创建Seurat 对象
2.1.1 基于10X标准输出文件（barcode, features,matrix）创建Seurat 对象
读取10X标准输出文件
# 将barcode, features,matrix三个文件，放入 sample1文件夹中
data <- Read10X(data.dir =samples1,gene.column = 1)
feature<-read.table("features.tsv.gz",head=F,sep="\t")
gene<-as.character(feature[match(rownames(data),feature$V1),"V2"])
创建Seurat对象
project<-CreateSeuratObject(counts =data, project = samples1)
2.1.2 基于GEO处理后的表达矩阵创建Seurat 对象
如果你下载到的数据是gene-sampleID的表达矩阵，那么这个数据集是别人已经处理好了的。你可以直接使用，或者自己重新降维，聚类，细胞类群鉴定等等。（以 GSE132465为例）
datadir = "GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix.txt"
metadata = "GSE132465_GEO_processed_CRC_10X_cell_annotation.txt"
outdir = "./Dimention"
nfeature<-100
ncount<-200
mito<-10
ribo<-100
red<-10
## read gene-sampleID matrix and metadata ---------------
data<-fread(datadir,sep = "\t", header = T,check.names = F)
data<-data.frame(data,check.names = F, row.names = 1)
metadata<-fread(metadata,sep = "\t", header = T,check.names = F)
metadata<-data.frame(metadata,check.names = F, row.names = 1)
- 根据表达矩阵创建SeuratObject
（在生信分析中，这种对象有很多，比如你们可能会用到的CellChatObject, 最好是搞清楚他的组成和使用）
## Create SeuratObject ------------------------------------
project = CreateSeuratObject(counts = data,metadata = metadata)
project = AddMetaData(object = project, metadata = metadata)
2.1.3 数据清洗
table(grepl("^RP[SL][[:digit:]]",rownames(project)))
project[["mito.per"]] = PercentageFeatureSet(project, pattern = "^MT-")
project[["ribo.per"]] = PercentageFeatureSet(project, pattern = "^RP[SL][[:digit:]]")
project[["redcell.per"]] = PercentageFeatureSet(project, pattern = "^HB-")
project <- subset(project,
subset = nFeature_RNA >= nfeature &
nCount_RNA >= ncount &
mito.per <= mito &
ribo.per <= ribo &
redcell.per<=red )
数据保存，这个很重要哦，方便以后对结果的补充和复现呀。
2.2 数据降维聚类的流程
主要流程为：log：Normalize --> FindVariableFeatures --> ScaleData --> RunPCA --> FindNeighbors --> FindClusters --> RunTSNE/RunUMAP
2.2.1 log : NormalizeData
project<- NormalizeData(object = project, scale.factor = 10000)
2.2.2 找差异基因
project<- FindVariableFeatures(object = project, selection.method = "vst", nfeatures = 2000, mean.cutoff = c(0.1, 8),dispersion.cutoff=c(1, Inf))
2.2.3 标准化
project <- ScaleData(project,verbose = FALSE,features=rownames(project))
2.2.4 PCA
project <- RunPCA(project, npcs = 10, verbose = FALSE)
2.2.5 构建图
project <- FindNeighbors(project, reduction = "pca", dims = 1:10)
2.2.6 聚类
project <- FindClusters(project, resolution = 0.8)
resolution ：聚类分辨率，可以调整该参数以调节合适的聚类数目。
2.2.7 降维 tsne /umap
tsne
project <- RunTSNE(object = project, reduction = "pca",dims = 1:10)
umap
project <- RunUMAP(object = project, reduction = "pca",dims = 1:10,umap.method = "umap-learn", check_duplicates = FALSE)
- umap.method = "umap-learn" 基于python 的umap 包
-
执行之前，先执行 -- py_module_available(module = 'umap')
umap.png
对降维后的细胞可视化，即可看到降维后聚类的细胞分布情况。
尽量选择不同类别边界清晰的一组参数进行后续的细胞类群鉴定。
- 可以通过调节 resolution,改变cluster (本次聚类中，resolution=0.4 RNA_snn_res.0.4， 聚成21 类)
DimPlot(object=project,dims = c(1, 2), reduction ="umap", group.by = c("seurat_clusters"),label = TRUE)
保存聚类后的seurat对象为project.rds 用于后续细胞类群鉴定。
三，细胞类群鉴定
在这里，以UMAP图为基础，进行细胞类群鉴定
3.1 根据已有marker 鉴定细胞类型
首先需要自行整理出markergene.txt文件
思路：将marker gene map到UMAP图上，综合确定细胞类群。
获取数据中存在的marker gene
gene<-read.table("markergene.txt",sep="\t",header=T,check.names=F)
data<-project@assays$RNA@data
target_gene<-intersect(gene[,2],rownames(data))
将marker gene map到UMAP图上
DefaultAssay(project) <- "RNA"
subdata <- lapply(1:length(target_gene),function(tg){
sub = data.frame(barcode=colnames(data),
exp=data[target_gene[tg],],
Gene=target_gene[tg],
celltype=pdf_name)
umap = as.data.frame(project@reductions$umap@cell.embeddings)
sub <- cbind(sub,umap)
return(sub)
})
subdata<-do.call('rbind',subdata)
ggplot(subdata,
aes(x=UMAP_1, y=UMAP_2, color=exp)) +
geom_point(size=0.5,alpha=1)+
theme_bw()+facet_wrap(~Gene,ncol=col_num)+theme(legend.position="right",panel.grid=element_blank(),axis.text=element_blank(),axis.title=element_text(size=rel(2),color="black"),strip.background=element_rect(fill="white",colour="white"),axis.line =element_line(colour="black",size=),axis.ticks.length = unit(0,"cm"),strip.text = element_text(face="bold", size=rel(2)),plot.title=element_text(size=rel(3),color="black",hjust=0.5))+
scale_color_gradient(low="#cee7f4", high="red")
先初步判定cluster 4,5,11为B cell
这个就很明显了，cluster 3,13,14,15,17, 为 上皮细胞，cluster 9,先打个问号吧。
cluster 6,7 认为是 Myeloid cell
cluster21 为 mast cell
cluster 1,2,16,18,20 为T cell, cluster 9也有可能是T cell?
cluster 12 为 , Endothelial cell
好了，总结一下，到目前为止鉴定出来的cluster分别是
- B cells: cluster 4,5,11
- T cells: cluster 1,2,16,18,20
- epithelial cells: cluster 3,13,14,15,17
- endothelial cells:cluster 12
- Myeloid cells: cluster 6,7
- Mast cells: cluster 21
一共21个cluster, 现在不确定的就剩下cluster 9，8，10啦。
3.2 根据差异基因鉴定细胞类型
除了去map已知文献报道的marker gene, 我们还可以用生信分析的基本功--计算差异基因，再根据差异基因去网站（eg. CellMarker:http://bio-bigdata.hrbmu.edu.cn/CellMarker/）上查询所属的细胞类型，综合做出判断。
在CellMarker上，查询COL1A1。
发现很有可能是成纤维细胞，于是再次查询文献报道过的Fibroblast marker gene, 最后综合判断cluster 8,10 为Fibroblast cells.
我们在cluster9中并没有看到明显的独立的cluster特有的差异基因，综合上面的结果，判定cluster 为T cell
cluster19 太小了，好像刚刚漏掉了，因为该样本是肠癌样本，cluster19map 上了肠神经胶质细胞的marker , 同时，cluster19特异的gene 查询到属于Schwann cell，属于神经胶质细胞的一种，所以判定cluster19为 Enteric glial cells。
好啦，耗时大半天，终于确定细胞类型啦，我们来再次总结一下：
- B cells: cluster 4,5,11
- T cells: cluster 1,2,9,16,18,20
- Epithelial cells: cluster 3,13,14,15,17
- Endothelial cells:cluster 12
- Myeloids: cluster 6,7
- Mast cells: cluster 21
- Fibroblast: cluster 8,10
- Enteric glial cells: 19
由于我们下载到的GEO数据集是已经给出细胞类型鉴定结果的，所以我们来看一下作者的结果，还是有一点点差异的，不过大的细胞类群鉴定结果几乎相同。
作者的细胞类型如下：
- 个人认为，细胞类群鉴定本身不是一件特别精确的事情，他需要同时综合样本来源，常见细胞类群的关键基因的生物学知识，以及基于生信分析手段得到的差异基因的表达情况等来判定，总之能自圆其说就好，不必针对个别cluster锱铢必较啦。当然，如果你们有更好的方法，欢迎留言哦。
