003_10X基础下游分析流程

第一步是导入包

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(dplyr)
library(stringr)
library(magrittr)
setwd('~/Desktop/SingleCellDeepLearning/单细胞/10X/')

读取数据

这里注意一下，这个代码只读取了一套数据，如果是多样本处理的话要先进行样本整合

# 读取数据
# filtered_feature_bc_matrix 下面应该要有三个文件
# "barcodes.tsv.gz" # 这里注释了每一个细胞的 barcode
# "features.tsv.gz" 这里注释了每一个探针对应的基因名
# "matrix.mtx.gz"  行是探针 ID，列是细胞 ID
sce <- Read10X('test1/filtered_feature_bc_matrix/')   
# min.cells 某一个基因至少在多少个基因中表达
# min.features 某个细胞至少表达多少个基因
sce<- CreateSeuratObject(counts=sce, 
                          project = "10X_PBMC",
                          min.cells =3,
                          min.features =150)

计算线粒体 DNA 占比和重组 DNA 的比例，然后加到数据中

# 查看线粒体 dna 和 细胞重组 DNA
rownames(sce)[grepl('^mt-',rownames(sce),ignore.case = T)] 
rownames(sce)[grepl('^Rp[sl]',rownames(sce),ignore.case = T)]

sce[["percent.mt"]] = PercentageFeatureSet(sce, pattern = "^MT-")
head(sce[["percent.mt"]][,1])# 计算比例
rb.genes <- rownames(sce)[grep("^RP[SL]",rownames(sce))]
# 提取 count 数据
C<-GetAssaysce(object = sce, slot = "counts")
# 计算细胞重组 DNA在细胞中的比例
percent.rb <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
sce <- AddMetasce(sce, percent.rb, col.name = "percent.rb")

观察数据分布情况，主要是根据 4 个指标，rna 总数，特征个数，线粒体比例，重组 dna 比例

自行选择合理的阈值进行删选。
这一步也可以不做

plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.rb")
plot3 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2,plot3))
# 根据图 1 可以横坐标大于 50000 的 RNAcount 为异常
# 根据图 2 可以看出纵坐标大于 65000 左右的为异常
VlnPlot(sce, features = c(
  "percent.mt",
  "percent.rp",
  "nFeature_RNA",
  "nCount_RNA"), ncol = 4)
# 查看分布
# 利用分布合理的筛选数据
# 注意,数据筛选一定要合理
# 最高不可以超过 25%
sce[['percent.mt']] %>% summary()
sce[['percent.rb']] %>% summary()
sce[['nFeature_RNA']] %>% summary()
sce[['nCount_RNA']] %>% summary()

分析变化大的特征，将变化最大的特征提取并可视化

# 根据需求进行数据筛选
# 可以不用筛选的
sce = subset(sce,subset=(percent.mt < 25 &percent.rb < 55 &nFeature_RNA < 5000 & nCount_RNA<50000))
#数据标准化 
#归一化 去除样本/细胞效应， NormalizeData，默认方法 “LogNormalize”
sce = NormalizeData(sce, normalization.method =  "LogNormalize", 
              scale.factor = 10000)
GetAssay(sce,assay = "RNA")
#寻找高变基因，FindVariableFeatures，默认方法”vst”
#这些基因在一些细胞中高表达，
#在另一些细胞中低表达，这一步用FindVariableFeatures函数来执行。
#默认情况下，每个dataset返回2,000个基因。这些基因将用于下游分析，如PCA
sce<- FindVariableFeatures(sce,selection.method = "vst", nfeatures = 2000)
top20<-head(VariableFeatures(sce),20)# 提取差异最大的 top20 基因
plot1<- VariableFeaturePlot(sce)
plot2 <-LabelPoints(plot = plot1, points = top20, repel = TRUE)
plot2

PCA 降为

注意一下，这里有两种降为方式，第一种是根据差异最大的降为，还有一种是根据所有的特征降为

#去除基因效应，ScaleData
#用一个线性变换(“缩放scaling”)。这是在降维前（例如PCA）一个标准的预处理步骤，用
#ScaleData: （1）shift每个基因的表达，使细胞间的平均表达为0 （2）缩放每个基因的表达，使细胞间的差异为1。
#这一步给予下游分析同等的权重，这样那些非常高表达的基因就不会掩盖其他基因的变化，
# （根据所有的特征降维）
all.genes <- rownames(sce)
sce<-ScaleData(sce,all.genes)  #选择缩放所有基因，默认是top2000高变基因ScaleData(DATA)
sce<-RunPCA(sce,features = VariableFeatures(object =sce))
p1<-JackStraw(sce, num.replicate = 100)
p1<- ScoreJackStraw(p1, dims = 1:20)
JackStrawPlot(p1, dims = 1:20)  ##P值<0.05
ElbowPlot(sce)   ##明显拐点
DimPlot(sce, reduction = "pca")
#（根据差异的特征降为）
sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 
VizDimLoadings(sce, dims = 1:2, reduction = "pca")
DimHeatmap(sce, dims = 1:15, cells = 200, balanced = TRUE)
ElbowPlot(sce)# 查看各个主成分的方差
DimPlot(sce, reduction = "pca")
#确定有统计意义的主成分，需要几个主成分进行分析？一般20-30左右

TSNE

# 方差基本上在 15 以后不会变化
# 所以提取前 15 个
#细胞聚类
#首先我们在PCA空间里根据欧氏距离构建一个KNN图，并在任意两个细胞间要确定它们的边缘权重，这个过程#用FindNeighbors函数，这里的input就是前面我们定义的dataset的维数。
sce<- FindNeighbors(sce, dims = 1:20)

#再用FindClusters函数，该函数有一个“分辨率”的参数，该参数设置下游聚类的“粒度”，值
#越高，得到的聚类数越多。这个参数设置在0.4-1.2之间，
#对于3千个左右的单细胞数据通常会得到
#比较好的结果。对于较大的数据集，最佳分辨率通常会增加。
sce<- FindClusters(sce, resolution = 0.5) 
#Seurat提供了几种非线性的降维技术，如tSNE和UMAP。这些算法的目标是在低维空间中将相似的#细胞放在一起。上面所确定的基于图的集群中的单元应该在这些降维图上共同定位。作为
#UMAP和tSNE的input，我们建议使用与聚类分析的输入相同的PCs作为输入。
sce<-RunTSNE(sce,dims.use = 1:20)  ##tsne降维
sce<-RunUMAP(sce,dims = 1:20)  ##umap降维
##可视化
p1<-DimPlot(sce, reduction = "tsne")
p1
p2<-DimPlot(sce, reduction = "umap")
p2
CombinePlots(plots =list(p1, p2))
#以上便是利用Seurat进行降维分群的基本流程

set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
table(sce@meta.data$seurat_clusters,sce@meta.data$orig.ident)

phe=data.frame(cell=rownames(sce@meta.data),
               cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne')
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')

差异基因识别

#Seurat可以帮助你定义cluster的差异表达，默认情况下，它定义单个cluster的阳性和阴性的
# marker（与其他细胞群比较）。FindAllMarkers可以执行这个过程。
cluster1.markers <- FindMarkers(sce, ident.1 = 1, min.pct = 0.25)
#1簇与其他簇的差异，默认只返回上调的基因（即高表达），参数可改，only.pos = TRUE

#寻找cluster5与cluster0和3之间的差异marker
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)

#寻找每一个cluster里与其他所有细胞相比之后的差异marker:运行较久,test.use = "wilcox"
#检验方法默认wilcox，此外还有其他方法
pbmc.markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

##可视化marker基因，VlnPlot （展示在不同cluster里的表达可能分布），
#FeaturePlot（根据tSNE和PCA结果可视化基因表达）是最常用的可视化方法。你也可以用其他一
#些方法，例如：RidgePlot, CellScatter, DotPlot。
VlnPlot(sce, features = c("MS4A1", "CD79A"))
FeaturePlot(sce, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ"))
#还可以用DoHeatmap来画细胞和基因的热图，在这里我们画每一个cluster里的top10的marker。
top20 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
DoHeatmap(sce, features = top20$gene) + NoLegend()


new.cluster.ids <- c("Naive CD4 T", "monocyte", "CD8+ T", "unk", "B cell", "unk", "unk","unk","B cell","unk","unk","unk","unk")
names(new.cluster.ids) <- levels(sce) 
sce<- RenameIdents(sce, new.cluster.ids)
DimPlot(sce, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

如果您觉得我的文章对您有帮助，请点赞+关注，可以的话打个赏奖励一杯星巴克(～￣(OO)￣)ブ

Best Regards,  
Yuan.SH;
School of Basic Medical Sciences,  
Fujian Medical University,  
Fuzhou, Fujian, China.  
please contact with me via the following ways:  
(a) e-mail :yuansh3354@163.com