Seurat使用
2020-11-18 本文已影响0人
谢京合
这里我就用示例数据讲解啦。可以在GSE数据库下载你所关注的病症的数据集。
下载链接“https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE106273”
需要保证在这个目录下E:/10×单细胞矩阵练习/GSE106273_RAW_TEST/，
包含有三个文件（barcodes.tsv/genes.tsv/matrix.mtx）,并且命名要严格按照括弧当中所示的进行。
先把这些包都装上。安装方法自己看着办，一般是常规的方法，如果常规的方法不行就自己去官网查阅。
library(dplyr)
library(Seurat)
library(patchwork)
library(Matrix)
library(export)
然后按照seurat官网出示的教程来开始。
里面具体操作的含义有些其实不太懂，但是没关系呀，得到自己想要的结果就可以了。
### 构建Seurat对象
pbmc.data <- Read10X(data.dir = "F:/Seurat/")###这里会报错，这里的路径名字真的不能是中文！！！绝了。
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc", min.cells = 3, min.features = 200)
pbmc
### 预处理与质量控制
#计算线粒体RNA比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data, 5)
#使用violin plot查看QC指标分布，如图1所示
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
graph2ppt(file="1.violinplot_QC.pptx", width=12, aspectr=1.5)
#查看feature-feature之间关系，如图2所示
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
graph2ppt(file="2.feature_feature.pptx", width=12, aspectr=1.5)
#去除表达基因数量少于500或多于6000的细胞，去除线粒体RNA比例高于40%的细胞.这个具体标准的筛选，根据之前violin plot来确定。
pbmc <- subset(pbmc, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 40)
dim(pbmc)
上述步骤完成之后，就会得到你后期用于筛选的数据集，你可以把它保存一下，方便下次直接导入。
saveRDS(pbmc, file = "F:/Seurat/subset.rds")
接着就是筛选基因、降维分析和分群啦。
### 校正表达值
#校正细胞的表达值，将count数量根据表达总量矫正之后，乘以10000，并进行对数转换
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#以上代码的计算原理如下
#pbmc@assays$RNA@data <-log(t(t(pbmc@assays$RNA@counts)*10000/colSums(pbmc@assays$RNA@counts) + 1))
#以上是默认参数.
#下列代码实现相同功能
#pbmc <- NormalizeData(pbmc)
### 鉴定高变异基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#查看基因的平均表达与标准方差之间的关系，如图3所示
plot1 <- VariableFeaturePlot(pbmc)
top10 <- head(VariableFeatures(pbmc), 10)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
graph2ppt(file="3.top2000genes.pptx", width=10, aspectr=1.5)
### 归一化表达值
#对基因在细胞间的表达量进行线性变换，使平均值为0，方差为1，去除高表达基因的影响
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
#以上代码的计算原理如下:
#pbmc@assays$RNA@scale.data <- t(scale(t(pbmc@assays$RNA@data)))
#scale是个十分耗时的过程，为缩短计算时间，ScaleData默认仅对高变异基因进行scale；不论是对全部基因或仅对高变异基因进行scale，并不影响下游的PCA和聚类过程
#ScaleData同样可以用来去除其他因素（如细胞周期状态或线粒体RNA比例）带来的影响
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
dim(pbmc)
### 线性降维
#在scale.data上运行PCA，默认使用的基因为高变异基因
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#观察PC上重要基因
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
#在PC1与PC2上表示基因，如图4所示，这里展示的实际上是细胞进行PCA降维时每个基因的系数，即PCA feature loading，根据PCA算法，有pbmc@reductions$pca@cell.embeddings ＝ t(pbmc@assays$RNA@scale.data[VariableFeatures(pbmc),]) %*% pbmc@reductions$pca@feature.loadings
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
graph2ppt(file="4.PC1nPC2_value.pptx", width=9, aspectr=1.5)
#细胞的PCA降维，如图5所示
DimPlot(pbmc, reduction = "pca")
graph2ppt(file="5.PCA降维.pptx", width=9, aspectr=1.5)
#查看每个细胞中单个基因在不同PC维度上的得分，如图6所示
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
graph2ppt(file="6.单个PC维度_heatmap.pptx", width=9, aspectr=1.5)
DimHeatmap(pbmc, dims = 1:9, cells = 500, balanced = TRUE)
graph2ppt(file="7-1.1到9个PC维度_heatmap.pptx", width=12, aspectr=1.5)
DimHeatmap(pbmc, dims = 10:18, cells = 500, balanced = TRUE)
graph2ppt(file="7-2.9到18个PC维度_heatmap.pptx", width=12, aspectr=1.5)
DimHeatmap(pbmc, dims = 19:27, cells = 500, balanced = TRUE)
graph2ppt(file="7-3.18到27个PC维度_heatmap.pptx", width=12, aspectr=1.5)
###这里可以保存数据,这里最开始写错了，不是这里的数据用于doublefinder，而是后面的，不过既然写了，先放着吧。
saveRDS(pbmc, file = "F:/Seurat/PCA.rds")
确定PCA维度数量
### 确定PCA维度数量
#根据JackStrawPlot来确定下一步降维作用的PC数量，如图７所示。P值越小说明PC越重要，应该纳入下一步分析
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)#这里dims = 1:20表示看前20个PC的P值，因为第一次算出来的结果真的惊人。而且这个dim最大只能是20.
JackStrawPlot(pbmc, dims = 1:20) #用JackStrawPlot绘图，绘制前10个PC的P值分布曲线
graph2ppt(file="8.PCA维度_JackStraw.pptx", width=12, aspectr=1.5)
#由图7可以看出，在PC10-PC12之后，PC的P值出现显著降低
#或根据ElbowPlot确定PC数量，如图8所示。标准差较大的PC所解释的数据变异性更大，应当纳入下一步分析
ElbowPlot(pbmc)
graph2ppt(file="9.PCA维度_ElbowPlot.pptx", width=12, aspectr=1.5)
#由图8可看出，在PC9-PC10后出现标准差数据拐点。综上，下游分析中选择了前10个PC
接着，是细胞聚类，细胞聚类这里很魔性的，人为控制的因素占主导。换成人话就是：具体分几个细胞群，由你来控制。
但是也不能控制的太随意，还是需要些指标来参考滴。这里我就提供clustree这个包吧。
### 细胞聚类
#这里的这个dim表示你需要聚类的类别数量，这里官方提供是10，但是具体操作过程中考虑实际情况选大一些。这里我选20.
pbmc <- FindNeighbors(pbmc, dims = 1:20)
#这里的resolution参数可以修改细胞类型的多少。resolution越大细胞的类型越多，所以可以尝试。0.01-1群细胞，0.05-2群细胞，0.1-4群细胞，0.5-8群细胞。
pbmc <- FindClusters(pbmc, resolution = 0.05)
##或者上述的resolution还可以用另外一种方式将其划分，便于找到最优的分群数目。
pbmc <- FindClusters(pbmc, resolution = c(0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1))
##查看上述不同的resolution的时候，聚类的情况。
pbmc@meta.data %>% View()
#install.packages('clustree')
library(clustree)
clustree(pbmc@meta.data, prefix = "RNA_snn_res.")
graph2ppt(file="10.clustree.pptx", width=9, aspectr=1.5)
#根据上述的结果，找到最佳的分群resolution的值为0.3。所以这里需要再跑一边即可。
pbmc <- FindClusters(pbmc, resolution = 0.3)
然后就是降维之后的绘图呀。
### 非线性降维（UMAP/tSNE）
#UMAP降维
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap")
graph2ppt(file="11.umap.pptx", width=9, aspectr=1.5)
#tSNE降维
pbmc <- RunTSNE(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "tsne")
graph2ppt(file="12.tsne.pptx", width=9, aspectr=1.5)
###这里保存数据，用于鉴定singlelets和doublelets
saveRDS(pbmc, file = "F:/Seurat/ForDoublets.rds")
差异表达分析，里面具体的内容自己慢慢去看哈，都写得很清楚。
### 差异表达分析
# find all markers of cluster 1# 找到cluster1当中的所有marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
cluster3.markers <- FindMarkers(pbmc, ident.1 = 3, min.pct = 0.25)
head(cluster3.markers, n = 5)
#使用不同的假设检验方法（test.use = "roc"）来进行差异表达分析
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
# find all markers distinguishing cluster 5 from clusters 0 and 3#找到cluster5当中和clusters0和clusters3能够区分的marker。
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# 找到每一个cluster当中的marker，并且只展示阳性的marker。
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
##上述步骤的原理如下：
#pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
head(pbmc.markers, n = 10)
##将这个数据保存一下，用于亚群鉴定得时候用。
write.table(pbmc.markers,file = 'pbmc.markers.txt',sep = '\t')
#差异基因可视化，此外还可以通过RidgePlot, CellScatter, DotPlot等进行展示，这里可以每个亚群筛选一个，也可以根据需要。
##选择每个cluster的前三个基因。此外，这里的这个名字识别的是最后一列gene的名字。
VlnPlot(pbmc, features = c("AGR2", "TFF2", "HMOX1", "HMGN1", "FABP1", "BEX2"))
graph2ppt(file="13.clusterMarker_VlnPlot.pptx", width=20, aspectr=1.5)
# 使用原始count绘制
VlnPlot(pbmc, features = c("AGR2", "TFF2", "HMOX1", "HMGN1", "FABP1", "BEX2"), slot = "counts", log = TRUE)
graph2ppt(file="14.clusterMarker_VlnPlot_count.pptx", width=20, aspectr=1.5)
#差异基因在细胞分区图上的展示。
FeaturePlot(pbmc, features = c("AGR2", "TFF2", "HMOX1", "HMGN1", "FABP1", "BEX2"))
graph2ppt(file="15.clusterMarker_FeaturePlot.pptx", width=15, aspectr=1.5)
#差异基因可视化RidgePlot
RidgePlot(pbmc, features = c("AGR2", "TFF2", "HMOX1", "HMGN1", "FABP1", "BEX2"))
graph2ppt(file="16.clusterMarker_RidgePlot.pptx", width=18, aspectr=1.5)
#差异基因可视化CellScatter,参数还没搞懂。
#CellScatter(pbmc, features = c("AGR2", "UCHL1", "HMOX1", "MAGEA4", "FABP1", "TFF2", "CXCL1", "BEX2","SOX2"))
#差异基因可视化DotPlot
DotPlot(pbmc, features = c("AGR2", "TFF2", "HMOX1", "HMGN1", "FABP1", "BEX2"))
graph2ppt(file="17.clusterMarker_DotPlot.pptx", width=12, aspectr=1.5)
#热图显示。每个cluster按照avg_logFC筛选排名前10的，进行热图展示。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
graph2ppt(file="18.markersHeatmap.pptx", width=9, aspectr=1.5)
saveRDS(pbmc, file = "F:/Seurat/ClustersPlot_beyond.rds")
接着就是细胞分类之后的鉴定啦。
### 鉴定细胞类型
#根据已知细胞类型的marker gene对各个cluster的细胞进行命名，如图14所示
new.cluster.ids <- c("1", "2", "3")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
graph2ppt(file="鉴定细胞类型.pptx", width=9, aspectr=1.5)
### 保存结果
saveRDS(pbmc, file = "../output/final.rds")
