您的位置 首页 生物信息学

R|解决100个常见生信小问题(1-5)

前言:昨晚听了师兄开的博客熊言熊语,印象最深的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

运行结果:

image.png

再看看下面这个代码的输出和上面的有什么区别?

# 加载包
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

image.png

为什么出现了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?#

作者: 陈玫君

现任广州再生医学实验室研究实习员,苏州大学生物科学专业本科毕业,研究内容主要涉及分子生物学以及单细胞测序相关分析,曾实习于中科院上海植生所。



发表评论

您的电子邮箱地址不会被公开。

联系我们

联系我们

(44)07934433023

在线咨询: QQ交谈

邮箱: info@bioengx.org

关注微信
微信扫一扫关注我们

微信扫一扫关注我们

关注微博
返回顶部
Designed by

best down free | web phu nu so | toc dep 2017