“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”整合不同平臺的單細胞數據之旅
  

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的单细胞测序分析之小技巧之for循环批量处理数据和出图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。