前言
本次主要使用Seurat比较STAR与Cellranger的输出结果，只会进行简单的聚类工作。
数据读入
> library(Seurat)
> library(dplyr)
> library(magrittr)
> library(gtools)
> library(stringr)
> library(Matrix)
> setwd("D://data/ScRNAcode")
##读入STAR数据
> matrix.dir="STAR/"
> barcode.path <- paste0(matrix.dir,"barcodes.tsv")
> features.path <- paste0(matrix.dir,"features.tsv")
> matrix.path <- paste0(matrix.dir, "matrix.mtx")
> STARmatrix <- readMM(file = matrix.path)
> feature.names = read.delim(features.path,
+ header = FALSE,
+ stringsAsFactors = FALSE)
> barcode.names = read.delim(barcode.path,
+ header = FALSE,
+ stringsAsFactors = FALSE)
> colnames(STARmatrix) = barcode.names$V1
> rownames(STARmatrix) = feature.names$V2
> STARmatrix<-as.matrix(STARmatrix)
> STARmatrix[1:6,1:6]
##创建seurat对象
##创建STAR的Seurat对象
> STAR <- CreateSeuratObject(STARmatrix,project = "zsz",
min.cells = 3, min.features = 200)
##创建Cellranger的Seurat对象
> dir="cellranger/"
> counts <- Read10X(data.dir = dir)
> RANGER <- CreateSeuratObject(counts, project = "zsz",
min.cells=3, min.features = 200)
数据比较
> ##数据比较
> dim(STARmatrix)
[1] 33538 2048
> dim(counts)
[1] 33538 2112
> dim(STAR)
[1] 13350 2046
> dim(RANGER)
[1] 13314 2105
> fivenum(apply(STARmatrix,1,function(x) sum(x>0)))
MIR1302-2HG IGFBPL1 AL008723.2 FXYD7 MT-CO1
0 0 0 29 2047
> fivenum(apply(counts,1,function(x) sum(x>0)))
MIR1302-2HG FAM221B LINC01638 DGKE MT-CO1
0 0 0 29 2108
> pdf("box.pdf",height = 9,width = 9)
> boxplot(apply(STARmatrix,1,function(x) sum(x>0) ),main = "STAR",col = "lightgray")
> boxplot(apply(counts,1,function(x) sum(x>0) ),main = "Cellranger",col = "lightgray")
> dev.off()
RStudioGD
2
> pdf("hist.pdf",height = 9,width = 9)
> hist(apply(STARmatrix,2,function(x) sum(x>0) ),col = "lightgray",
+ breaks=20,xlim=c(0,4000),ylim=c(0,800),
+ labels=F,main="STAR",
+ xlab="genes",ylab="cells")
> abline(v=median(apply(STARmatrix,2,function(x) sum(x>0))),col='red')
> hist(apply(counts,2,function(x) sum(x>0) ),col = "lightgray",
+ breaks=20,xlim=c(0,4000),ylim=c(0,800),
+ labels=F,main="Cellranger",
+ xlab="genes",ylab="cells")
> abline(v=median(apply(counts,2,function(x) sum(x>0))),col='red')
> dev.off()
根据箱线图，直方图和矩阵的基本信息，可以看到STAR与cellranger的结果差距很小
质量控制与聚类比较
> pdf("qc.pdf",height = 9,width = 9)
> VlnPlot(STAR,
+ features = c("nFeature_RNA", "nCount_RNA"),
+ pt.size = 0.1,
+ ncol = 2)
> VlnPlot(RANGER,
+ features = c("nFeature_RNA", "nCount_RNA"),
+ pt.size = 0.1,
+ ncol = 2)
> dev.off()
RStudioGD
2
> STAR<-subset(STAR,subset=nFeature_RNA>500 & nFeature_RNA<2000)
> STAR <- NormalizeData(STAR, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> table(STAR@meta.data$orig.ident)
zsz
1896
> STAR <- FindVariableFeatures(STAR, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> top10 <- head(VariableFeatures(STAR), 10)
> top10
[1] "S100A9" "S100A8" "IGLC2" "IGKC" "LYZ" "IGLC3" "CCL3" "NFKBIA" "PTGDS" "S100A12"
> RANGER<-subset(RANGER,subset=nFeature_RNA>500 & nFeature_RNA<2000)
> RANGER<- NormalizeData(RANGER, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> table(RANGER@meta.data$orig.ident)
zsz
1895
> RANGER <- FindVariableFeatures(RANGER, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> top10 <- head(VariableFeatures(RANGER), 10)
> top10
[1] "S100A9" "S100A8" "IGLC2" "IGKC" "LYZ" "CCL3" "NFKBIA" "IGLC3" "PTGDS" "S100A12"
> scale.genes <- rownames(STAR)
> STAR <- ScaleData(STAR, features = scale.genes)
> STAR <- RunPCA(STAR, features = VariableFeatures(STAR))
> plot1 <- ElbowPlot(STAR, ndims=30, reduction="pca")
> scale.genes <- rownames(RANGER)
> RANGER <- ScaleData(RANGER, features = scale.genes)
> RANGER <- RunPCA(RANGER, features = VariableFeatures(RANGER))
> plot2 <- ElbowPlot(RANGER, ndims=30, reduction="pca")
> plotc <- plot1+plot2
> ggsave("pca.pdf", plot = plotc, width = 8, height = 4)
> STAR <- FindNeighbors(STAR, dims = 1:10)
> STAR <- FindClusters(STAR, resolution = 0.8)
> table(STAR@meta.data$seurat_clusters)
0 1 2 3 4 5 6 7 8 9
368 310 290 217 184 129 119 112 92 75
> metadata <- STAR@meta.data
> cell_cluster <-data.frame(cell_ID=rownames(metadata),
cluster_ID=metadata$seurat_clusters)
> STAR <- RunUMAP(STAR, dims = 1:20)
> embed_tsne <- Embeddings(STAR, 'umap')
> plot1 = DimPlot(STAR, reduction = "umap" ,label = "T", pt.size = 1,label.size = 4)
> RANGER <- FindNeighbors(RANGER, dims = 1:10)
> RANGER <- FindClusters(RANGER, resolution = 0.8)
> table(RANGER@meta.data$seurat_clusters)
0 1 2 3 4 5 6 7 8 9
375 300 290 212 195 128 122 108 91 74
> metadata <- RANGER@meta.data
> cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
> RANGER <- RunUMAP(RANGER,n.neighbors = 30,dims = 1:20)
> embed_umap <- Embeddings(RANGER, 'umap')
> plot2 = DimPlot(RANGER, reduction = "umap" ,label = "T", pt.size = 1,label.size = 4)
> plotc <- plot1+plot2
> ggsave("umap.pdf", plot = plotc, width = 8, height = 4)
可以看到两种分析方法umap与tsne聚类效果不太相同，但是基本的聚类与分群是一致的
结语
就结果而言，两种分析方法的结果不完全相同，但也基本一致，本次笔记中用到的数据和代码已上传github，在ScRNAseq_code/compare文件夹下，大家需要的可以下载试试
转载请注明：周小钊的博客>>>单细胞实战(4)：STAR与cellranger结果比较
