轻松一挖就节约10万经费

2022-06-08 20:07:06 浏览数 (1)

大家的国自然都不是大风刮来的，节约纳税人的钱也是咱们科研工作者的基本素质。最就看到朋友圈有人分享了发表在《Journal for ImmunoTherapy of Cancer》杂志的文章：《Genetic characteristics involving the PD-1/PD-L1/L2 and CD73/A2aR axes and the immunosuppressive microenvironment in DLBCL》，虽然写的是使用了多种测序技术，如下所示：

  • 全外显子测序
  • 靶向深度测序
  • RNA-Seq
  • 单细胞转录组测序

但是实际上它里面的单细胞转录组测序 是来源于公共数据集GSE182434的4例病人样本（DLBCL002 、DLBCL007、 DLBCL008、DLBCL111），也就是说他们省下来了4个病人的肿瘤单细胞转录组费用，哪怕是按照一两年前的均价2.5万的单个10x费用，也算是省下来了10万经费！

省下来了10万经费

这个公共数据集GSE182434的单细胞表达量矩阵和其细胞类型注释都是公开的，所以很容易挑选子集。首先自己去下载公共数据集GSE182434的两个txt.gz文件，存放在data文件夹里面，然后复制粘贴下面的代码，如下所示：

代码语言：javascript复制
rm(list=ls())
options(stringsAsFactors = F)
getwd()
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)

#倒入数据
#直接是一个matix数据。
library(data.table) 
dat=fread( "data/GSE182434_raw_count_matrix.txt.gz",data.table = F)  
dim(dat) 
dat[1:4,1:4]
rownames(dat)=dat[,1]
dat=dat[,-1]
dat[1:4,1:4] 
annotation = fread( "data/GSE182434_cell_annotation.txt.gz",data.table = F)
annotation[1:4,1:4]
table(annotation$Patient,annotation$CellType) 
#先筛选一下数据，在创建sec对象
#筛选条件①、四个样本 ：DLBCL002 、DLBCL007、 DLBCL008、DLBCL111
#筛选条件②、CD8 细胞
sample = c("DLBCL002","DLBCL007","DLBCL008","DLBCL111")
CellType = "T cells CD8"
annotation = annotation[annotation$Patient %in% sample,]
annotation = annotation[annotation$CellType %in% CellType,]
dim(annotation)  #4个样本，一共5千个CD8细胞
dat = dat[,annotation$ID]
dim(dat) 
#创建Seurat object
library(Seurat)
sce <- CreateSeuratObject(counts = dat, 
                          project = "sce", 
                          min.cells = 3, 
                          min.features = 200)
sce 
# 这里省略了把  annotation 信息添加到 sce 里面的代码，我相信你肯定是可以写出来。

这样就构建好了自己的单细胞转录组seurat对象了，接下来就是对这个常规的降维聚类分群。因为是来源于4个病人的数据，所以需要harmony处理一下，也可以看处理前后的效果哦：

代码语言：javascript复制

#简单画一个DimPlot，判断是否要harmony
sce = NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce = FindVariableFeatures(sce)
sce = ScaleData(sce, vars.to.regress = c("nFeature_RNA", "percent_mito"))
sce = RunPCA(sce, npcs = 20)
sce = RunTSNE(sce, npcs = 20)
sce = Runtsne(sce, dims = 1:10)
sce = RunUMAP(sce, dims = 1:10)
sce = FindNeighbors(sce, dims = 1:20, k.param = 60, prune.SNN = 1/15)  #这个dims要和PCA的npcs对应起来。
DimPlot(sce,reduction = "umap",label=T,group.by = "sample") 
ggsave(filename="3-cluster/no_harmony_DimPlot.pdf")


#需要 harmony 去批次 因为是主动筛选4个样本出来的
library(harmony)
sce <- RunHarmony(sce, "sample")
names(sce@reductions)
sce = RunUMAP(sce,  dims = 1:15, reduction = "harmony")
sce = RunTSNE(sce, npcs = 20,reduction = "harmony")
sce = FindNeighbors(sce, reduction = "harmony", dims = 1:15) 
DimPlot(sce,reduction = "umap",label=T,group.by = "sample") 
ggsave(filename="3-cluster/harmony_DimPlot.pdf")


###设置分辨率
#分辨率分类
#文中是分了6群
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1,0.15, 0.2, 0.3, 0.4,0.45,0.5,0.8,0.9,1)) {
  sce=FindClusters(sce, resolution = res, algorithm = 1)
}
colnames(sce@meta.data)
apply(sce@meta.data[,grep("RNA_snn_res",colnames(sce@meta.data))],2,table)

#0.8
## (2)可视化高低分辨率的分群情况
p1_dim = plot_grid(ncol = 4, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.01")   
                     ggtitle("louvain_0.01"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.05")   
                     ggtitle("louvain_0.05"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.1")   
                     ggtitle("louvain_0.1"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.2")   
                     ggtitle("louvain_0.2"))
p1_dim
ggsave(plot = p1_dim, filename="3-cluster/Dimplot_diff_resolution_low.pdf",width = 14,height = 3)
p1_dim = plot_grid(ncol = 5, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.3")   
                     ggtitle("louvain_0.3"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.5")   
                     ggtitle("louvain_0.5"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.8")   
                     ggtitle("louvain_0.8"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.9")   
                     ggtitle("louvain_0.9"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.1")   
                     ggtitle("louvain_1"))
p1_dim
ggsave(plot = p1_dim, filename = "3-cluster/Dimplot_diff_resolution_high.pdf",width = 17.5,height = 3)
## (3)聚类树
p2_tree = clustree(sce@meta.data, prefix = "RNA_snn_res.")
p2_tree
ggsave(plot = p2_tree, filename="3-cluster/Tree_diff_resolution.pdf",width = 10,height = 11)

#### 4. 因为文章中聚类数为7）,所以接下来分析，按照分辨率为0.4
colnames(sce@meta.data)
sel.clust = "RNA_snn_res.0.4"
sce <- SetIdent(sce, value = sel.clust)
table(sce@active.ident)

最基础的往往是降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释

后面的一些可视化（一些基因的表达量小提琴图，细胞比例条形图），我们就直接看学徒的图表复现吧，代码真的是超级简单。

下面是学徒的笔记

1、原文文献：

《Genetic characteristics involving the PD-1/PD-L1/L2 and CD73/A2aR axes and the immunosuppressive microenvironment in DLBCL》

2、数据集：GSE182434

（其中4例样本（DLBCL002 、DLBCL007、 DLBCL008、DLBCL111）

3、研究目的：

研究，在功能失调型的CD8 T细胞里面 ，PD-1和A2aR这两个基因的表达情况。

以确定：PD-1和A2aR的表达，是否会影响到CD8 T的功能正常表达

4、分析和复现流程：

①、去批次：

因为是从整个数据集中筛选4个样本出来，人为选定的。所以需要去掉批次。

可以看到。是需要去批次的

没有去批次

harmony去批次的

②、分群：

最后分群选择的resolution=0.4 。分为7群。

原文

复现

③、确定功能失调的CD8 细胞亚群

利用CD8A、GZMB、CTLA4、TIGIT、LAG3 这5个基因的表达情况来确定

可以确定的是：原文的1群，就是复现的0群.也就是功能失调的CD8 T细胞

原文

复现

④、比例图

可以看出，虽然比例上是有差异的。

但是总体趋势是一样的。

都是DLBCL008这个样本的功能失调的CD8 T细胞 占的比例最高

原文

复现

⑤、PD-1、A2aR 的表达情况

也是和原文趋势一致。

结论就是：

DLBCL002在功能失调型的CD8 T细胞不表达PD-1和A2aR，

DLBCL007和DLBCL111仅表达PD-1，

DLBCL008同时表达PD-1和A2aR，其含有大量的功能失调的CD8 T细胞

原文

复现

写在文末

本文关心的是功能失调的CD8 T细胞，其实T细胞耗竭也是有多个阶段，比如文章：《Developmental Relationships of Four Exhausted CD8 T Cell Subsets Reveals Underlying Transcriptional and Epigenetic Landscape Control Mechanisms》提到的 ：

代码语言：javascript复制
 TexProg1 = c("Tcf7","Myb","Il7r","Sell","Cxcr5","Icos","Cd28","ISGs","Irf7","Oas1","Stat1")
 TexProg2 = c("Mki67","Anxa2d","Itgb7","Alcam","Top2a") #细胞周期 
 TexInt = c("Grzma","Grzmb","Prf1","Klrk1","Cx3cr1","Tbx21","Zeb2","Id2","Prdm1")# 细胞毒效应相关基因
 TexTerm = c("Pdcd1","Lag3","Tigit","Cd244","Entpd1","Cd101","Cd38","Zap70","Nfatc1") # 抑制性受体

学徒作业

去下载公共数据集GSE182434的两个txt.gz文件，同样的使用我们的代码筛选例病人样本（DLBCL002 、DLBCL007、 DLBCL008、DLBCL111），走我们的降维聚类分群代码，看看我这上面列出来的Four Exhausted CD8 T Cell Subsets相关基因是否有规律。

axes txt

0 人点赞