您的位置 首页 生物信息学

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

来源: https://www.yuque.com/docs/share/f52fe166-d1f4-469a-8a4a-7d1a05567fc5?# 发布: 2020-04-15 7,297阅读 84评论

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

701赞



相关文章

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注

评论列表(84)

联系我们

联系我们

(44)07934433023

在线咨询: QQ交谈

邮箱: info@bioengx.org

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

微信扫一扫关注我们

关注微博
返回顶部