来吧综艺网

图片快播三级电影快播三级电影

弁言

Hello小伙伴们宇宙好，我是生信手段树的小学徒”我才不吃蛋黄“。继胃癌单细胞复现系列罢了之后，经过了顿然的休整，我又转头了。之前每一期的推文，我齐有追究准备。推文发出之后，也仔细阅读了每一条留言。至极感谢宇宙的饱读舞。

接下来的一段时间里，我将再次开启一个新的学徒共享系列，给宇宙系统整理肺腺癌单细胞测序的代码。

著作的具体内容，二由真诚已在生信菜鸟团公众号发布，请宇宙移步阅读：单细胞测序+空间转录组描述从癌前病变到浸润性肺腺癌的动态演变。

本系列包括但不限于以下内容：数据下载与读取；质控和去批次；降维聚类；分群谛视；相反分析；富集分析；亚群细分；inferCNV;拟时序分析；细胞通信；SCENIC转录因子分析等。

为了加深对代码的调治，除了阅读推文，不错加入单细胞月更群，也不错相助不雅看生信手段树B站系列造就视频。

图片

如若你是的单细胞外行，推选报名干预生信马拉松课程（生信初学&数据挖掘线上直播课10月班），让你快速掌执生信数据挖掘手段。

本系列推文展望12篇傍边，宽恕宇宙不竭追更，包涵公众号，点赞+指摘+保藏+在看，您的饱读舞将是我更新的能源，同期宽恕宇宙品评指正。

数据搞定经由

最初铲除系统环境变量，加载R包：

rm(list=ls())options(stringsAsFactors = F) getwd()setwd("../GSE189357-LUAD-scRNA-ST")source('scRNA_scripts/lib.R')library(Seurat)library(ggplot2)library(clustree)library(cowplot)library(data.table)library(dplyr)

当我们在建造数据框的时候，R谈话将会默许把字符型（character）当成因子（factor）。stringsAsFactors = F意味着，“在读入数据时，遭逢字符串之后，不将其改换为factors，仍然保留为字符串样式”。

在R谈话中，stringsAsFactors是一个全局选项，它截止着当读取数据时，字符串类型的列是否自动改换为因子(factors)类型。因子类型在R中是一种特殊的数据类型，用于默示分类数据。

在R的早期版块中，默许情况下，读取数据时字符串会被改换为因子。然则，因子类型有其特定的属性和活动，这在某些情况下可能不是用户所期许的。举例，因子类型会将数据视为分类变量，何况会存储数据的独一值动作里面整数编码，这可能会增多内存使用量，何况在进行某些类型的数据分析时可能会引起浑浊。

在单细胞数据分析中，数据经常包含广宽的基因抒发信息，这些信息是以字符串形状存在的。在进行数据分析之前，莳植stringsAsFactors = F的观念经常是为了：

幸免无用要的内存遽然：将字符串改换为因子类型会增多内存的使用，因为因子需要存储一个里面的整数编码和对应的独一值列表。在搞定大边界的单细胞数据集时，这可能会成为一个问题。

保持数据的一致性：在数据分析过程中，保持数据类型的一致性是很紧迫的。如若某些操作期许数据是字符串类型，而数据被自动改换为因子，这可能会导致造作或者不一致的末端。

擢升代码的可读性和可儿戴性：明确地莳植stringsAsFactors = F不错让代码的读者知说念数据不会被自动改换为因子，这有助于调治代码的意图和活动。

幸免因子类型带来的潜在问题：因子类型在某些统计分析中可能会引起问题，比如在进行机器学习或者数据挖掘时，因子类型的数据可能需要特殊的搞定技艺被正确使用。

因此，在读取单细胞数据之前莳植stringsAsFactors = F是一种风雅的编程践诺，它有助于确保数据以正确的形状被搞定，何况不错减少后续分析中可能出现的问题。

step1:导入数据

数据存储在“GSE189357_RAW”文献夹中：

dir='GSE189357_RAW/' fs=list.files('GSE189357_RAW/'，'^GSM')fslibrary(tidyverse)samples=str_split(fs，'_'，simplify = T)[，1]

搞定数据，将原始文献别离整理为barcodes.tsv.gz，features.tsv.gz和matrix.mtx.gz到各自的文献夹。批量将文献名改为 Read10X() 函数概况识别的名字：

if(F){lapply(unique(samples)，function(x){  # x = unique(samples)[1]  y=fs[grepl(x，fs)]  folder=paste0("GSE189357_RAW/"， paste(str_split(y[1]，'_'，simplify = T)[，1:2]， collapse = "_"))  dir.create(folder，recursive = T)  #为每个样本创建子文献夹  file.rename(paste0("GSE189357_RAW/"，y[1])，file.path(folder，"barcodes.tsv.gz"))  #重定名文献，并移动到相应的子文献夹里  file.rename(paste0("GSE189357_RAW/"，y[2])，file.path(folder，"features.tsv.gz"))  file.rename(paste0("GSE189357_RAW/"，y[3])，file.path(folder，"matrix.mtx.gz"))})}

自行创建函数function，内置Read10X()函数和CreateSeuratObject()函数，批量读取文献原始文献创建Seurat对象：

dir='GSE189357_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)}) 

这里我们赢得的sceList实际上包含了9个样本的Seurat对象，接着，我们需要使用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)#> [1] "counts.GSM5699777_TD1" "counts.GSM5699778_TD2" "counts.GSM5699779_TD3"#> [4] "counts.GSM5699780_TD4" "counts.GSM5699781_TD5" "counts.GSM5699782_TD6"#> [7] "counts.GSM5699783_TD7" "counts.GSM5699784_TD8" "counts.GSM5699785_TD9"

此时，诚然我们仍是完成了Seurat对象的创建，然则不错看到，九个样本仍然有9个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.alldim(sce.all[["RNA"]]$counts )

在归拢Seurat和layers后，终于赢得了一个完满的Seurat对象”sce.all“。我们不错稽查”sce.all“里面的一些信息，以此来查验数据是否完满。

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)# fivenum(sce.all$nFeature_RNA)# table(sce.all$nFeature_RNA>800) # sce.all=sce.all[，sce.all$nFeature_RNA>800]# sce.all

在胜利构建Seurat对象”sce.all“后，我们还需要给样本添加meta.data分组信息，以便后续作念不同分组之间的对比以及索要亚组后进行进一步分析。最初，我们稽查现存的meta.data信息有哪些：

library(stringr)phe = sce.all@meta.datatable(phe$orig.ident)phe$patient = str_split(phe$orig.ident，'[_]'，simplify = T)[，2]table(phe$patient)

按照患者源流，肺腺癌（LUAD）从原位腺癌（AIS）证明到微浸润性腺癌（MIA）和随后的浸润性腺癌（IAC）

phe$sample = phe$patienttable(phe$sample)phe$sample = gsub("TD5|TD7|TD8"， "AIS"， phe$sample) phe$sample = gsub("TD3|TD4|TD6"， "MIA"， phe$sample) phe$sample = gsub("TD1|TD2|TD9"， "IAC"， phe$sample) table(phe$sample)sce.all@meta.data = phe

step2: QC质控

质控（quality control， QC）的观念是为了去除质料较差细胞，低质料的细胞会酿成特有的亚群，使分群末端变得复杂；在主要素分析过程中，前几个主要要素将拿获质料相反而不是生物学相反，从而裁减降维成果。因此，低质料的细胞可能会导致下流分析出现误导性末端。为了幸免上述情况的发生，我们需要不才游分析启动前排裁撤这些低质料细胞。QC主要的标的有nCount_RNA(每个细胞的UMI数量)和nFeature_RNA（每个细胞中检测到的基因数量）以及"percent_mito"（默示细胞中线粒体基因的比例）这三个标的。此外，还不错纳入”percent_ribo“（核糖体基因比例）和”percent_hb“（红血细胞基因比例）两个标的。细胞的UMI分子数和基因数反应细胞的质料，数量太低可能是细胞碎屑，太高则可能是两个或多个细胞结团的情况；线粒体基因高抒发的细胞，可能是处于凋一火景色或者裂解景色的细胞；核糖体基因高抒发的细胞，细胞内出现RNA降解时；血红卵白基因高抒发的细胞经常为红细胞，而红细胞自己是不含有细胞核的，且捎带的基因少，其捎带的基因与疾病以及生物发育等过程莫得太大的相干，是以不错平直剔裁撤。在这里，我们要诓骗PercentageFeatureSet函数别离筹划每个细胞的"percent_mito"，”percent_ribo“和”percent_hb“，具体可见scripts代码

sp='human'dir.create("./1-QC")setwd("./1-QC")# 如若过滤的太狠，就需要去修改这个过滤代码source('../scRNA_scripts/qc.R')sce.all.filt = basic_qc(sce.all)print(dim(sce.all))print(dim(sce.all.filt))##细胞减少了极少setwd('../')getwd()fivenum(sce.all.filt$percent_ribo)table(sce.all.filt$nFeature_RNA> 5)

qc函数如下：

basic_qc <- function(input_sce){  #筹划线粒体基因比例  mito_genes=rownames(input_sce)[grep("^MT-"， rownames(input_sce)，ignore.case = T)]   print(mito_genes) #可能是13个线粒体基因  #input_sce=PercentageFeatureSet(input_sce， "^MT-"， col.name = "percent_mito")  input_sce=PercentageFeatureSet(input_sce， features = mito_genes， col.name = "percent_mito")  fivenum(input_sce@meta.data$percent_mito)    #筹划核糖体基因比例  ribo_genes=rownames(input_sce)[grep("^Rp[sl]"， rownames(input_sce)，ignore.case = T)]  print(ribo_genes)  input_sce=PercentageFeatureSet(input_sce，  features = ribo_genes， col.name = "percent_ribo")  fivenum(input_sce@meta.data$percent_ribo)    #筹划红血细胞基因比例  Hb_genes=rownames(input_sce)[grep("^Hb[^(p)]"， rownames(input_sce)，ignore.case = T)]  print(Hb_genes)  input_sce=PercentageFeatureSet(input_sce，  features = Hb_genes，col.name = "percent_hb")  fivenum(input_sce@meta.data$percent_hb)    #可视化细胞的上述比例情况  feats <- c("nFeature_RNA"， "nCount_RNA"， "percent_mito"， "percent_ribo"， "percent_hb")  feats <- c("nFeature_RNA"， "nCount_RNA")  p1=VlnPlot(input_sce， group.by = "orig.ident"， features = feats， pt.size = 0， ncol = 2) +  NoLegend()  p1   w=length(unique(input_sce$orig.ident))/3+5;w  ggsave(filename="Vlnplot1.pdf"，plot=p1，width = w，height = 5)    feats <- c("percent_mito"， "percent_ribo"， "percent_hb")  p2=VlnPlot(input_sce， group.by = "orig.ident"， features = feats， pt.size = 0， ncol = 3， same.y.lims=T) +  scale_y_continuous(breaks=seq(0， 100， 5)) + NoLegend()  p2   w=length(unique(input_sce$orig.ident))/2+5;w  ggsave(filename="Vlnplot2.pdf"，plot=p2，width = w，height = 5)    p3=FeatureScatter(input_sce， "nCount_RNA"， "nFeature_RNA"， group.by = "orig.ident"， pt.size = 0.5)  ggsave(filename="Scatterplot.pdf"，plot=p3)    #字据上述标的，过滤低质料细胞/基因  #过滤标的1:最少抒发基因数的细胞&最少抒发细胞数的基因  # 一般来说，在CreateSeuratObject的时候仍是是进行了这个过滤操作  # 如若后期看到了我方的单细胞降维聚类分群末端很诡异，就不错回偏激来看质料截留步调  # 先走默许经由即可  if(F){ selected_c <- WhichCells(input_sce， expression = nFeature_RNA > 500) selected_f <- rownames(input_sce)[Matrix::rowSums(input_sce@assays$RNA$counts > 0 ) > 3] input_sce.filt <- subset(input_sce， features = selected_f， cells = selected_c) dim(input_sce)  dim(input_sce.filt)   }    input_sce.filt =  input_sce    # par(mar = c(4， 8， 2， 1))  # 这里的C 这个矩阵，有极少大，不错商酌随抽样   C=subset(input_sce.filt，downsample=100)@assays$RNA$counts  dim(C)  C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100  most_expressed <- order(apply(C， 1， median)， decreasing = T)[50:1]    pdf("TOP50_most_expressed_gene.pdf"，width=14)  boxplot(as.matrix(Matrix::t(C[most_expressed， ]))， cex = 0.1， las = 1，  xlab = "% total count per cell"，  col = (scales::hue_pal())(50)[50:1]，  horizontal = TRUE)  dev.off()  rm(C)    #过滤标的2:线粒体/核糖体基因比例(字据上头的violin图)  selected_mito <- WhichCells(input_sce.filt， expression = percent_mito < 25)  selected_ribo <- WhichCells(input_sce.filt， expression = percent_ribo > 3)  selected_hb <- WhichCells(input_sce.filt， expression = percent_hb < 1 )  length(selected_hb)  length(selected_ribo)  length(selected_mito)    input_sce.filt <- subset(input_sce.filt， cells = selected_mito)  #input_sce.filt <- subset(input_sce.filt， cells = selected_ribo)  input_sce.filt <- subset(input_sce.filt， cells = selected_hb)  dim(input_sce.filt)    table(input_sce.filt$orig.ident)     #可视化过滤后的情况  feats <- c("nFeature_RNA"， "nCount_RNA")  p1_filtered=VlnPlot(input_sce.filt， group.by = "orig.ident"， features = feats， pt.size = 0， ncol = 2) +  NoLegend()  w=length(unique(input_sce.filt$orig.ident))/3+5;w   ggsave(filename="Vlnplot1_filtered.pdf"，plot=p1_filtered，width = w，height = 5)    feats <- c("percent_mito"， "percent_ribo"， "percent_hb")  p2_filtered=VlnPlot(input_sce.filt， group.by = "orig.ident"， features = feats， pt.size = 0， ncol = 3) +  NoLegend()  w=length(unique(input_sce.filt$orig.ident))/2+5;w   ggsave(filename="Vlnplot2_filtered.pdf"，plot=p2_filtered，width = w，height = 5)   return(input_sce.filt)   }

可视化

图片

Vlnplot1

图片

Vlnplot2

图片

Scatterplot

图片

TOP50_most_expressed_gene

图片

Vlnplot1_filtered

图片

Vlnplot2_filteredstep3: harmony整合多个单细胞样品

细胞筛选之后，我们需要进行harmony整合。第一期我们在创建总的Seurat对象时，使用了merge函数对多个Seurat进行了浅近的归拢。merge仅仅按照行和列进行了归拢，并未对数据进行其他搞定。

在拿到下流单细胞矩阵前，样本阅历了多个实验步调的搞定，时间、搞定东说念主员、试剂以及技艺平台等变量齐会成为混杂因素。以上因素搀杂到一齐，就会导致数据产生批次效应（batch effect）。为了尽可能幸免混杂因素，我们不错严格把控测序的技艺经由，同期也需要不才游分析中进行过后挽救（算法去批次）。当今单细胞测序常用的去批次算法有Harmony，CCA，RPCA，FastMNN，scVI等。在这里，我们使用Harmony进行演示。

set.seed(10086)table(sce.all.filt$orig.ident)if(T){  dir.create("2-harmony")  getwd()  setwd("2-harmony")  source('../scRNA_scripts/harmony.R')  # 默许 ScaleData 莫得添加"nCount_RNA"， "nFeature_RNA"  # 默许的  sce.all.int = run_harmony(sce.all.filt)  setwd('../')}

harmony函数如下：

run_harmony <- function(input_sce){    print(dim(input_sce))  input_sce <- NormalizeData(input_sce，   normalization.method = "LogNormalize"，  scale.factor = 1e4)   input_sce <- FindVariableFeatures(input_sce)  input_sce <- ScaleData(input_sce)  input_sce <- RunPCA(input_sce， features = VariableFeatures(object = input_sce))  seuratObj <- RunHarmony(input_sce， "orig.ident")  names(seuratObj@reductions)  seuratObj <- RunUMAP(seuratObj，  dims = 1:15，   reduction = "harmony")    # p = DimPlot(seuratObj，reduction = "umap"，label=T )   # ggsave(filename='umap-by-orig.ident-after-harmony'，plot = p)    input_sce=seuratObj  input_sce <- FindNeighbors(input_sce， reduction = "harmony"，  dims = 1:15)   input_sce.all=input_sce    #莳植不同的分辨率，不雅察分群成果(禁受哪一个？)  for (res in c(0.01， 0.05， 0.1， 0.2， 0.3， 0.5，0.8，1)) { input_sce.all=FindClusters(input_sce.all， #graph.name = "CCA_snn"，  resolution = res， algorithm = 1)  }  colnames(input_sce.all@meta.data)  apply(input_sce.all@meta.data[，grep("RNA_snn"，colnames(input_sce.all@meta.data))]，2，table)    p1_dim=plot_grid(ncol = 3， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.0.01") +    ggtitle("louvain_0.01")， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.0.1") +    ggtitle("louvain_0.1")， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.0.2") +    ggtitle("louvain_0.2"))  ggsave(plot=p1_dim， filename="Dimplot_diff_resolution_low.pdf"，width = 14)    p1_dim=plot_grid(ncol = 3， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.0.8") +    ggtitle("louvain_0.8")， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.1") +    ggtitle("louvain_1")， DimPlot(input_sce.all， reduction = "umap"， group.by = "RNA_snn_res.0.3") +    ggtitle("louvain_0.3"))  ggsave(plot=p1_dim， filename="Dimplot_diff_resolution_high.pdf"，width = 18)    p2_tree=clustree(input_sce.all@meta.data， prefix = "RNA_snn_res.")  ggsave(plot=p2_tree， filename="Tree_diff_resolution.pdf")  table(input_sce.all@active.ident)   saveRDS(input_sce.all， "sce.all_int.rds")  return(input_sce.all)  }

可视化

图片

Tree_diff_resolution

图片

Dimplot_diff_resolution_low

图片

Dimplot_diff_resolution_high结语

本期我们整理了肺腺癌单细胞数据集GSE189357个的10X样式的单细胞测序数据，在批量读取了10X文献后，进行了归拢并胜利构建Seurat对象，在此基础上将患者的临床信息添加到meta.data分组信息中。在此基础上走Seurat V5步调经由，对Seurat对象进行QC质控、并诓骗harmony整合去批次、并按步调经由进行降维聚类分群。下一期，我们将进行单细胞测序数据分析最紧迫（个东说念主觉得）的一步：细胞谛视。下一期我们不见不散！

图片

本站仅提供存储办事，扫数内容均由用户发布，如发现存害或侵权内容，请点击举报。