前言:昨晚听了师兄开的博客熊言熊语,印象最深的Jimmy, 以及徐洲更,果子的两期,顺带做了一枚DIY的发卡。除去Jimmy的有趣经历,印象最深的就是两期里都提到的,要以任务为导向去学习编程。我的师姐们也想学习编程,因此我想给她们整理一些小任务,帮助她们更好地入门R。我也不知道能不能坚持到整理出100个,只希望师姐们不会催我。
1. 基因ID转换(Gene ID Conversion)
首先来看看下面这串基因ID,你认识吗?
TP53
P53
BCC7
LFS1
BMFS5
TRP53
ENSG00000141510
7157
但是如果我告诉你这些都是同一个基因,你会惊讶吗?但是其实这里面只包含来三种Gene ID类型,第一个是SYMBOL(),第二个到第六个是它的别名,ENS开头的是ENSEMBL ID, 最后那个纯数字的是ENTREZ ID(目前国际上最权威的Gene ID编号)。
TP53(SYMBOL)
ENSG00000141510(ENSEMBL)
7157(ENTREZ)
不同数据库来源有不同的Gene ID,你是不是要疯了呢?但是,远不止这些。此外,还有Uniprot, HGNC等等类型(其他的用得很少,请自行了解,一般就是上面提到的三种之间互相转换)。
背景知识讲完了,下面我们来开始做任务了。
这里需要用到Y叔的clusterProfiler和相应的物种注释包(human:org.Hs.eg.db; mouse: org.Mm.eg.db)。
由于这些包都在Bioconductor上(Bioconductor是基因组数据分析相关的软件包仓库,也有的会提交到CRAN上),无法用install.packages()
安装,所以用Bioconductor官方的安装方法即可。
# 安装相应的包 if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("org.Hs.eg.db") BiocManager::install("clusterProfiler")
安装tips:
安装时发现直接运行上面的命令无法顺利安装。需要先安装Rcpp,GO.db, DO.db。故可先运行下面的命令。
install.packages("Rcpp") if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("GO.db") BiocManager::install("DO.db")
这里需要使用clusterProfiler的bitr函数转换基因ID,需要调用标准注释库org.db类型的org.Hs.eg.db(包含了各色基因ID信息)
# 加载包 library("clusterProfiler") library("org.Hs.eg.db") genelist1 <- c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4") genelist2 <- bitr(genelist1, fromType = "SYMBOL", toType = "ENTREZID", org.Hs.eg.db) genelist2
运行结果:
再看看下面这个代码的输出和上面的有什么区别?
# 加载包 library("clusterProfiler") library("org.Hs.eg.db") genelist1 <- c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4","P53") genelist2 <- bitr(genelist1, fromType = "SYMBOL", toType = "ENTREZID", org.Hs.eg.db) genelist2
为什么出现了Warning呢?这是因为后面加的P53和TP53是同一个基因,而且SYMBOL格式的是TP53,P53就不被识别了。
如有错误,敬请指正。
2. 数据分类(within)
第二个问题是这样的,有这样一个数据表(data frame)expn, 第一列是genename,第二列是表达水平值,实际数据分析中,我们要面对的可能是几千甚至几万个基因,表达水平也有很多形式,这里只是创建了一个很简单的小数据表。
expn <- data.frame(genename=c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4"),exp=c(50,200,35,300,80,90)) expn
输出:这就是我们的要用来操作的数据表了。
## genename exp ## 1 TP53 50 ## 2 YY1 200 ## 3 DNMT1 35 ## 4 SOX7 300 ## 5 CTCF 80 ## 6 KLF4 90
以50和100为阈值,划分高中低表达组,并新建一列保存结果。
expn.new <- within(expn,{ Type <- NA Type[expn$exp >= 100] <- "High" Type[expn$exp < 100 & expn$exp > 50] <- "Medium" Type[expn$exp <= 50] <- "Low" }) expn.new
输出:
## genename exp Type ## 1 TP53 50 Low ## 2 YY1 200 High ## 3 DNMT1 35 Low ## 4 SOX7 300 High ## 5 CTCF 80 Medium ## 6 KLF4 90 Medium
数据分类在实际数据分析中很常见,如画图的时候需要根据不同类别标记不同的颜色。
3.如何对多个重复的表达值求均值?(行操作:rowMeans, rowSums)
第二个问题中数据表只有一个表达值,实际中我们会遇到重复,那么如何应对有重复值的数据呢,那就是要求均值了。一般情况下,
expn <- data.frame(genename=c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4"),rp1=c(50,200,35,300,80,90),rp2=c(45,210,40,290,85,88),rp3=c(52,205,39,309,87,93)) expn
输出:使用的数据表
## genename rp1 rp2 rp3 ## 1 TP53 50 45 52 ## 2 YY1 200 210 205 ## 3 DNMT1 35 40 39 ## 4 SOX7 300 290 309 ## 5 CTCF 80 85 87 ## 6 KLF4 90 88 93
查看2-4列,即我们要求均值的数据。
expn[2:4]
输出:
## rp1 rp2 rp3 ## 1 50 45 52 ## 2 200 210 205 ## 3 35 40 39 ## 4 300 290 309 ## 5 80 85 87 ## 6 90 88 93
求均值:
expn$Mean <- rowMeans(expn[2:4]) expn
## genename rp1 rp2 rp3 Mean ## 1 TP53 50 45 52 49.00000 ## 2 YY1 200 210 205 205.00000 ## 3 DNMT1 35 40 39 38.00000 ## 4 SOX7 300 290 309 299.66667 ## 5 CTCF 80 85 87 84.00000 ## 6 KLF4 90 88 93 90.33333
表达谱不会需要求和的,但是我懒得再创建新的数据表了,将就着练习吧。
对2-4列数据求和:
expn$Sum <- rowSums(expn[2:4]) expn
## genename rp1 rp2 rp3 Mean Sum ## 1 TP53 50 45 52 49.00000 147 ## 2 YY1 200 210 205 205.00000 615 ## 3 DNMT1 35 40 39 38.00000 114 ## 4 SOX7 300 290 309 299.66667 899 ## 5 CTCF 80 85 87 84.00000 252 ## 6 KLF4 90 88 93 90.33333 271
目前提到了rowMeans, rowSums,其实还有colMeans, colSums,不是经常会用到,自己去看吧。
4. 我们来写个for loop吧
第三个问题中提到了求均值,但是实际中往往还需要求标准差(Standard Error, SD)。在R语言中,求SD的函数就是sd()
。这里和第三个问题一样,需要创建一个小的数据表。
expn <- data.frame(genename=c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4"),rp1=c(50,200,35,300,80,90),rp2=c(45,210,40,290,85,88),rp3=c(52,205,39,309,87,93)) expn
输出:
## genename rp1 rp2 rp3 ## 1 TP53 50 45 52 ## 2 YY1 200 210 205 ## 3 DNMT1 35 40 39 ## 4 SOX7 300 290 309 ## 5 CTCF 80 85 87 ## 6 KLF4 90 88 93
这里涉及到行名的问题。创建data.frame的时候命名了列名,但是没有命名行名,这里使用rownames()
进行命名,把genename作为行名(需要注意的是,用于行名的数据不允许出现重复值)。然后新加了SD列,并填充空值(NA)。
然后就是for循环对每个基因的三个重复(2-4列)的表达值求SD,并加到SD列中。
rownames(expn) <- expn$genename expn$SD <- NA for(i in rownames(expn)){ expn[i,]$SD <- sd(expn[i,2:4]) } expn
输出:如果想比较一下,可以在for loop前打印一下expn,可以看到SD列都是NA。
## genename rp1 rp2 rp3 SD ## TP53 TP53 50 45 52 3.605551 ## YY1 YY1 200 210 205 5.000000 ## DNMT1 DNMT1 35 40 39 2.645751 ## SOX7 SOX7 300 290 309 9.504385 ## CTCF CTCF 80 85 87 3.605551 ## KLF4 KLF4 90 88 93 2.516611
一个简单的for loop就是长这样了。
5. 再来看看if else吧
第五个问题是用for loop套if else结构实现第二问中的效果。
创建和第二问中相同的data.frame。
expn <- data.frame(genename=c("TP53", "YY1","DNMT1","SOX7","CTCF","KLF4"),exp=c(50,200,35,300,80,90)) expn
输出:
## genename exp ## 1 TP53 50 ## 2 YY1 200 ## 3 DNMT1 35 ## 4 SOX7 300 ## 5 CTCF 80 ## 6 KLF4 90
然后就是一样的,先写一个for loop,再来写if else。下面涉及的if…else if …else是if...else if...else
语句的基本语法结构,这只是if else中的一种形式,也是我们分析中最常需要的一种。
rownames(expn) <- expn$genename expn$Type <- NA for(i in rownames(expn)){ if(expn[i,]$exp >= 100){ expn[i,]$Type <- "High" } else if(expn[i,]$exp <= 50){ expn[i,]$Type <- "Low" } else{ expn[i,]$Type <- "Medium" } } expn
输出:结果和第二问中一致。
## genename exp Type ## TP53 TP53 50 Low ## YY1 YY1 200 High ## DNMT1 DNMT1 35 Low ## SOX7 SOX7 300 High ## CTCF CTCF 80 Medium ## KLF4 KLF4 90 Medium
如果想要阅读更舒服,可前往语雀原文链接:https://www.yuque.com/docs/share/f52fe166-d1f4-469a-8a4a-7d1a05567fc5?#