发布日期:2024-11-06 07:19 点击次数:167
图片软件开发资讯软件开发资讯
引子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')fslibrary(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)。
app接着,软件开发价格咱们需要使用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 resultLayerData(sce.all, assay = "RNA", layer = "counts")#望望并吞前后的sce变化sce.allsce.all <- JoinLayers(sce.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(sce.all@meta.data, 10)table(sce.all$orig.ident) length(sce.all$orig.ident)
咱们不错检讨每个样本的细胞数量及总的细胞数量:
图片
4.添加meta.data分组信息在收效构建Seurat对象”sce.all“后,咱们还需要给样本添加meta.data分组信息,以便后续作念不同分组之间的对比以及索取亚组后进行进一步分析。当先,咱们检讨现存的meta.data信息有哪些:
library(stringr)phe = sce.all@meta.datatable(phe$orig.ident)View(phe)
图片
不错看到,现存的信息仅有:orig.ident(样品名),nCount_RNA(每个细胞的UMI数量)和nFeature_RNA(每个细胞中检测到的基因数量)。但原文中包含有样本的患者开头,组织开头、篡改部位等信息,这些信息不错通过样本编号进行分裂。因此咱们不错垄断文本措置函数”str_split“、”gsub“对患者编号进行措置,并添加以上信息到meta.data。
龙头分析:历史同期第181期龙头分别开出号码:01→01→05,龙头开出比较密集,去年同期龙头号码上升了4个点位,与去年龙头相比,今年第181期预计龙头转向下降,参考号码04。
函数 str_split 用于拆分字符串:
phe$group = str_split(phe$orig.ident,'[_]',simplify = T)[,2]
添加篡改部位分组信息
phe$sample = phe$orig.identphe$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.identtable(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 = gsub("GSM5004185_LN2", "Patient5", phe$patient) phe$patient = gsub("GSM5004187_P1", "Patient6", phe$patient) phe$patient = gsub("GSM\\S+", "Patient4", phe$patient) table(phe$patient)
文中还有一类分群为宽泛NT,原发PT,篡改M,用ifelse函数添加
phe$tissue <- ifelse(phe$orig.ident %in% c("GSM5004180_PT1","GSM5004181_PT2","GSM5004182_PT3"),"PT", ifelse(phe$orig.ident %in% c("GSM5004183_NT1"),"NT","M"))table(phe$tissue)sce.all@meta.data = pheView(phe)
图片
临了不错看到,咱们收效添加了系数的样本信息。至此,数据下载、整理、读取扫尾,接下来不错运转走卑劣的Seurat V5模范进程。结语本期咱们对文献选录进行了简要回首,下载了GSE163558胃癌数据集10个样本的10X形貌的单细胞测序数据,并对文献进行了整理,在批量读取了10X文献后,进行了并吞并收效构建Seurat对象,在此基础上将患者的临床信息添加到meta.data分组信息中。下一期,咱们将在此基础上走Seurat V5模范进程,对Seurat对象进行QC质控、并垄断harmony整合去批次、并按模范进程进行降维聚类分群。
图片
本站仅提供存储奇迹,系数内容均由用户发布,如发现存害或侵权内容,请点击举报。