概述
細(xì)胞聚類：基于基因表達(dá)信息物咳，將表達(dá)譜相近的細(xì)胞聚為一類蹄皱，表達(dá)差別大的細(xì)胞彼此分開。
Seurat使用
image.png
library(Seurat)
library(ggplot2)
library(dplyr)
options(stringsAsFactors=F)
##第1步：讀入表達(dá)矩陣
data <- Read10X("../Data/filtered_gene_bc_matrices/")
##第2步：創(chuàng)建Seurat對(duì)象
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步：細(xì)胞質(zhì)控
#1. 每個(gè)細(xì)胞檢測(cè)的基因數(shù)目
#2. 每個(gè)細(xì)胞測(cè)序的UMI總數(shù)
#3. 每個(gè)細(xì)胞的線粒體基因比例
#其中压鉴，1.和2.在創(chuàng)建Seurat對(duì)象時(shí)候已經(jīng)完成相關(guān)計(jì)算
train[["percent.mt"]] <- PercentageFeatureSet(train, pattern = "^MT-") ## 計(jì)算線粒體基因比例,pattern為匹配
##對(duì)于人：線粒體基因均以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)
##對(duì)于小鼠：線粒體基因均以mt-開頭
##細(xì)胞質(zhì)控信息存儲(chǔ)在train@meta.data
head(train@meta.data) ###nCOUNTrna【每個(gè)細(xì)胞檢測(cè)出的RNA的數(shù)量】 nfeture rna【每個(gè)細(xì)胞檢測(cè)出基因表達(dá)的數(shù)量】
image.png
image.png
## 細(xì)胞質(zhì)控信息可視化
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))
## 細(xì)胞質(zhì)控指標(biāo)相關(guān)性分析
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)
##細(xì)胞過濾击蹲，去掉不符合質(zhì)控標(biāo)準(zhǔn)的細(xì)胞
train <- subset(train, subset = percent.mt < 10 & nFeature_RNA >= 250 & nFeature_RNA < 3000)
train ##取子集的意思
## 第4步：數(shù)據(jù)歸一化
train <- NormalizeData(train, assay = "RNA", normalization.method = "LogNormalize", scale.factor = 10000)
##歸一化后的數(shù)據(jù)存儲(chǔ)在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步：高變基因鑒定和可視化【在細(xì)胞之間表達(dá)變量比較大的基因】
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步：表達(dá)量標(biāo)度化
train <- ScaleData(train,features=rownames(train))
##scaleData后的信息存儲(chǔ)在train[["RNA"]]@scale.data里面
head(train[["RNA"]]@scale.data[,1:5])
image.png
## 第7步：PCA降維分析
train <- RunPCA(
object=train,
features=VariableFeatures(train),
npcs=50,
)
##檢測(cè)前5個(gè)主成分的前5個(gè)特征基因（Positive和Negative各5個(gè)）
print(train[["pca"]], dims = 1:5, nfeatures = 5)
##可視化前2個(gè)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"))
##基于前兩個(gè)PC的細(xì)胞分布散點(diǎn)圖
pdf("train.pca.dimplot.pdf")
DimPlot(train, reduction = "pca")
dev.off()
ggsave("train.pca.dimplot.png",DimPlot(train, reduction = "pca"))
image.png
##前15個(gè)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
##注意： 為了提高計(jì)算速度类咧，可以改:num.replicate = 20
train <- JackStraw(train, reduction="pca",num.replicate = 20,prop.freq=0.01) ##reduction="pca"表示降維的方法是PCA谴咸，num.replicate = 20：表示抽樣計(jì)算進(jìn)行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步：細(xì)胞聚類
train <- FindNeighbors(train, dims = 1:10)
train <- FindClusters(train, resolution = 0.5)
##提取各個(gè)細(xì)胞的聚類結(jié)果
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對(duì)數(shù)據(jù)進(jìn)行降維及可視化
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對(duì)數(shù)據(jù)進(jìn)行降維及可視化
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步：標(biāo)記基因鑒定和可視化
markers <- FindAllMarkers(train,logfc.threshold=0.5,test.use="wilcox",min.pct=0.25,only.pos=TRUE)
head(markers)
##排序，將同一個(gè)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)
## 標(biāo)記基因可視化
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螟炫，散點(diǎn)圖
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對(duì)象
saveRDS(train, file = "train.rds")
image.png
