生信学习——R语言小作业-中级（附详细答案解读）

写在前面——自从上次做了初级的题目之后，就一直在看这个中级的题目。因为中间有事耽搁了许久，所以间隔了很多天才做完。虽然按照视频和百度磕磕绊绊的把这个题目写完了，但是脑子还是一团浆糊。知道代码是干嘛的，但是不知道为什么要这么做。革命尚未成功，同志仍需努力…

题目原文：http://www.bio-info-trainee.com/3750.html
视频教程：https://www.bilibili.com/video/BV1cs411j75B?p=13
优质答案：https://www.jianshu.com/p/e15ee2cd3174

注意：如果library(…)报错的话，是因为没有安装包，需要install.packegs(…)安装对应的包。

1. 请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)。

ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

# 清空Environment中的变量，也可以直接点小扫帚清空
rm(list = ls())

# 避免把字符串项当成因子
options(stringsAsFactors = FALSE)

a <- read.table("practice/e1.txt")

# 思路：先得到egSYMBOL和egENSEMBL数据框
library(org.Hs.eg.db)
g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)

# 保留a中的v1小数点前面的部分，并将其赋给a的ensembl_id
a$ensembl_id <- unlist(lapply(a$V1, function(x){
   
  strsplit(x, "[.]")[[1]][1]
}))

# 把a和g2e通过ensembl融合
tmp <- merge(a, g2e, by = "ensembl_id")
# 最后再根据gene_id进行融合
result <- merge(tmp, g2s, by = "gene_id")

2. 根据R包hgu133a.db找到下面探针对应的基因名(symbol)。

1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_at

rm(list = ls())
options(stringsAsFactors = FALSE)
b <-  read.table("practice/e2.txt")

# 安装hgu133a.db包
# 若出错请查看 https://blog.csdn.net/narutodzx/article/details/119378949
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("hgu133a.db")

library(hgu133a.db)
ids <- toTable(hgu133aSYMBOL)
head(ids)

colnames(b) <- "probe_id"
result <- merge(ids, b, by="probe_id")

​

3. 找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图。想想如何通过 ggpubr 进行美化。

rm(list = ls())
options(stringsAsFactors = FALSE)

# BiocManager::install("CLL")
# BiocManager::install("hgu95av2.db")

# 隐藏包在加载的时候显示的信息
suppressPackageStartupMessages(library(CLL))
# 加载数据
data(sCLLex)
sCLLex

# 获得表达矩阵
# 用expr()提取assayData信息
exprSet <- exprs(sCLLex) 
# 获得临床信息
# 用pData()提取phenoData信息
pd <- pData(sCLLex)

library(hgu95av2.db)
ids <- toTable(hgu95av2SYMBOL)
head(ids)

# 在ids中搜索TP53

# 绘图
boxplot(exprSet['1939_at',] ~ pd$Disease) #signal
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)