ixxmu
commented
1 year ago 给你一份批量数据处理的学习笔记，感受一下编程的魅力！ by 果子学生信
果子荐读
编程带给我最大的震撼应该就是批量处理。他让一些数据分析从不能到能，从慢到快。这是我获益最多的技能。
有了这个技能，我才感到自己无可替代。因为批量是个通用技能，但是批量的那个动作是个性化需求。
而个性化的需求最终只能自己实现。
数据分析，数据挖掘没有尽头，第一，数据源源不断，第二，认知持续提升。在这其中，拥有一个傍身的技能是行走江湖的基本要求。
因此，我总结了常用的批量流程，写成了教程，也录制了视频。
视频课程: R语言中的批量操作，裂变你的技能。
教程发布后，很多朋友迅速学习，立马用到自己的课题中去。
我邀请了两位朋友，给大家分享一下学习感受。
听讲座，看分享，并不能立竿见影地提高技能，但他能够拓宽你的眼界，告诉你努力的底线和上限。
今天是第一位，是Cloudy，2年的老朋友了。
使用幕布梳理思路，使用markdown写作，使用图床保存图片，这个习惯很有辨识度啊。
熊确实影响了不少人。
以下是正文
最近学了果子老师推出的GZ07课程，略有所获，下面分享下自己的学习心得。
一、课程框架与主要内容
这个课程总共提供了10个脚本，学习视频时长约150min。课程的主体是关于R语言中批量操作的一些代码编写技巧，然后课程提供了丰富的例子对如何使用循环和函数批量操作数据进行了实践。
其实如果单从批量操作学习的角度，个人感觉到脚本07生存分析那里课程其实就挺完整的了，基本上多过几遍能学习到for循环到基本函数编写的基本技能。但是随着课程学习的循序渐进，我发现从脚本08、脚本09、脚本10及附加脚本pancorplot才是课程的精华所在。
这三个脚本使用XENA数据库的泛癌数据进行推演运算，里面涵盖的知识点包括且不限于：批量相关性分析、统计结果的提取、条件判断及ggplot2、for循环嵌套函数、函数嵌套函数、R语言控制结构、ggplot2的图形语法等。基本是课程前面课程讲述的所有知识点的一次总结和综合应用。
二、课程中令我印象深刻的函数们
数据库最后4个脚本里面几个重点的函数是十分值得学习的：
函数1：pancor
功能： 计算任意两个基因在泛癌数据中的相关性
pancor <- function(gene1,gene2,data=splitdata){
do.call(rbind,lapply(data, function(x){
dd=cor.test(as.numeric(x[,gene1]),as.numeric(x[,gene2]),type="pearson")
data.frame(type=x$type[1],cor=dd$estimate,p.value=dd$p.value )
}))}
函数2：pancorBatch
功能： 批量计算两个基因在泛癌数据中的相关性并按相关性系数进行筛选相关的癌种
pancorBatch <- function(gene1,gene2,splitdata){
## 一对基因批量分析
pancordata = do.call(rbind,lapply(splitdata, function(x){
dd <- cor.test(as.numeric(x[,gene1]),as.numeric(x[,gene2]),type="pearson")
data.frame(type=x$type[1],cor=dd$estimate,p.value=dd$p.value )
}))
## 提取数据
subdata = pancordata[pancordata$p.value < 0.05 ,]
## 汇总数据
data.frame(gene1=gene1,
count= nrow(subdata),
avercor= mean(subdata$cor),
positive= sum(subdata$cor > 0),
engitive= sum(subdata$cor < 0))
}
这个函数初看还是有点高能的，因为里面是function+function的结构，但其实认真读进去会发现其实很好理解，splitdata是已知的泛癌数据使用split函数切割后的数据，下位的函数传递的每一个x都是泛癌数据的表达矩阵。
所以，其实只有两个参数是未知的，两个基因，gene1和gene2，只要设定好这两个实际的参数，整个函数就可以启动了。也因为泛癌数据涉及33种癌症，所以基本上，无论gene1和gene2输入了啥，出来的pancordata都是一个33行，3列的矩阵。大概长这样：
后面，subdata是兴趣基因对在所有癌种中的相关性筛选出来的p值小于0.05的结果。实际上还是两个基因在泛癌中的相关性。因为输出结果是gene1，所以其实函数默认了我们的兴趣基因是gene2。举个例子，我设定gene1为ESR1，那么单次运算出来的结果就如下图：
count是癌种的数目，avecor是平均的相关性系数，positive和engitive是正相关和负相关的个数。
注意这里是用相关性系数＞0或者＜0去筛选的，自己做的时候可以进一步筛选
函数3：pancorplot
pancorplot <- function(data=plotdf,anotate="none"){
require(ggplot2)
require(ggrepel)
options(scipen = 2)
p=ggplot(plotdf,aes(-log10(p.value),cor))+
geom_point(data=subset(plotdf, plotdf$p.value >= 0.05),size=6,fill="grey",alpha=0.6,shape = 21,colour="black",stroke = 1.5)+
geom_point(data=subset(plotdf, plotdf$p.value < 0.05 & plotdf$cor>=0),size=6,fill="red",alpha=0.6,shape = 21,colour="black",stroke = 1.5)+
geom_point(data=subset(plotdf, plotdf$p.value < 0.05 & plotdf$cor< 0),size=6,fill="blue",alpha=0.6,shape = 21,colour="black",stroke = 1.5)+
scale_y_continuous(expand = c(0,0),limits = c(-1.1,1.1),breaks = seq(-1,1,0.2))+
scale_x_log10(limits = c(0.01, 1000),breaks = c(0.01,0.1,10,1000))+
geom_hline(yintercept = 0,size=1.5)+
geom_vline(xintercept = -log10(0.05),size=1.5)+
labs(x=bquote(-log[10]~italic("P")),y="Pearson correlation (r)")+
theme(axis.title=element_text(size=20),
axis.text = element_text(face = "bold",size = 16),
axis.ticks.length=unit(.4, "cm"),
axis.ticks = element_line(colour = "black", size = 1),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1.5),
plot.margin = margin(1, 1, 1, 1, "cm"))
index = anotate == "none"
index = index[1]
if(index){
p
}else{
p + geom_point(data=subset(plotdf, plotdf$type%in% anotate),
size=6,shape = 21,colour="green",stroke = 2)+
geom_label_repel(data=subset(plotdf, plotdf$type %in% anotate),
aes(label=type),col="black",size=6,
box.padding = 10,
arrow=arrow(angle = 30, length = unit(0.25, "inches"),
ends = "first", type = "closed"),
segment.size=1,
segment.color = "red")
}
}
出的图大概长这样：
这个函数本质上是点图，前面写的思路和ggplot2绘制火山图是如出一辙的（嗯，火山图本质上也是点图），就是不断在提取数据的子集，然后在画布上打点、上色，最后加上线条。
但是它的构建逻辑很好，后面那个标签部分是值得摸索和学习的。画出来的图打的标签都能令人基本满意，如果不满意，解决方法果子老师也提供了，在GZ01的更新课程中，使用export包导出成ppt自己改到满意为止。("▔□▔)/ 这个包的链接在这：
当然，esquisse包也是个不错的选择，看个人喜好和习惯。
之后课程中脚本09和脚本10继续整活，前面我们已知:
①pancor函数：能让我们计算兴趣基因和任意一个其他基因在泛癌数据中的相关性
②pancorBatch函数：能够让我们获取任意兴趣基因和任意一个其他基因在33种癌症中的相关性情况。
③pancorplot函数：能够可视化pancorBatch 函数的结果进行可视化。
那么如果我们想知道兴趣基因和所有其他基因在泛癌数据中的整体相关性情况，怎么整？对，启动批量，lapply或者future.apply往上面一套就出来了。如果你想让lapply也带个进度条，那pblapply包了解下。
图一个一个画不过瘾，for循环+ggsave函数安排上，噌噌噌，你的文件夹里都是图。
通常情况下，“套娃行为”会让人困惑，但是编程里面这样的套娃却让人感到十足的快乐，多少有一种不经意间就解放双手，做到超级多的事情满足感！
如果觉得这样的需求很需要，但是代码还是看的云里雾里，那么没有问题，课程里面果子老师授人以渔也授人以鱼了。直接套用脚本08/09后面现成脚本，把他们粘贴复制，创建个新脚本，导入自己的数据，以上的需求基本上也就分分钟解决了。什么叫用1%的时间，解决99%的问题呀（战术后仰）╮(╯▽╰)╭。
诚如果子老师最后总结，GZ07是可以和可以和GZ01课程对接的，先使用同样的方法，筛选在GETx数据库中和我们手头兴趣基因表达呈正相关的基因（或者负相关）。然后再泛癌数据中再操作一遍，筛选和我们兴趣基因表达呈负相关（或者正相关）的基因。这样的基因对可能因为癌症的发生表达调控发生改变，可以作为以后研究的idea。
三、学完之后自己也整个函数玩一下
GZ07可以当成需求的实现，也可以作为对R语言编程的一次系统入门或者复习回顾，学完之后我也尝试了整了个函数玩玩，代码如下：
library(tidyverse)
##生成示例数据
set.seed(3)
data.frame(gene=letters[1:10],value1=1:10,
value2=1:10,value3=1:10,
value4=1:10,value5=1:10,
value6=1:10) %>%
column_to_rownames("gene") %>%
t() %>%
as.data.frame() %>%
mutate(group=rep(LETTERS[1:2],each=3)) %>%
mutate(kt=rnorm(6),bt=sample(seq(1,6,0.9),6,T)) %>%
select(group,everything())->dat
###构建函数
Compare_means<- function(gene, dat) {
# 方差齐性检验
var_check <- var.test(dat[, gene] ~ group, data = dat)
# 如果方差是0或者无法计算为缺失值时，数据会无法执行正态分布检验，所以直接设置返回缺失值
if (var_check$estimate == 0 | is.na(var_check$estimate)) {
return(data.frame(gene = gene, statistic = NA, pvalue = NA, method = NA))
}
# 正态性检验
shap <- tapply(dat[, gene], dat$group, shapiro.test)
# 如果两组有任意一组数据符合正态分布且两组数据对比方差齐则执行t检验，否则执行 wilcox.test
#shapiro.test：p值大于0.05符合正态分布
#var.test：p值大于0.05符合方差齐性
# 判断执行t检验还是秩和检验并返回统计量、p值
if (all(c(shap[[1]]$p.value, shap[[2]]$p.value) > 0.05 && var_check$p.value > 0.05)) {
test <- t.test(dat[, gene] ~ group, data = dat)
res <- data.frame(gene = gene, statistic = test$statistic, pvalue = test$p.value, method = "t.test")
} else {
test <- wilcox.test(dat[, gene] ~ group, data = dat, exact = F)
res <- data.frame(gene = gene, statistic = test$statistic, pvalue = test$p.value, method = "wilcox.test")
}
return(res)
}
do.call(rbind,lapply(names(dat)[-1], Compare_means,dat))
示例数据输出结果
gene statistic pvalue method
1 a NA NA <NA>
2 b NA NA <NA>
3 c NA NA <NA>
4 d NA NA <NA>
5 e NA NA <NA>
6 f NA NA <NA>
7 g NA NA <NA>
8 h NA NA <NA>
9 i NA NA <NA>
10 j NA NA <NA>
t kt -0.04193745 0.9686219 t.test
W bt 4.00000000 1.0000000 wilcox.test
我用这个小函数实现了在自己的数据中批量对基因进行t检验或者秩和检验，想看看这样分析基因在癌旁和肿瘤中的表达和差异分析结果会有什么不同，有没有哪些可以注意的地方。函数加上future.apply，5分钟不到在我的小电脑上2w+个基因就算完了，一和差异基因对比，还真是有一些不同的地方，也带来了一点新的启发。
当然在写的过程中其实也遭遇了一些列的报错和调试，其实也是挺麻烦的。报错的点出现在数据的第5000多个基因那里，因为表达量问题，无法执行正态检验。还有因为数据方差是0报错的，总之一个小小的函数也踩了不少坑。我后来是通过先将函数改写成for循环，让循环打印迭代次数才捕捉到这些个报错然后加判断语句绕过它们的。
哈哈，也期待果子老师在后续课程更新可以分享一些关于函数书写、调试、报错debug的经验和知识~("▔□▔)/
function的报错最终不会返回任何结果，只有控制台的ERROR。对比for循环，如果我们创建了容器，那么至少报错之前的结果也是可以存下来的。这大概也是两者的一个不同的点吧。
批量会为数据分析带来无限可能，for循环是基础，lapply的隐式循环+向量化编程会让操作提速提速再提速~
分享到这就告一段落吧，祝各位大佬眼明心亮，在GZ系列课程中都学有所获，早日走上数据分析的无忧之路~( ^_^ )/~~拜拜
https://mp.weixin.qq.com/s/OXiKZrpFyEchzsfcaP95jQ