胃癌单细胞数据集GSE163558复现（一）：数据下载、整理与读取

前言

Hello小伙伴们大家好，我是生信技能树的小学徒”我才不吃蛋黄“。接下来的一段时间里，将由我开启一个新的学徒分享系列，给大家系统整理单细胞测序的代码。此系列包括但不限于以下内容： 数据下载与读取；质控和去批次；降维聚类；分群注释；差异分析；富集分析；拟时序分析；细胞通讯；CopyKAT 。

基本囊括了单细胞测序常规分析的所有方法，特别适合新手及入门者系统学习。因此，欢迎大家持续追更，关注公众号，点赞+评论+收藏+在看，您的鼓励将是我更新的动力，提前谢谢大家。

1.文献简介

本系列复现的文献题目为“Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing”。通讯作者是浙江大学的范骁辉教授，于2022年发表在Clin Transl Med杂志（IF=10.6）。

摘要：

目的：
通过单细胞测序分析揭示胃癌(GC)及其转移的生物学特性。

方法:
主要是收集了6例患者共10个新鲜组织标本（包括原发肿瘤、癌旁组织和不同器官或组织的转移瘤）进行了单细胞测序技术。并使用组织学分析和Bulk转录数据集进行了验证。

结果：

• 恶性上皮亚群细胞与侵袭特征、腹膜内转移倾向、上皮-间充质转化诱导的肿瘤干细胞表型或休眠样特征有关。基于TCGA胃癌队列的KM生存分析结果显示，前三个恶性上皮亚群相关基因的是GC患者预后的风险因子。

• 免疫和基质细胞表现出细胞异质性，并创造了亲肿瘤和免疫抑制的微环境。

• 基于淋巴来源的耗竭CD8+T细胞的20个基因signature可预测胃癌淋巴转移，该结果在并在TCGA队列中得到了验证。

• 除恶性肿瘤细胞，还发现了一个内皮细胞亚群、粘膜相关不变T细胞、T细胞样B细胞、浆细胞样树突状细胞、巨噬细胞、单核细胞和中性粒细胞可能有助于与 细胞毒/耗尽的CD8+T细胞 和/或 自然杀伤(NK)细胞相互作用（HLA-E-KLRC1/KLRC2），NKG2A(KLRC1)可能是胃癌治疗新的靶点。

• CD8+T细胞中PD-1的表达可预测胃癌患者对PD-1抑制剂的临床反应。

结论：
本研究对胃癌原发肿瘤和器官特异性转移的异质性微环境提供了深入的认识，为准确的诊断和治疗提供了支持。

以上便是本文的简介，接下来我们进入数据分析部分，开始下载并读取数据。

2.数据下载与整理

文章的单细胞测序数据存放于GEO数据库，编号为GSE163558（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163558）。本数据集收录了6例胃癌患者共10个新鲜组织标本（包括原发肿瘤、癌旁组织和不同器官或组织的转移瘤）：

可以看到，数据的文件格式为10X标准文件。10X标准文件包含cellranger上游比对分析产生的barcodes.tsv.gz，features.tsv.gz和matrix.mtx.gz 3个文件，分别代表细胞标签(barcode)、基因ID(feature)、表达数据（matrix）。

在读取文件之前，我们需要先加载R包：

rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(data.table)
library(dplyr)

下载数据之后，开始对原始文件进行处理，将原始文件分别整理为barcodes.tsv.gz，features.tsv.gz和matrix.mtx.gz并放到到各自的文件夹中。代码如下：

setwd("")  
dir='GSE163558_RAW/' 
fs=list.files('GSE163558_RAW/','^GSM')
fs
library(tidyverse)
samples=str_split(fs,'_',simplify = T)[,1]

##处理数据，将原始文件分别整理为barcodes.tsv.gz，features.tsv.gz和matrix.mtx.gz到各自的文件夹
#批量将文件名改为 Read10X()函数能够识别的名字
lapply(unique(samples),function(x){
  # x = unique(samples)[1]
  y=fs[grepl(x,fs)]
  folder=paste0("GSE163558_RAW/", paste(str_split(y[1],'_',simplify = T)[,1:2], collapse = "_"))
  dir.create(folder,recursive = T)
  #为每个样本创建子文件夹
  file.rename(paste0("GSE163558_RAW/",y[1]),file.path(folder,"barcodes.tsv.gz"))
  #重命名文件，并移动到相应的子文件夹里
  file.rename(paste0("GSE163558_RAW/",y[2]),file.path(folder,"features.tsv.gz"))
  file.rename(paste0("GSE163558_RAW/"




    
,y[3]),file.path(folder,"matrix.mtx.gz"))
})

3.数据读取与合并

整理好10X标准文件后，使用Read10X()函数对这三个标准文件进行整合，得到稀疏表达矩阵（行为基因、列为细胞，dgCMatrix格式）。在稀疏表达矩阵”tmp“的基础上，使用CreateSeuratObject函数构建Seurat对象。多个样本就需要对多个文件批量读取，在这里我们使用了lapply函数（亦可使用for循环）。

dir='GSE163558_RAW/'
samples=list.files( dir )
samples 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  tmp = Read10X(file.path(dir,pro )) 
  if(length(tmp)==2){
    ct = tmp[[1]] 
  }else{ct = tmp}
  sce =CreateSeuratObject(counts =  ct ,
                          project =  pro  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)#返回创建的Seurat对象，将其存储在sceList中。
}) 
View(sceList)

这里我们得到的sceList实际上包含了10个样本的Seurat对象，

查看其中一个：

PT1 <- sceList[1]
View(PT1)

可以看到，这是一个包含各类元素的Seurat V5对象（V5版本的assays对象下面多出了layers的结构）（详情见之前推文 https://mp.weixin.qq.com/s/2dtIS1qd0tPM1dQRCKptAQ ）。

接着，我们需要使用Seurat包的merge函数，将十个Seurat合并成一个Seurat对象。

do.call(rbind,lapply(sceList, dim))
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids = samples  ) 
names(sce.all@assays$RNA@layers)

此时，虽然我们已经完成了Seurat对象的创建，但是可以看到，十个样本仍然有10个layers。如果不进一步处理，后续在提取counts时数据不完整，分析会一直出错。因此我们需要使用JoinLayers函数对layers进行合并。

sce.all[["RNA"]]$counts 
# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
#看看合并前后的sce变化
sce.all
sce.all all)
sce.all

这是整合之前的：

这是整合之后的：

可以看到在合并Seurat和layers后，终于得到了一个完整的Seurat对象”sce.all“。我们可以查看”sce.all“内部的一些信息，以此来检查数据是否完整。

dim(sce.all[["RNA"]]$counts )
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head([email protected], 10)
table(sce.all$orig.ident) 
length(sce.all$orig.ident)

我们可以查看每个样本的细胞数量及总的细胞数量：

4.添加meta.data分组信息

在成功构建Seurat对象”sce.all“后，我们还需要给样本添加meta.data分组信息，以便后续做不同分组之间的对比以及提取亚组后进行进一步分析。
首先，我们查看现有的meta.data信息有哪些：

library(stringr)
phe = [email protected]
table(phe$orig.ident)
View(phe)

可以看到，现有的信息仅有：orig.ident(样品名),nCount_RNA(每个细胞的UMI数目)和nFeature_RNA（每个细胞中检测到的基因数量）。但原文中包含有样本的患者来源，组织来源、转移部位等信息，这些信息可以通过样本编号进行区分。因此我们可以利用文本处理函数”str_split“、”gsub“对患者编号进行处理，并添加以上信息到meta.data。

函数 str_split 用于拆分字符串：

phe$group = str_split(phe$orig.ident,'[_]',simplify = T)[,2] 

添加转移部位分组信息

phe$sample = phe$orig.ident
phe$sample = gsub("GSM\\d+_PT\\d+", "GC", phe$sample) 
phe$sample = gsub("GSM\\d+_LN\\d+", "GC_lymph_metastasis", phe$sample) 
phe$sample = gsub("GSM\\d+_O1", "GC_ovary_metastasis", phe$sample) 
phe$sample = gsub("GSM\\d+_P1", "GC_peritoneum_metastasis", phe$sample) 
phe$sample = gsub("GSM\\d+_Li\\d+", "GC_liver_metastasis", phe$sample) 
phe$sample = gsub("GSM\\d+_NT\\d+", "adjacent_nontumor", phe$sample) 
table(phe$sample)

添加患者来源信息

phe$patient = phe$orig.ident
table(phe$patient)
phe$patient = gsub("GSM5004180_PT1|GSM5004188_Li1", "Patient1", phe$patient) 
phe$patient = gsub("GSM5004181_PT2|GSM5004183_NT1", "Patient2", phe$patient) 
phe$patient = gsub("GSM5004186_O1", "Patient3", phe$patient) 
phe$patient