GZ07课程学习笔记
前段时间买了果子学生信的第7个课程,GZ_07,是关于批量操作的。最开始学习do.call这个函数的使用也是从果子老师的公众号推文里面学到的。
大概有两个小时的视频,一边学习,一边操作,花了半天的时间。总的来说还是蛮有收获的,总结如下:
1. 关于for循环的使用
“能够处理一次,就能够处理多次”。这是精髓。只要能够用代码实现一次的操作,比如计算相关性,生存分析,绘图等,就可以通过for循环进行多次批量的操作。,把需要实现的目标写下来,一个步骤一个步骤的实现。然后再通过for循环实现批量的操作。
比如,将HALLMARKS的条目写成我们喜欢的没有下划线的那种大小写的形式:Tnfa signaling pathway。
第一步:将下划线改为空格;
第二步:将前面的HALLMARK去掉;
第三部:改成首字母大写,其他小写的形式。
代码:
1 | rm(list = ls()) |
课程里面还讲到,容器的概念。就是使用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 | rm(list = ls()) |
最后,通过绘图展示相关性。
1 | ### 画图展示全貌 |
得到课程封面的这个棒棒糖图:
future_lapply是future.apply包中的函数,该R包可以实现apply函数的并行计算,而且代码相对于其他并行计算的方法简单得多。
lapply到future_lapply其实很简单:
1 | library(future.apply) |
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,速度就慢了单个基因在33个癌种中与哪些基因具有较好的相关性。
总结
R语言的学习,或者说生信分析的学习,模仿是很重要的一点。学习前辈的优秀代码,思维方式,然后用到自己的项目中,就变成自己的内容。这个课程还是很有收获的。
记得最开始学习文献的项目是我写的这篇文献分享:读文献-大规模全基因组基因芯片数据集中的布尔推断网络。由于开始学习R语言不久,找来作者的所有相关的文章,找代码,因为不懂不理解搞不定,中间放下好几次,最后批量复现成功,一次代码跑完花了将近6个小时。虽然代码很烂,但是得出结果,得出那些批量的数据与图的那一刻,是激动的,是满心欢喜的。感谢自己的坚持。