概述
细胞聚类：基于基因表达信息，将表达谱相近的细胞聚为一类，表达差别大的细胞彼此分开。
Seurat使用
image.png
library(Seurat)
library(ggplot2)
library(dplyr)
options(stringsAsFactors=F)
##第1步：读入表达矩阵
data <- Read10X("../Data/filtered_gene_bc_matrices/")
##第2步：创建Seurat对象
train <- CreateSeuratObject(counts = data, project="train",min.cells = 3, min.features = 200)
train
expr_matrix <- train[["RNA"]]@counts
head(expr_matrix[,1:5])
write.table(expr_matrix,file="train.UMI.counts.xls",col.names=T,row.names=T,quote=F,sep="\t")
## 第3步：细胞质控
#1. 每个细胞检测的基因数目
#2. 每个细胞测序的UMI总数
#3. 每个细胞的线粒体基因比例
#其中，1.和2.在创建Seurat对象时候已经完成相关计算
train[["percent.mt"]] <- PercentageFeatureSet(train, pattern = "^MT-") ## 计算线粒体基因比例,pattern为匹配
##对于人：线粒体基因均以MT-开头(MT-ND1, MT-ND2, MT-CO1, MT-CO2, MT-ATP8, MT-ATP6, MT-CO3, MT-ND3, MT-ND4L, MT-ND4, MT-ND5, MT-ND6, MT-CYB)
##对于小鼠：线粒体基因均以mt-开头
##细胞质控信息存储在train@meta.data
head(train@meta.data) ###nCOUNTrna【每个细胞检测出的RNA的数量】 nfeture rna【每个细胞检测出基因表达的数量】
image.png
image.png
## 细胞质控信息可视化
pdf("train.cellqc.pdf")
VlnPlot(train, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()
ggsave("train.cellqc.png",VlnPlot(train, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3))
## 细胞质控指标相关性分析
plot1 <- FeatureScatter(train, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(train, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
pdf("train.cellqc.scatter.pdf")
plot1 + plot2
dev.off()
ggsave("train.cellqc.scatter.png",plot1 + plot2)
##细胞过滤，去掉不符合质控标准的细胞
train <- subset(train, subset = percent.mt < 10 & nFeature_RNA >= 250 & nFeature_RNA < 3000)
train ##取子集的意思
## 第4步：数据归一化
train <- NormalizeData(train, assay = "RNA", normalization.method = "LogNormalize", scale.factor = 10000)
##归一化后的数据存储在train[["RNA"]]@data里面
head(train[["RNA"]]@data[,1:5])
write.table(train[["RNA"]]@data,file="train.normdata.xls",quote=F,sep="\t",col.names=T,row.names=T)
## 第5步：高变基因鉴定和可视化【在细胞之间表达变量比较大的基因】
train <- FindVariableFeatures(train, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(train), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(train)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
pdf("train.hvg.pdf")
plot1 + plot2
dev.off()
ggsave("train.hvg.png",plot1 + plot2)
image.png
## 第6步：表达量标度化
train <- ScaleData(train,features=rownames(train))
##scaleData后的信息存储在train[["RNA"]]@scale.data里面
head(train[["RNA"]]@scale.data[,1:5])
image.png
## 第7步：PCA降维分析
train <- RunPCA(
object=train,
features=VariableFeatures(train),
npcs=50,
)
##检测前5个主成分的前5个特征基因（Positive和Negative各5个）
print(train[["pca"]], dims = 1:5, nfeatures = 5)
##可视化前2个PC的top30基因
pdf("train.pca.vizdim.pdf")
VizDimLoadings(train, dims = 1:2, nfeatures = 30, reduction = "pca")
dev.off()
ggsave("train.pca.vizdim.png",VizDimLoadings(train, dims = 1:2, nfeatures = 30, reduction = "pca"))
##基于前两个PC的细胞分布散点图
pdf("train.pca.dimplot.pdf")
DimPlot(train, reduction = "pca")
dev.off()
ggsave("train.pca.dimplot.png",DimPlot(train, reduction = "pca"))
image.png
##前15个PC的热图
pdf("train.pca.heatmap.pdf")
DimHeatmap(train, dims = 1:15, cells = 500, balanced = TRUE)
dev.off()
ggsave("train.pca.heatmap.png",DimHeatmap(train, dims = 1:15, cells = 500, balanced = TRUE))
## JackStraw图
#JackStraw :Determine statistical significance of PCA scores
##注意： 为了提高计算速度，可以改:num.replicate = 20
train <- JackStraw(train, reduction="pca",num.replicate = 20,prop.freq=0.01) ##reduction="pca"表示降维的方法是PCA，num.replicate = 20：表示抽样计算进行20次，freq=0.01每次抽样的比例为0.01
#ScoreJackStraw : Compute Jackstraw scores significance
train <- ScoreJackStraw(train, dims = 1:20)
#visualization for comparing the distribution of p-values for each PC with a uniform distribution (dashed line).
pdf("train.JackStrawPlot.pdf")
JackStrawPlot(train, dims = 1:15)
dev.off()
ggsave("train.JackStrawPlot.png",JackStrawPlot(train, dims = 1:15))
image.png
image.png
pdf("train.ElbowPlot.pdf")
ElbowPlot(train)
dev.off()
ggsave("train.ElbowPlot.png",ElbowPlot(train))
image.png
## 第8步：细胞聚类
train <- FindNeighbors(train, dims = 1:10)
train <- FindClusters(train, resolution = 0.5)
##提取各个细胞的聚类结果
cellcluster <- train@meta.data
cellcluster$cellid <- rownames(cellcluster)
cellcluster <- subset(cellcluster,select=c("cellid","seurat_clusters"))
write.table(cellcluster,file="train.cell.cluster.xls",quote=F,col.names=T,row.names=F,sep="\t")
## 采用TSNE对数据进行降维及可视化
train <- RunTSNE(object=train, dims=1:10)
pdf("train.tsne.pdf")
DimPlot(object = train,reduction="tsne")
dev.off()
ggsave("train.tsne.png",DimPlot(object = train,reduction="tsne"))
## 采用UMAP对数据进行降维及可视化
train <- RunUMAP(object = train, dims = 1:10)
pdf("train.umap.pdf")
DimPlot(object = train,reduction="umap")
dev.off()
ggsave("train.umap.png",DimPlot(object = train,reduction="umap"))
image.png
##第9步：标记基因鉴定和可视化
markers <- FindAllMarkers(train,logfc.threshold=0.5,test.use="wilcox",min.pct=0.25,only.pos=TRUE)
head(markers)
##排序，将同一个cluster的marker gene排在一起
#markers <- markers %>% group_by(cluster)
write.table(markers,file="train.cellmarker.xls",sep="\t",row.names=F,col.names=T,quote=F)
## 标记基因可视化
top2 <- markers %>% group_by(cluster) %>% top_n(n = 2, wt=avg_log2FC)
#1, 热图
pdf("train.marker.heatmap.pdf")
DoHeatmap(train, features = unique(top2$gene)) + NoLegend()
dev.off()
ggsave("train.marker.heatmap.png",DoHeatmap(train, features = unique(top2$gene)) + NoLegend())
top1 <- markers %>% group_by(cluster) %>% top_n(n = 1, wt=avg_log2FC)
#2，小提琴图
pdf("train.marker.vlnplot.pdf")
VlnPlot(train, features = unique(top1$gene))
dev.off()
ggsave("train.marker.vlnplot.png",VlnPlot(train, features = unique(top1$gene)))
#3，散点图
pdf("train.marker.featureplot.umap.pdf")
FeaturePlot(train, features = unique(top1$gene),reduction="umap")
dev.off()
ggsave("train.marker.featureplot.umap.png",FeaturePlot(train, features = unique(top1$gene),reduction="umap"))
#4，气泡图
pdf("train.marker.dotplot.pdf")
DotPlot(object = train, features = unique(top1$gene)) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
dev.off()
ggsave("train.marker.dotplot.png",DotPlot(object = train, features = unique(top1$gene)) + theme(axis.text.x = element_text(angle = 45, hjust = 1)))
##第10步：保存Seurat对象
saveRDS(train, file = "train.rds")
image.png
网友评论