logo
搜索
我的图书馆

发文章

发文工具

其他工具

留言交流
搜索
分享

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

 向向n38atyvda4 2022-10-09 发布于广东

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

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

主要分析流程

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

一，数据处理 -- 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为例）


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)
  • 根据表达矩阵创建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)

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)

  • umap.method = "umap-learn" 基于python 的umap 包

  • 执行之前，先执行 -- py_module_available(module = 'umap')


    umap.png

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


metadata.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)

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分别是

  • 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./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

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

  • 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
UMAPplot_by_celltype.png

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

作者的细胞类型如下：


cell_type.png
  • 个人认为，细胞类群鉴定本身不是一件特别精确的事情，他需要同时综合样本来源，常见细胞类群的关键基因的生物学知识，以及基于生信分析手段得到的差异基因的表达情况等来判定，总之能自圆其说就好，不必针对个别cluster锱铢必较啦。当然，如果你们有更好的方法，欢迎留言哦。

    本站是提供个人知识管理的网络存储空间，所有内容均由用户发布，不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息，谨防诈骗。如发现有害或侵权内容，请点击一键举报。
    转藏 分享 献花（0）

    来自： 向向n38atyvda4 > 《单细胞测序》

    举报/认领

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多