单细胞测序生信分析

单细胞RNAseq常规分析流程（10X）

2021-05-24  本文已影响0人  QXPLUS

降维聚类，类型鉴定，差异分析

近年来，单细胞测序的成本逐渐降低，使得大众研究从组织水平转移到单个细胞水平的分析成为可能。本文主要是简单介绍一下SingleCellRNAseq（ScRNA）分析的基本流程，适用于大多数单细胞数据的前期分析，旨在介绍分析的流程和用到的方法，具体细节还需要感兴趣的你们深入挖掘。

主要分析流程

  1. 原始数据处理
  2. 细胞降维与聚类
  3. 细胞类群鉴定
  4. 基因差异表达分析
  5. 基因功能注释，细胞通讯，通路差异分析等

一，数据处理 -- Cellranger

mkdir cellranger_output
cd cellranger_output
cellranger count \
--id= sample1 \
--transcriptome= "refdata-cellranger-GRCh38-3.0.0/" \
--fastqs= "../rawdata/sample1/"  \
--r2-length = 98

二，细胞降维与聚类 -- 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为例）


precesseddata.png 
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)

## 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)

Normalize.png

2.2.2 找差异基因

project<- FindVariableFeatures(object = project, selection.method = "vst", nfeatures = 2000, mean.cutoff = c(0.1, 8),dispersion.cutoff=c(1, Inf))

FindVariableFeatures.png

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)

FindNeighbors.png

2.2.6 聚类

project <- FindClusters(project, resolution = 0.8)
resolution ：聚类分辨率，可以调整该参数以调节合适的聚类数目。

FindClusters.png

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)

对降维后的细胞可视化，即可看到降维后聚类的细胞分布情况。
尽量选择不同类别边界清晰的一组参数进行后续的细胞类群鉴定。


metadata.png

DimPlot(object=project,dims = c(1, 2), reduction ="umap", group.by = c("seurat_clusters"),label = TRUE)

UMAP_by_Sample.png
tSNE_by_Sample.png

保存聚类后的seurat对象为project.rds 用于后续细胞类群鉴定。

三，细胞类群鉴定

在这里，以UMAP图为基础，进行细胞类群鉴定

3.1 根据已有marker 鉴定细胞类型

首先需要自行整理出markergene.txt文件

markergene.txt.png

思路：将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")
UMAP.png Bmarker.png

先初步判定cluster 4,5,11为B cell

Epithelial cell marker.png

这个就很明显了，cluster 3,13,14,15,17, 为 上皮细胞，cluster 9,先打个问号吧。

myeloid cellmarker.png

cluster 6,7 认为是 Myeloid cell


Mast cell marker.png

cluster21 为 mast cell

T cell marker.png

cluster 1,2,16,18,20 为T cell, cluster 9也有可能是T cell?

Endothelial cell marker.png

cluster 12 为 , Endothelial cell

好了，总结一下，到目前为止鉴定出来的cluster分别是

3.2 根据差异基因鉴定细胞类型

除了去map已知文献报道的marker gene, 我们还可以用生信分析的基本功--计算差异基因，再根据差异基因去网站（eg. CellMarker:http://bio-bigdata.hrbmu.edu.cn/CellMarker/）上查询所属的细胞类型，综合做出判断。

cluster8-deg.png
在CellMarker上，查询COL1A1。
CellMarker-COL1A1.png
发现很有可能是成纤维细胞，于是再次查询文献报道过的Fibroblast marker gene, 最后综合判断cluster 8,10 为Fibroblast cells.
Fibroblast marker.png

我们在cluster9中并没有看到明显的独立的cluster特有的差异基因，综合上面的结果，判定cluster 为T cell


Cluster9-deg.png

cluster19 太小了，好像刚刚漏掉了，因为该样本是肠癌样本，cluster19map 上了肠神经胶质细胞的marker , 同时，cluster19特异的gene 查询到属于Schwann cell，属于神经胶质细胞的一种，所以判定cluster19为 Enteric glial cells。


Enteric glial cells.png
Cluster19-deg.png sina.blog.png

好啦，耗时大半天，终于确定细胞类型啦，我们来再次总结一下：

UMAPplot_by_celltype.png

由于我们下载到的GEO数据集是已经给出细胞类型鉴定结果的，所以我们来看一下作者的结果，还是有一点点差异的，不过大的细胞类群鉴定结果几乎相同。

作者的细胞类型如下：


cell_type.png
上一篇下一篇

猜你喜欢

热点阅读