前段时间买了果子学生信的第7个课程，GZ_07，是关于批量操作的。最开始学习do.call这个函数的使用也是从果子老师的公众号推文里面学到的。

大概有两个小时的视频，一边学习，一边操作，花了半天的时间。总的来说还是蛮有收获的，总结如下：

1. 关于for循环的使用

“能够处理一次，就能够处理多次”。这是精髓。只要能够用代码实现一次的操作，比如计算相关性，生存分析，绘图等，就可以通过for循环进行多次批量的操作。，把需要实现的目标写下来，一个步骤一个步骤的实现。然后再通过for循环实现批量的操作。

比如，将HALLMARKS的条目写成我们喜欢的没有下划线的那种大小写的形式：Tnfa signaling pathway。

第一步：将下划线改为空格；

第二步：将前面的HALLMARK去掉；

第三部：改成首字母大写，其他小写的形式。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
rm(list = ls())
load(file = "data/hallmarks_gene_set.Rdata")
#1.把下划线替为空格
library(stringr)#我喜欢使用stringr来处理字符
i= 1
test= as.character(hallmarks[i,1])
test = str_replace_all(test, '_',' ')
#2.去掉Hallmark
test = str_remove(test,"HALLMARK ")
#3. 首字母大写，其他改为小写
test = str_to_sentence(test)
test

#写for循环
hallmarks$term = as.character(hallmarks$term)
for (i in 1:nrow(hallmarks)) {
  print(i)
  test = hallmarks[i,1]
  test = str_replace_all(as.character(test), '_', ' ')
  test = str_remove(test, 'HALLMARK ')
  hallmarks[i,1] = str_to_sentence(test)
}

课程里面还讲到，容器的概念。就是使用for循环，最后输出的大部分情况是数据框，这时候我们先要建立一个新的容器，把这些需要输出的内容都写进这个数据框中。

for循环讲到的第二个例子：统计这些基因在显著富集通路中的次数，也是很有实际意义的，有时候需要看某个基因富集在哪些通路，或者看某个基因在显著富集通路中出现频次，就可以用到这里提到的方法。

2. 关于函数

函数方面没有太多需要回顾的。如果一些代码重复使用超过三次，就值得把它封装成函数来使用。写函数的思路跟for循环是一样的，首先将一次的操作写下来，然后写成for循环或者写成函数。只是要知道怎么把for循环中的i编程函数的变量。

函数可以有多个输入的变量，而输出的话默认是最后一个代码的结果。如果想要输出多个结果，就需要用到return，比如return(c(a,b))。

3. lapply和future_lapply

写函数的目的是可以让我们的代码简洁，使用方便可重复。对于R语言来说另外一个好处就是可以结合apply系列的函数来使用，运行的效率通常比for循环要快。使用lapply相比for循环而言，不需要事先定义一个存放数据的容器，而是直接返回list。

首先看lapply函数，它是输入向量和函数，输出list，然后可以通过do.call将list转变为矩阵。

在这里，举了一个例子是单基因和免疫细胞浸润分数的相关性（课程里面的每一个例子都很实用，确实是这样！）。首先，需要两组数据，一个是基因在样本中的表达量，另一个是免疫细胞在样本中的免疫浸润分数的数据，两者样本是一致的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
rm(list = ls())
load(file = "data/expr_immu_BRCA.Rdata")
### 基因KLF5和免疫浸润的相关性分析
gene <- "KLF5"
y <- as.numeric(expr_data[,gene])

#1. 写一个函数，用来计算单个基因的表达量和免疫浸润分数的相关性
mycor = function(x){#x是指免疫浸润数据的某一类细胞,y是上面定义的基因的表达矩阵
  dd = cor.test(as.numeric(immu_data[,x]),y,method ="spearman",exact=FALSE)
  data.frame(cell=x,cor=dd$estimate,p.value=dd$p.value)
}
#2. 使用lapply，返回list数据
lapplylist = lapply(colnames(immu_data),mycor)
#3. 使用docall转换list
cor_data <- do.call(rbind,lapplylist)

最后，通过绘图展示相关性。

1
2
3
4
5
6
7
8
9
10
11
12
13
### 画图展示全貌
library(dplyr)
library(ggplot2)
cor_data %>% 
  filter(p.value <0.05) %>% 
  ggplot(aes(cor,forcats::fct_reorder(cell,cor)))+
  geom_segment(aes(xend=0,yend=cell))+
  geom_point(aes(col=p.value,size=abs(cor)))+
  scale_colour_gradientn(colours=c("#7fc97f","#984ea3"))+
  #scale_color_viridis_c(begin = 0.5, end = 1)+
  scale_size_continuous(range =c(2,8))+
  theme_bw()+
  ylab(NULL)

得到课程封面的这个棒棒糖图：

图1 在乳腺癌中KLF5与免疫细胞浸润分数的相关性图

future_lapply是future.apply包中的函数，该R包可以实现apply函数的并行计算，而且代码相对于其他并行计算的方法简单得多。

lapply到future_lapply其实很简单：

1
2
3
4
library(future.apply)
plan(multisession)
#然后将lapply改为future_lapply即可
lapplylist = future_lapply(colnames(immu_data),mycor)

4. 其他的一些例子

批量操作的内容基础基本上就是上面这些。然后，课程通过其他几个实用的例子来进行实践：

  • 单基因的相关性分析

  • 如何整理从TCGA下载的原始数据

    这一点之前我都不敢自己操作，看到metadata的那些数据就乱。都是通过GDCquery函数直接下载的，因为下载完后就直接整理成dataframe了。现在，我也不怕自己组装TCGA的数据了。

  • 批量生存分析

  • 两个基因在33个癌种中的相关性以及绘图

    在这里，学到了一个函数split，可以将33个癌种的数据裂解开，以方便进行批量的操作。在不知道split之前，按照我的理解，我可能会先建立一个33个癌种的字符向量，然后再使用for循环进行计算，但是这种方式就很难转换为lapply的写法，从而在速度上大大的变慢了。

    如果是单纯计算两个基因的相关性，可以这么写：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    #通过管道的写法，计算两个基因的相关性
    rm(list=ls())
    load(file = "data/pancancer_mRNA_exprSet.Rdata")
    data <- mRNA_exprSet
    Result = data %>%
          group_by(type) %>%
          summarize(
            pvalue = cor.test(as.numeric(METTL3),as.numeric(SETD2))$p.value,
            cor =  cor.test(as.numeric(METTL3),as.numeric(SETD2))$estimate
        )
    #但是这种方法可能没有办法转换为函数，而且计算了两次cor.test，速度就慢了

    图2 两个基因在不同癌种中的相关性

  • 单个基因在33个癌种中与哪些基因具有较好的相关性。

总结

R语言的学习，或者说生信分析的学习，模仿是很重要的一点。学习前辈的优秀代码，思维方式，然后用到自己的项目中，就变成自己的内容。这个课程还是很有收获的。

记得最开始学习文献的项目是我写的这篇文献分享：读文献-大规模全基因组基因芯片数据集中的布尔推断网络。由于开始学习R语言不久，找来作者的所有相关的文章，找代码，因为不懂不理解搞不定，中间放下好几次，最后批量复现成功，一次代码跑完花了将近6个小时。虽然代码很烂，但是得出结果，得出那些批量的数据与图的那一刻，是激动的，是满心欢喜的。感谢自己的坚持。