• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)单细胞数据分析神器-Seurat

      (未测试)单细胞数据分析神器-Seurat

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2020年2月28日
      • 评论 0评论

      近年来,单细胞技术日益火热,并且有着愈演愈烈的趋势。在2015年至2017年,甚至对某细胞群体或组织进行单细胞测序,解析其细胞成分就能发一篇CNS级别的文章。近两三年,单细胞技术从最开始的基因组,转录组测序,发展成现在的单细胞DNA甲基化,单细胞ATAC-seq等等。测序手段也从早期的10X Genomics、 Drop-seq等,发展为现在的多种多样个性化的方法。研究内容更不仅仅局限于解析细胞群体的成分,而是向研究细胞功能和生物学特性发展。今天小编向大家简单一个实用并且易上手的单细胞数据分析软件——Seurat,大家躺在床上为国家做贡献的同时也能get新技能。

      介绍一下今天的主角,Seurat是由New York Genome Center, Satija Lab开发的单细胞数据分析集成软件包。其功能不仅包含基本的数据分析流程,如质控,细胞筛选,细胞类型鉴定,特征基因选择,差异表达分析,数据可视化等等。同时也包括一些高级功能,如时序单细胞数据分析,不同组学单细胞数据整合分析等。今天,小编以官网中提供的单细胞基因表达数据为例,为大家简单介绍一下Seurat软件包中的基础分析流程,希望能够抛砖引玉,祝大家在科研的道路上越走越远。

      第一步,数据集导入

      在本教程中,我们将分析从10X基因组学免费获得的外周血单个核细胞(PBMC)数据集,来源于Illumina NextSeq 500测得的2700个单细胞转录组数据。首先,我们需要把数据集存储成Seurat可识别的数据格式,

      读入的数据可以是一个矩阵,行代表基因,列代表细胞。

      library(dplyr)
      library(Seurat)
      # Load the PBMC dataset
      pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
      # Initialize the Seurat object with the raw (non-normalized data).
      pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
      pbmc
      ## An object of class Seurat 
      ## 13714 features across 2700 samples within 1 assay
      ## Active assay: RNA (13714 features)

      数据导入成功以后,我们可以看到pbmc对象中包含了一个13714(基因数)X  2700(细胞数)的矩阵,其实在数据导入的时候,数据集中测到的少于200个基因的细胞(min.features = 200),和少于3个细胞覆盖的基因(min.cells = 3),就已经被过滤掉了。

      第二步,数据质控

      质控的参数主要有两个:1.每个细胞测到的unique feature数目(unique feature代表一个细胞检测到的基因的数目,可以根据数据的质量进行调整)2.每个细胞检测到的线粒体基因的比例,理论上线粒体基因组与核基因组相比,只占很小一部分。所以线粒体基因表达比例过高的细胞会被过滤。

      在质控前,我们需要目测一下数据集的基本情况。

      # The [[ operator can add columns to object metadata. This is a great place to stash QC stats
      pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
      # Visualize QC metrics as a violin plot
      VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

      代码运行后会得到下图所示结果,nFeature_RNA代表每个细胞测到的基因数目,nCount代表每个细胞测到所有基因的表达量之和,percent.mt代表测到的线粒体基因的比例。

      这里,我们过滤掉上图所示的离群点,并对基因的表达量进行log转换。

      第三步,鉴定细胞间表达量高变的基因(feature selection)

      这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因,用于后续鉴定细胞类型,我们使用默认参数,即“vst”方法选取2000个高变基因。

      pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
      # Identify the 10 most highly variable genes
      top10 <- head(VariableFeatures(pbmc), 10)
      # plot variable features with and without labels
      plot1 <- VariableFeaturePlot(pbmc)
      plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
      CombinePlots(plots = list(plot1, plot2))

      代码运行后会得到下图所示的结果。

      第四步,细胞分类

      1) 分类前首先要对数据集进行降维

      #Scaling the data
      all.genes <- rownames(pbmc)
      pbmc <- ScaleData(pbmc, features = all.genes)
      #Perform linear dimensional reduction
      pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
      # Examine and visualize PCA results a few different ways
      print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
      ## PC_ 1 
      ## Positive:  CST3, TYROBP, LST1, AIF1, FTL 
      ## Negative:  MALAT1, LTB, IL32, IL7R, CD2 
      ## PC_ 2 
      ## Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 
      ## Negative:  NKG7, PRF1, CST7, GZMB, GZMA 
      ## PC_ 3 
      ## Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
      ## Negative:  PPBP, PF4, SDPR, SPARC, GNG11 
      ## PC_ 4 
      ## Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 
      ## Negative:  VIM, IL7R, S100A6, IL32, S100A8 
      ## PC_ 5 
      ## Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY 
      ## Negative:  LTB, IL7R, CKB, VIM, MS4A7

      官网中提供了不同的方法可视化降维的结果,此处不再赘述。感兴趣的读者可以自行探索降维的结果,代码如下。 

      VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
      DimPlot(pbmc, reduction = "pca")
      DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)

      2) 定义数据集的“维度”

      这个小步骤就比较关键了,我们需要选择出主成分的数目,用于后续细胞分类。值得注意的是,这里定义的“维度”并不代表细胞类型的数目,而是对细胞分类时需要用到的一个参数。

      # NOTE: This process can take a long time for big datasets, comment out for expediency. More
      # approximate techniques such as those implemented in ElbowPlot() can be used to reduce
      # computation time
      pbmc <- JackStraw(pbmc, num.replicate = 100)
      pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
      JackStrawPlot(pbmc, dims = 1:15)
      ElbowPlot(pbmc)

      JackStraw(左图)和Elbow(右图)都可以决定数据的“维度”。但是Elbow比较直观,我们选择Elbow结果进行解读。可以看到,主成分(PC)7到10之间,数据的标准差基本不在下降。所以我们需要在7到10之间进行选择,为了尊重官网的建议,我们选取10,即前10个主成分用于细胞的分类。

      3) 细胞分类

      pbmc <- FindNeighbors(pbmc, dims = 1:10)
      pbmc <- FindClusters(pbmc, resolution = 0.5)
      
      ## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
      ## 
      ## Number of nodes: 2638
      ## Number of edges: 96033
      ## 
      ## Running Louvain algorithm...
      ## Maximum modularity in 10 random starts: 0.8720
      ## Number of communities: 9
      ## Elapsed time: 0 seconds

      这里我们设置了dims = 1:10 即选取前10个主成分来分类细胞。分类的结果如下,可以看到,细胞被分为9个类别。

      # Look at cluster IDs of the first 5 cells
      head(Idents(pbmc), 5)
      ## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG
      ## 1 3 1 2 6
      ## Levels: 0 1 2 3 4 5 6 7 8

      4) 可视化分类结果

      TSNE和UMAP两种方法经常被用于可视化细胞类别。

      # If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
      pbmc <- RunUMAP(pbmc, dims = 1:10)
      
      # note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters
      DimPlot(pbmc, reduction = "umap")
      
      saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

      第五步,提取各个细胞类型的marker gene

      利用 FindMarkers 命令,可以找到找到各个细胞类型中于其他类别的差异表达基因,作为该细胞类型的生物学标记基因。其中ident.1参数设置待分析的细胞类别,min.pct表示该基因表达数目占该类细胞总数的比例

      # find all markers of cluster 1
      cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
      head(cluster1.markers, n = 5)

      利用 DoHeatmap 命令可以可视化marker基因的表达

      top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
      DoHeatmap(pbmc, features = top10$gene) + NoLegend()

      第六步,探索感兴趣的基因

       Seurat提供了小提琴图和散点图两种方法,使我们能够方便的探索感兴趣的基因在各个细胞类型中的表达情况

      VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

      我们能够看到,MS4A1和CD79A两个基因在细胞群体3中特异性表达。

      FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

      这种展示方法把基因表达量映射到UMAP结果中,同样可以直观的看到基因表达的特异性。

      第七步,利用先验知识定义细胞类型

      通过对比我们鉴定的marker gene与已发表的细胞类型特意的基因表达marker,可以定义我们划分出来的细胞类群。

       最后,给我们定义好的细胞类群加上名称

      new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
          "NK", "DC", "Platelet")
      names(new.cluster.ids) <- levels(pbmc)
      pbmc <- RenameIdents(pbmc, new.cluster.ids)
      DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

      以上就是小编学习到的基本的单细胞转录组数据分析流程,整理出来供大家参考,本文主要参考了Seurat官网给出的单细胞转录组数据分析教程,若文中有小编解读错误之处,还请广大读者批评指正。

      Seurat官网地址

      https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      (未测试)R语言数据如何高效输出到文件中?
      2020年2月28日

      下一篇文章

      (未测试)再也不用费力调整GO富集分析的输出图片了
      2020年2月28日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?