前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >单细胞测序分析之小技巧之for循环批量处理数据和出图

单细胞测序分析之小技巧之for循环批量处理数据和出图

作者头像
生信宝典
发布2022-01-18 20:18:41
6650
发布2022-01-18 20:18:41
举报
文章被收录于专栏:生信宝典

在进行单细胞转录组测序分析中,我们发现比如样本较多或者需要大量出图的时候,我一开始就是大量手动一个一个的出图,但回头想想,这样的操作模式不都是一样的嘛,直接用for循环不就搞定啦!

基础

首先我们讲点for循环的基础知识及举个小栗子!

for循环基本结构如下:

代码语言:javascript
复制
for(变量 in 值){}

也就是说当变量在值的范围内将执行中括号内的操作。是不是非常简单?

我们举个栗子:

比如我们想要计算一个向量中偶数的个数:

代码语言:javascript
复制
x <- c(2,5,3,9,8,11,6)
count <- 0
for (val in x) {
if(val %% 2 == 0)  count = count+1
}
print(count)
代码语言:javascript
复制
[1] 3

在上面的示例中,由于向量x具有7个元素,因此循环迭代了7次。

在每次迭代中,val取x的对应元素的值。

我们使用了一个计数器来计算x中的偶数。我们可以看到x包含3个偶数。

进阶

比如我们现在有两个患者的鼻腔样本,然后我们进行单细胞测序后,cellranger后我们在filtered中分别生成了3个文件:barcodes,features和matrix。我懒呀,我想万一我有好多个样本怎么办,不如用一个for循环来搞定!

于是我的文件就成了这个样子:

代码语言:javascript
复制
batch_list=list("P2","P3")
batch_data_list=list("P2"=1,"P3"=1)
for( i in 1:length(batch_list))
{
  print(batch_list[[i]])
  s_object=Read10X(paste("~/Input_files/",batch_list[[i]],sep=""))
  s_object=CreateSeuratObject(counts =s_object, min.cells = 0, min.features = 400, project = "P23")
  s_object[["percent.mt"]] <- PercentageFeatureSet(s_object, pattern = "^MT-")
  s_object <- subset(s_object, subset = nFeature_RNA >100 & nFeature_RNA <8000 & percent.mt <10)
  s_object@meta.data[, "run"] <- batch_list[i]
  s_object=NormalizeData(s_object)
  batch_data_list[[i]]=FindVariableFeatures(s_object, selection.method = "vst", nfeatures =5000)
}

那么我们仔细看一下刚才发生了什么,我们首先把我们的“P2”和“P3”设置为list,然后在for循环中我们分别进行了读取数据,提取线粒体基因比例,QC筛选,在metadata中添加新的一列,进行归一化并计算高变基因。最后将P2和P3合并在一个list中。

这时候一定会有好同志问这样一个问题,为什么在batch_data_list=list(“P2”=1,”P3”=1)中将P2和P3都赋值为1,这时候我们不妨不对其进行设置,使batch_data_list=list(“P2”,”P3”),我们会看见下图中的P2会消失哦!

在我们使用seurat中的FindAllMarkers()得到每个cluster的高变基因后,我也同时得到了一个csv表,可是我觉得太不直观了,于是我现在要循环出一些不同clusters的vlnplot,我应该怎么办呢?嗨,循环起来呀!

代码语言:javascript
复制
clustersss <-
  list(
    "0",
    "1",
    "2",
    "3",
    "4",
    "5",
    "6",
    "7",
    "8",
    "9",
    "10",
    "11",
    "12",
    "13",
    "14",
    "15",
    "16",
    "17",
    "18",
    "19",
    "20",
    "21",
    "22"
  )

for (i in clustersss) {
  for (m in 1:nrow(run.combined.markers)) {
    if (run.combined.markers["cluster"][m, ] == i) {
      filename <- paste(run.combined.markers$gene[m], 'VlnPlot.pdf', sep = '_')

      p <-
        VlnPlot(object = run.combined,
                features = c(run.combined.markers$gene[m]))

      print(p)
      ggsave(p,
             filename = paste(i, run.combined.markers$gene[m], 'VlnPlot.pdf', sep = '_'))
      dev.off()
    }
  }
}

我给解释一下上面的内容,首先我们把我们的cluster设为list,i代表cluster,m代表run.combined.marker的排序,使用两个for循环进行嵌套,最后在保存文件时将cluster+基因名+vlnplot结合进行保存。

每次看见这样出图我都特别有成就感,,,,哈哈哈哈,快have a try!

其实也可以写一个apply版的,获得所有plotList,再用patchworkcowplot进行拼图。

代码语言:javascript
复制
plotMarker <- function(cluster, run.combined) {
  for (m in 1:nrow(run.combined.markers)) {
    if (run.combined.markers["cluster"][m,] == cluster) {
      filename <-
        paste(run.combined.markers$gene[m], 'VlnPlot.pdf', sep = '_')

      p <-
        VlnPlot(object = run.combined,
                features = c(run.combined.markers$gene[m]))

      ggsave(p,
             filename = paste(i, run.combined.markers$gene[m], 'VlnPlot.pdf', sep = '_'))

      return(p)
    }
  }
}

plotList = lapply(clustersss, plotMarker, run.combined = run.combined)
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-06-21,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信宝典 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 基础
  • 进阶
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档