对20个R语言习题的解答和思考之一
国内学习生信分析，肯定逃不开“生信技能树”的资源，那就翻翻历朝资源，从最基本的20个R语言习题开始吧。虽然已经有官方提供的参考答案，网上也已经有很多前辈分享各种坑和知识点，但只有自己亲身实践，才是自己的。
- 准备工作
rm(list=ls(all=TRUE))
setwd("~/Code/R20") # 或者建立Project
- 安装一些R包：
数据包： ALL, CLL, pasilla, airway 软件包：limma，DESeq2，clusterProfiler 工具包：reshape2 绘图包：ggplot2
不同领域的R包使用频率不一样，在生物信息学领域，尤其需要掌握bioconductor系列包。
一般两种方式：
# 从CRAN下载
install.packages('[package-name]')
# 从bioconductor下载
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install('[package-name]')
# 比如：安装"hgu95av2.db"：需要先安装Biostrings，后安装org.Hs.eg.db，再安装hgu95av2.db
其他：
为提高从CRAN或bioconductor源镜像下载稳定性，按目前网络情况建议改成西湖大学的源，清华源也如此推荐 ,修改方法：
镜像源配置文件之一是 .Rprofile（mac下位于 ~/.Rprofile, Windows 下位于 X:\Programs\R\R-4.4.2\etc\Rprofile.site）。在文末添加如下语句：
options("repos"= c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"));
options(BioC_mirror="<https://mirrors.westlake.edu.cn/bioconductor/>")
-
了解ExpressionSet对象，比如CLL包里面就有data(sCLLex) ，找到它包含的元素，提取其表达矩阵(使用exprs函数)，查看其大小
ExpressionSet类似于AnnData，一种数据集的封装标准。
library(CLL)
# 从sCLLex这个ExpressionSet数据集中提取必要信息
#ExpressionSet/SummarizedExperiment 类的数据容器，包含了基因表达矩阵及其注释
data("sCLLex")
exprSet = exprs(sCLLex) # 提取表达矩阵，数据结构：矩阵，12625行（基因）*22列（样本）
samples = sampleNames(sCLLex) # 提取样本清单，数据结构：列表，22个（样本）
pData = pData(sCLLex) # 提取pheno数据，数据结构：矩阵，22行（样本）*2列（SampleID + Disease）
group_list = as.character(pData[,2]) # 从pheno数据提取临床分组信息，数据结构：列表，22个（分组归类）
dim(exprSet) # 了解表达矩阵的维度，22个sample，12625个探针
- 了解 str,head,help函数，作用于 第二步提取到的表达矩阵
# 展示数据的structure，对于复合数据结构，展示内部元素结构、数据类型等等
str(sCLLex)
# 展示数据的前面部分的内容，R中head默认取6行
head(sCLLex)
# 展示数据的统计量，长度、数据类型、模式
summary(sCLLex)
- 安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后显示的那些变量
生信技能树曾对此挖掘过用法，Bioconductor系列之hgu95av2.db
library(hgu95av2.db)
ls("package:hgu95av2.db") # db包括36个映射数据，把芯片探针ID号转为其它36种主流ID号的映射
columns(hgu95av2.db) # 查看芯片注释包提供的注释信息
head(toTable(hgu95av2SYMBOL)) # 提取注释信息，这里选择GENESYMBOL，也可以选择其他ID进行注释（选择其他包）
- 理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID
“hgu95av2.db”是平台型注释包，使用Affymetrix公司的hgu95av2基因芯片得到的探针表达，都可以用这个数据包hgu95av2是芯片探针ID映射为基因在其他体系的ID或名字。
其中，hgu95av2SYMBOL是探针名probe_id与基因名symbol的映射数据库。并且通过数据探索可知，数据库总共有12625个probe id，设置了与SYMBOL映射关系的有11683个，因多id映射到同基因的关系，共映射到8776个基因SYMBOL
# ids是取hgu95av2SYMBOL核心数据内容形成table
ids <- toTable(hgu95av2SYMBOL)
library(dplyr)
filter(ids, symbol == "TP53")
# TP53对应三个基因探针。
# 探索hgu95av2SYMBOL。这个数据库记录了probe ID和GENE SYMBOL的映射关系
summary(hgu95av2SYMBOL)
# SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
# |
# | Lkeyname: probe_id (Ltablename: probes)
# | Lkeys: "1000_at", "1001_at", ... (total=12625/mapped=11683)
# |
# | Rkeyname: symbol (Rtablename: gene_info)
# | Rkeys: "A1BG", "A2M", ... (total=193329/mapped=8776)
# |
# | direction: L --> R
x <- hgu95av2SYMBOL
mapped_genes <- mappedkeys(x)
xx <- as.list(x[mapped_genes])
# length(x)=12625
# length(xx)=11683
if(length(xx) > 0) {
xx[1:5]
xx[[1]]
}
length(unique(xx)) # 8776
- 理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？
不管是Agilent芯片，还是Affymetrix芯片，设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb，因此，会出现多个探针对应一个基因。所以，可以取最大表达值的探针来作为基因的表达量。
# ids（hgu95av2SYMBOL）表有两列，一列probe_id，一列symbol
# probe_id值是唯一的，有11683个，symbol不是唯一的，唯一值有8776个
length(ids$symbol) # 11683
length(ids$probe_id) # 11683
length(unique(ids$probe_id))# 11683
length(unique(ids$symbol)) # 8776
# 将symbol出现次数汇总制表并排序（默认升序），可知有4个基因带8个探针，6个基因带7个探针
tail(sort(table(ids$symbol)),n=10L)
# ACTB CD44 ERBB3 GAPDH RBPMS TERF1 INPP4A MYB PTGER3 STAT1
# 7 7 7 7 7 7 8 8 8 8
table(sort(table(ids$symbol)))
# 1 2 3 4 5 6 7 8
# 6725 1444 451 103 27 16 6 4
plot(table(sort(table(ids$symbol))))
- 第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。
exprSet的rowname如果不出现在hgu95av2的probe_id列，就无法做相应的ID匹配。 匹配后可发现有942个探针无法找到SYMBOL名字匹配。
table(rownames(exprSet) %in% ids$probe_id)
n_exprSet <- exprSet[!(rownames(exprSet) %in% ids$probe_id), ]
dim(n_exprSet)
- 过滤表达矩阵，删除那1165个没有对应基因名字的探针。
exprSet <- exprSet[rownames(exprSet) %in% ids$probe_id, ]
dim(exprSet)
# 法2
length(hgu95av2SYMBOL)
# [1] 12625
probe_map <- hgu95av2SYMBOL
length(probe_map)
#全部的探针数目
# [1] 12625
probe_info <- mappedkeys(probe_map)
length(probe_info)
#探针与基因产生映射的数目
# [1] 11683
gene_info <- as.list(probe_map[probe_info])
# 转化为数据表
length(gene_info)
head(gene_info)
# [1] 11683
gene_symbol <- toTable(probe_map[probe_info])
# 从hgu95av2SYMBOL文件中，取出有映射关系的探针，并生成数据框给gene_symbol
head(gene_symbol)
# probe_id symbol
# 1 1000_at MAPK3
# 2 1001_at TIE1
# 3 1002_f_at CYP2C19
# 4 1003_s_at CXCR5
# 5 1004_at CXCR5
# 6 1005_at DUSP1
- 整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。
- 提示，理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
- 然后根据得到探针去过滤原始表达矩阵
by(data, INDICES, FUN)函数的典型用法，将data数据框或矩阵按照INDICES因子水平进行分组，然后对每组应用FUN函数。
本例应用by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] ) - 对表达矩阵exprSet进行分组，将同一个symbol所对应的多个探针分成不同的组，并对每组探针进行统计：计算每组中每行探针表达量的平均值（也就是每个探针在6个样本中表达量的均值rowMeans(x)），再取平均值最大的那个探针作为该symbol所对应的唯一探针 - 第二个参数ids$symbol定义了分组，将第一参数—exp表达矩阵分成了若干个小矩阵，每个小矩阵里存放着同一个symbol所对应的所有探针。 - 第三个参数是我们自己定义的函数：计算每个小矩阵中每行探针表达量的平均值（也就是每个探针在6个样本中表达量的均值rowMeans(x)），再取平均值最大的那个探针作为该symbol所对应的唯一探针which.max(rowMeans(x)) - 参考：https://cloud.tencent.com/developer/article/1431357
# 整理数据：
# ids（处理前）：hgu95av2的芯片相关所有数据，probe_id和symbol的对应关系
# ids（处理后）：确认仅在本例数据库中出现的probe
dim(ids)
ids=ids[match(rownames(exprSet),ids$probe_id),]
dim(ids)
# 处理exprSet，去除没有对应symbol命名的探针（在ids中没有对应probe_id的）
# 根据symbol将exprSet分组，并计算各组均值，取出均值最大行的行名（探针名）
tmp = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
head(probes)
# [1] "36512_at" "39463_at" "38434_at" "36332_at" "36185_at" "35761_at"
# 仅保留上述探针所在的行，剩8776行
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
- 把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。
rownames(exprSet) <- ids[match(rownames(exprSet), ids$probe_id), 2]
head(exprSet)
# 有几个管家基因
exprSet['GAPDH',]
boxplot(exprSet[,1])
boxplot(exprSet['GAPDH',])
exprSet['ACTB',]
浙公网安备 33010602011771号