单细胞测序分析之小技巧之for循环批量处理数据和出图
"harmony"整合不同平台的单细胞数据之旅生物信息学习的正确姿势
NGS系列文章包括NGS基础、转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这) 、ChIP-seq分析 ( ChIP-seq基本分析流程) 、单细胞测序分析 ( 重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程 (原理、代码和评述)) 、DNA甲基化分析、重测序分析、GEO数据挖掘( 典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集) 等内容。
在进行单细胞转录组测序分析中,我们发现比如样本较多或者需要大量出图的时候,我一开始就是大量手动一个一个的出图,但回头想想,这样的操作模式不都是一样的嘛,直接用for循环不就搞定啦!基础
首先我们讲点for循环的基础知识及举个小栗子!
for循环基本结构如下:for(变量 in 值){}
也就是说当变量在值的范围内将执行中括号内的操作。是不是非常简单?
我们举个栗子:
比如我们想要计算一个向量中偶数的个数:x <- c(2,5,3,9,8,11,6) count <- 0 for (val in x) { if(val %% 2 == 0) count = count+1 } print(count)[1] 3
在上面的示例中,由于向量x具有7个元素,因此循环迭代了7次。
在每次迭代中,val取x的对应元素的值。
我们使用了一个计数器来计算x中的偶数。我们可以看到x包含3个偶数。进阶
比如我们现在有两个患者的鼻腔样本,然后我们进行单细胞测序后,cellranger后我们在filtered中分别生成了3个文件:barcodes,features和matrix。我懒呀,我想万一我有好多个样本怎么办,不如用一个for循环来搞定!
于是我的文件就成了这个样子: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,我应该怎么办呢?嗨,循环起来呀!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 ,再用patchwork 或cowplot 进行拼图。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)这个Nature推荐的代码海洋竟然有文章作者上传的所有可重现性脚本,涉及单细胞、微生物组、转录组分析、机器学习等相关 10X单细胞测序分析软件:Cell ranger,从拆库到定量 NBT|45种单细胞轨迹推断方法比较,110个实际数据集和229个合成数据集 重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程 (原理、代码和评述) "harmony"整合不同平台的单细胞数据之旅