search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
122 stars 30 forks source link

十分钟完成上周的Science杂志的玉米单细胞文章的数据处理 #3911

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/YvhXhEV4DXRRU872XmI4Ag

ixxmu commented 1 year ago

十分钟完成上周的Science杂志的玉米单细胞文章的数据处理 by 生信技能树


我们有一个学徒培养环节，主打《单细胞数据处理交流》，全部的的数据分析，目前已经进行到70%，文末有参加交流会的腾讯会议直播方式！

下面12月份学徒的投稿

但是仍然是有部分学徒迟迟不肯开始，觉得单细胞过于复杂，为了最大化激励大家行动起来，我们甚至承诺大家只需放开手分享，提供带有数据集的文献，我们直接帮忙写代码，让学徒们去理解和解读。手把手教学

然后有一个小伙伴发在交流群的pdf文献，就是上周（2021年12月2日）的文章，Science在线发表的美国纽约大学Kenneth D. Birnbaum团队题为《Ground tissue circuitry regulates organ complexity in maize and Setaria》的研究论文，使用的是是Seurat的标准降维聚类分群流程，所以得到了编号为 0-20这21个单细胞亚群，如下所示：

编号为 0-20这21个单细胞亚群

可以看到，在Seurat的标准降维聚类分群流程后，研究者们就直接给了这些细胞亚群一个生物学名字，一般来说亚群命名都是需要特定的标记基因，来源于过去一百多年的生命科学领域蓬勃发展的文章。但是我看了看文章里面的：

特定的标记基因

感觉基因名字都稀奇古怪的，简单的一个编号而已，也不知道是哪个团队在规范和维护这样的基因名字。

下面开始10分钟复现这个Science杂志的玉米单细胞文章的数据处理

其单细胞数据集在 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173087 ，可以看到是9个样品：

GSM5259613 Zmroot_sc_read1/2_rep1
GSM5259614 Zmroot_sc_read1/2_rep2
GSM5259615 Zmroot_sc_read1/2_rep3
GSM5603496 Zmroot_sc_read1/2_rep4
GSM5603497 Zmroot_sc_read1/2_rep5
GSM5603498 Zmroot_sc_read1/2_rep6
GSM5603499 Zmroot_sc_read1/2_rep7
GSM5603500 Zmroot_sc_read1/2_rep8
GSM5603501 Zmroot_sc_read1/2_rep9

其中下面的 17.4 Mb 文件存储了不到5000个单细胞：

GSE173087_Maize_rootsc_counts.txt.gz 

大家需要简单下载它，然后复制粘贴我的30行左右的代码就可以运行啦：

library(Seurat)

library(data.table) 
# 读取   17.4 Mb 文件 文件 
ct=fread( 'GSE173087_Maize_rootsc_counts.txt.gz',data.table = F)
ct[1:4,1:4]
rownames(ct)=ct[,1]
ct=ct[,-1] 
## Initialize the Seurat object with the raw (non-normalized data).
maize <- CreateSeuratObject(counts = ct, 
                           project = "maize", 
                           min.cells = 3, 
                           min.features = 200)
maize 
## Normalizing the data
maize <- NormalizeData(maize, normalization.method = "LogNormalize", 
                      scale.factor = 10000)

maize <- NormalizeData(maize)

## Identify the 2000 most highly variable genes
maize <- FindVariableFeatures(maize, selection.method = "vst", nfeatures = 2000)

## In addition we scale the data
all.genes <- rownames(maize)
maize <- ScaleData(maize, features = all.genes)

maize <- RunPCA(maize, features = VariableFeatures(object = maize), 
               verbose = FALSE)
maize <- FindNeighbors(maize, dims = 1:10, verbose = FALSE)
maize <- FindClusters(maize, resolution = 0.5, verbose = FALSE)
maize <- RunUMAP(maize, dims = 1:10, umap.method = "uwot", metric = "cosine")
table(maize$seurat_clusters)
phe=maize@meta.data
save(phe,file = 'phe-by-basic-seurat.Rdata')

# maize.markers <- FindAllMarkers(maize, only.pos = TRUE, min.pct = 0.25,  logfc.threshold = 0.25, verbose = FALSE)
DimPlot(maize, reduction = "umap", group.by = 'seurat_clusters',
        label = TRUE, pt.size = 0.5) 

我上面的代码其实就是Seurat的标准降维聚类分群流程，可以看到是编号为0-18的单细胞亚群：

编号为0-18的单细胞亚群

接下来需要对它们进行已知的标记基因的可视化，代码也是超级的简单：


head(rownames(maize))
cg= c( 'ZmMYB46','ZmCO2','ZmLBD29','ZmLCR','ZmRHD6')
cg=c('d012714','d013033','d052952','d002191','d023651',
     'd035689','d051478','d017508','d046778','d043242',
     'd048131','d019985','d042035','d041611','d049606',
     'd025402','d011156','d008925','d022457','d046186','d026406')
library(ggplot2)
p1=DotPlot(maize, features =paste0('Zm00001',cg) ,
        group.by = 'seurat_clusters')  + coord_flip()
p2=DimPlot(maize, reduction = "umap", 
        label = TRUE, pt.size = 0.5) + NoLegend()
ggsave(p1,filename = 'markers.pdf')
ggsave(p2,filename = 'umap.pdf')
library(patchwork)
p1+p2
ggsave( filename = 'markers-and-umap.pdf',width = 15)

if(F){
  ## Assigning cell type identity to clusters
  new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T",
                       "FCGR3A+ Mono", "NK", "DC", "Platelet")
  names(new.cluster.ids) <- levels(maize)
  maize <- RenameIdents(maize, new.cluster.ids)
} 

出图 如下所示：

奇奇怪怪的基因名字

可以看到这些基因名字都是加上了一个前缀  Zm00001，我不知道为什么这个玉米物种的数据库的特殊性版本啊！

看到了这里如果你10分钟无法完成我们这个教程的两个核心图表

说明你不懂R语言，不知道如何去安装我教程里面的提到的R包，不知道如何复制粘贴运行代码。

当然了，出图并不意味着结束，我们确实是10分钟出来了这个Science杂志的玉米单细胞文章的数据处理的重要图表，可是我们对这个Science文献的理解完全没有任何增长。

期待大家加入我们的单细胞交流群啊，只需要你找到感兴趣的单细胞文章，而且愿意分享它，我们就给你写代码分析它！

你分享，我分析，强强联合！

每晚九点，持续100次！欢迎对单细胞感兴趣的小伙伴加入一起讨论学习以及 提高！（入会密码是 1024 ）

虽然说上面的代码都是复制粘贴即可运行，但是如果要更好地完成上面的图表，通常是需要掌握5个R包，分别是: scater,monocle,Seurat,scran,M3Drop，需要熟练掌握它们的对象，：一些单细胞转录组R包的对象 而且分析流程也大同小异：

  • step1: 创建对象
  • step2: 质量控制
  • step3: 表达量的标准化和归一化
  • step4: 去除干扰因素(多个样本整合)
  • step5: 判断重要的基因
  • step6: 多种降维算法
  • step7: 可视化降维结果
  • step8: 多种聚类算法
  • step9: 聚类后找每个细胞亚群的标志基因
  • step10: 继续分类

单细胞转录组数据分析的标准降维聚类分群，并且进行生物学注释后的结果。可以参考前面的例子：人人都能学会的单细胞聚类分群注释 ，我们演示了第一层次的分群。

如果你对单细胞数据分析还没有基础认知，可以看基础10讲：

写在文末

咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！见：急！计划招募100个单细胞爱好者，免费学全套单细胞降维聚类分群和生物学亚群注释