轻松一挖就节约10万经费
大家的国自然都不是大风刮来的，节约纳税人的钱也是咱们科研工作者的基本素质。最就看到朋友圈有人分享了发表在《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文件夹里面，然后复制粘贴下面的代码，如下所示：
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处理一下，也可以看处理前后的效果哦：
#简单画一个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》提到的 ：
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相关基因是否有规律。
腾讯云开发者
