前段时间买了果子学生信的第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个小时。虽然代码很烂,但是得出结果,得出那些批量的数据与图的那一刻,是激动的,是满心欢喜的。感谢自己的坚持。