单细胞测序流程(五)t-sne聚类分析和寻找marker基因

系列文章目录

单细胞测序流程(一)简介与数据下载

单细胞测序流程(二)数据整理

单细胞测序流程(三)质控和数据过滤——Seurat包分析,小提琴图和基因离差散点图

单细胞测序流程(四)主成分分析——PCA


 本期主讲内容——t-sne聚类分析和寻找marker基因

介绍:T-SNE是一种用于探索高维数据的非线性降维算法。它将多维数据映射到适合于人类观察的两个或多个维度。

单细胞测序的意义:​根本在于细胞的异质性。就是说,细胞与细胞之间存在个体差异性,即便是出于同一位置的细胞,也可能在基因表达等方面存在一些差异。


一、课前准备

之前所使用的数据

R语言的IDE

二、过程

聚类分析的目的是根据细胞的差别所进行细胞分群,两个亚群的互相靠近代表着两个亚群之间的相关。

寻找marker基因目的是寻找在每个亚群中的标志基因,即提到那个基因立马想到对应的亚群。

直接运行代码就可以出线结果

代码如下

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("singscore")
#BiocManager::install("GSVA")
#BiocManager::install("GSEABase")
#BiocManager::install("limma")

#install.packages("devtools")
#library(devtools)
#devtools::install_github('dviraran/SingleR')
library(limma)
library(Seurat)
library(dplyr)
library(magrittr)
setwd(目录")             #设置工作目录

#读取文件,并对重复基因取均值
rt=read.table("geneMatrix.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#将矩阵转换为Seurat对象,并对数据进行过滤
pbmc <- CreateSeuratObject(counts = data,project = "seurat", min.cells = 3, min.features = 50, names.delim = "_",)
#使用PercentageFeatureSet函数计算线粒体基因的百分比
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pdf(file="04.featureViolin.pdf",width=10,height=6)           #保存基因特征小提琴图
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 50 & percent.mt < 5)    #对数据进行过滤

#测序深度的相关性绘图
pdf(file="04.featureCor.pdf",width=10,height=6)              #保存基因特征相关性图
plot1 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt",pt.size=1.5)
plot2 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",,pt.size=1.5)
CombinePlots(plots = list(plot1, plot2))
dev.off()

#对数据进行标准化
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#提取那些在细胞间变异系数较大的基因
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 1500)
#输出特征方差图
top10 <- head(x = VariableFeatures(object = pbmc), 10)
pdf(file="04.featureVar.pdf",width=10,height=6)              #保存基因特征方差图
plot1 <- VariableFeaturePlot(object = pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
dev.off()

##PCA分析
pbmc=ScaleData(pbmc)                     #PCA降维之前的标准预处理步骤
pbmc=RunPCA(object= pbmc,npcs = 20,pc.genes=VariableFeatures(object = pbmc))     #PCA分析

#绘制每个PCA成分的相关基因
pdf(file="05.pcaGene.pdf",width=10,height=8)
VizDimLoadings(object = pbmc, dims = 1:4, reduction = "pca",nfeatures = 20)
dev.off()

#主成分分析图形
pdf(file="05.PCA.pdf",width=6.5,height=6)
DimPlot(object = pbmc, reduction = "pca")
dev.off()

#主成分分析热图
pdf(file="05.pcaHeatmap.pdf",width=10,height=8)
DimHeatmap(object = pbmc, dims = 1:4, cells = 500, balanced = TRUE,nfeatures = 30,ncol=2)#1:4是热图分成几个,ncol是分为几列。
dev.off()

#每个PC的p值分布和均匀分布
pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
pdf(file="05.pcaJackStraw.pdf",width=8,height=6)
JackStrawPlot(object = pbmc, dims = 1:20)
dev.off()
#
#
#
#
#这一章的代码从这里开始
##TSNE聚类分析
pcSelect=20
pbmc <- FindNeighbors(object = pbmc, dims = 1:pcSelect)                #计算邻接距离
pbmc <- FindClusters(object = pbmc, resolution = 0.5)                  #对细胞分组,优化标准模块化
pbmc <- RunTSNE(object = pbmc, dims = 1:pcSelect)                      #TSNE聚类
pdf(file="06.TSNE.pdf",width=6.5,height=6)
TSNEPlot(object = pbmc, do.label = TRUE, pt.size = 2, label = TRUE)    #TSNE可视化
dev.off()
write.table(pbmc$seurat_clusters,file="06.tsneCluster.txt",quote=F,sep="\t",col.names=F)

##寻找差异表达的特征
logFCfilter=0.5
adjPvalFilter=0.05
pbmc.markers <- FindAllMarkers(object = pbmc,
                               only.pos = FALSE,
                               min.pct = 0.25,
                               logfc.threshold = logFCfilter)
sig.markers=pbmc.markers[(abs(as.numeric(as.vector(pbmc.markers$avg_logFC)))>logFCfilter & as.numeric(as.vector(pbmc.markers$p_val_adj))<adjPvalFilter),]
write.table(sig.markers,file="06.markers.xls",sep="\t",row.names=F,quote=F)

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
#绘制marker在各个cluster的热图
pdf(file="06.tsneHeatmap.pdf",width=12,height=9)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
dev.off()

#绘制marker的小提琴图
pdf(file="06.markerViolin.pdf",width=10,height=6)
VlnPlot(object = pbmc, features = c("IGLL5", "MBOAT1"))
dev.off()

#绘制marker在各个cluster的散点图
pdf(file="06.markerScatter.pdf",width=10,height=6)
FeaturePlot(object = pbmc, features = c("IGLL5", "MBOAT1"),cols = c("green", "red"))
dev.off()

#绘制marker在各个cluster的气泡图
pdf(file="06.markerBubble.pdf",width=12,height=6)
cluster10Marker=c("MBOAT1", "NFIB", "TRPS1", "SOX4", "CNN3", "PIM2", "MZB1", "MS4A1", "ELK2AP", "IGLL5")
DotPlot(object = pbmc, features = cluster10Marker)
dev.off()


这张图可以看出这个基因在各个亚群之间的表达情况,颜色越深就代表表达越高,可以观察这个基因是不是差异基因。  

这个图是marker气泡图,它代表什么? 气泡的大小代表着基因在亚群中的占比颜色的深浅代表着她的表达程度,颜色越深表达程度越高。

 这个图就是利用T-SNE算法对所有的细胞进行分群,相互靠近的亚群代表着二者之间有联系。

这个图是基因在各个亚群中的小提琴图,一个黑点代表一个细胞,纵轴代表着表达水平,越高就代表着表达程度越高,

细胞亚群top20基因聚类热图,绘制这个图是为了进一步确认较为理想的分群结果,并确认比对找到可能作为亚群的特征基因。对每个亚群的top20的基因进行聚类分析观察结果。聚类分析不仅可以获得不同细胞亚群 top20 基因的表达模式,还可以判断同一亚群的基因是否能够聚集成类。因为这些同类的基因可能具有相似的功能,来自同一类型的细胞。


注意事项:因为这次的结果很多取决于之前的数据,所以必须要把上一届的内容也要用到,我已放进代码块中,如果之前内容已经取得,不想重复获取,只需将之前代码中的PDF开头的代码行到dev.off()代码行全部删掉,就不会将上一张的图形再绘制出来。

单细胞测序流程(五)t-sne聚类分析和寻找marker基因到这就结束了
下一章会讲解单细胞的细胞类型的注释
我所做的所有分析与教程的代码都会在我的个人公众号中,请打开微信搜索“生信学徒”进行关注,欢迎生信的研究人员和同学前来讨论分析。
ps:公众号刚刚建立比较简陋,但是该有的内容都不会少。

 

  • 10
    点赞
  • 60
    收藏
  • 打赏
    打赏
  • 12
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:数字20 设计师:CSDN官方博客 返回首页
评论 12

打赏作者

生信学徒

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值