生信学习宝藏库单细胞学习

Seurat包基本分析实战

2020-10-02  本文已影响0人  小贝学生信

一、背景知识

二、大致思路

(1)根据GEO号下载表达矩阵,以及meta信息;
(2)根据meta信息筛选tumor cell 的表达矩阵;
(3)根据表达矩阵构建Seurat对象,以及添加分组信息、质控;
(4)归一化→降维→聚类;
(5)组图

三、具体实现

1、下载数据

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465

1-1
1-2
1-3

2、导入数据

a <- read.table("GSE84465_GBM_All_data.csv.gz")
a[1:4,1:4]
dim(a) 
#[1] 23465  3589
#原始矩阵有3598个cell,23465个gene
b <- read.table("SraRunTable.txt",sep = ",",header = T)
dim(b)
b[1:4,1:4]
#[1] 3589   33

3、筛选目标类型cell

new.b <- data.frame(names=paste0("X",b$plate_id,".",b$Well),
#names即为a的colname
                    tissue=b$TISSUE, #tumor二分类
                    group=b$Patient_ID) #样本分类

new.b <- new.b[match(colnames(a),new.b$names),]
#使顺序一致,以便后期的筛选
table(new.b$names==colnames(a))
rownames(new.b) <- colnames(a)

index <- new.b$tissue=="Tumor"
a_filt <- a[,index]
dim(a_filt)
#挑选到2343个tumor cell 表达矩阵

b$Patient_ID有四类,说明这些细胞是从四个病人的病灶组织中取的。

4、构建Seurat对象,以及质控、可视化

group=data.frame(group=new.b[index,3],
                 row.names = colnames(a_filt))

library("Seurat")
scRNA = CreateSeuratObject(counts=a_filt,
                           meta.data = group,
                           min.cells = 3, min.features = 50)
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
#标记线粒体基因
head(scRNA@meta.data)
summary(scRNA@meta.data)
#结果显示没有线粒体基因,因此感觉过滤也没有意义
pctMT=5
scRNA <- subset(scRNA, subset = percent.mt < pctMT)
4-1

可视化三个图

#图A:观察不同组cell的counts、feature分布
col.num <- length(levels(scRNA@meta.data$group))
p1_1 <- VlnPlot(scRNA,
        features = c("nFeature_RNA","nCount_RNA"),
        group.by = "group",
        cols =rainbow(col.num)) +
  theme(legend.position = "none") +
  labs(tag = "A")
#图B:nCount_RNA与对应的nFeature_RNA关系
library(ggplot2)
p1_2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",
               group.by = "group",pt.size = 1.3) +
  labs(tag = "B")

5、挑选hvg高变基因

#标准化
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
#归一化
scRNA <- ScaleData(scRNA, features = (rownames(scRNA)))
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 1500) 
top10 <- head(VariableFeatures(scRNA), 10) 
plot1 <- VariableFeaturePlot(scRNA) 
p1_3 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) +
  theme(legend.position = c(0.1,0.8)) +
  labs(tag = "C")
#标记top10 hvg
library(patchwork)
p1_1 | p1_2 | p1_3
上三个子图

6、降维

#默认使用前2000个hvg降维处理
scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 
p2_1 <- DimPlot(scRNA, reduction = "pca", group.by="group")+
  labs(tag = "D")

#挑选主成分
scRNA <- JackStraw(scRNA,reduction = "pca")
scRNA <- ScoreJackStraw(scRNA,dims = 1:20)
p2_2 <- JackStrawPlot(scRNA,dims = 1:20, reduction = "pca") +
  theme(legend.position="bottom") +
  labs(tag = "E")
p2_3 <- ElbowPlot(scRNA, ndims=20, reduction="pca") 
#结果显示可挑选前20个pc
p2_1| (p2_2 | p2_3)
中三个子图

7、聚类

pc.num=1:20
scRNA <- FindNeighbors(scRNA, dims = pc.num) 
# dims参数,需要指定哪些pc轴用于分析;这里利用上面的分析,选择20
scRNA <- FindClusters(scRNA, resolution = 0.4)
table(scRNA@meta.data$seurat_clusters)
#  0   1   2   3   4   5   6   7   8   9  10  11  12 
#393 380 281 271 177 140 122 108  72  51  41  35  18 
scRNA = RunTSNE(scRNA, dims = pc.num)
p3_1 <- DimPlot(scRNA, reduction = "tsne") +
  labs(tag = "E")
diff.wilcox = FindAllMarkers(scRNA)
#大概2-3min
head(diff.wilcox)
library(tidyverse)
all.markers = diff.wilcox %>% select(gene, everything()) %>%
  subset(p_val<0.05 & abs(diff.wilcox$avg_logFC) > 0.5)
top20 = all.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
top20
#0-9共10个cluster,每个cluster选取top10高变基因
top20 = CaseMatch(search = as.vector(top20$gene), match = rownames(scRNA)) 
length(top20)
length(unique(sort(top20)))
#这里选取的是wilcox方法挑选的差异基因
p3_2 <- DoHeatmap(scRNA, features = top20, group.by = "seurat_clusters")

此处的结果也是与原文差别比较大的地方。均是对每个clust寻找top20 marker gene。但是原文使用的limma包识别,去重后仅有96个gene,而我自己尝试的或还有227个,相差比较大。
The differential analysis identified 8,025 marker genes. The top 20 marker genes of each cell cluster are displayed in the heatmap. A total of 96 genes are listed beside of the heatmap after omitting the same top marker genes among clusters. The colors from purple to yellow indicate the gene expression levels from low to high.

p3_1 | p3_2
下二个子图

8、总图合并

p <- (p1_1 | (p1_2 | p1_3) ) /
  ((p2_1| p2_2 | p2_3) /
     (p3_1 | p3_2))
ggsave("my_try.pdf", plot = p, width = 15, height = 18) 
my_try.pdf
上一篇下一篇

猜你喜欢

热点阅读