刘小泽写于2018.8.26，今天先做8道题
练习题来自生信技能树jimmy，http://www.bio-info-trainee.com/3409.html
1 安装R包
数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite(c("ALL","CLL", "pasilla", "airway")) #数据包
biocLite(c("limma","DESeq2", "clusterProfiler")) #软件包
install.packages("reshape2") #工具包
install.packages("ggplot2") #绘图包
2 了解ExpressionSet对象
CLL包中有data(sCLLex)，是一个表达芯片数据对象，其中包含许多信息
sCLLex数据包
第一行的ExpressionSet就是表达矩阵，查看它使用exprs ，用dim 查看矩阵大小
suppressPackageStartupMessages(library(CLL)) #隐藏具体加载信息
data("sCLLex")
sCLLex
e=exprs(sCLLex) #提取出来表达矩阵，赋给e
dim(e) #查看表达矩阵大小
#结果可以看到，包含了12625个探针，22个样本
str(e) # 查看结构
head(e) # 查看前6行
sampleNames(sCLLex) #查看样本编号
varMetadata(sCLLex) #查看所有表型变量
featureNames(sCLLex)[1:100] #查看基因芯片编码
使用str查看对象的结构，使用head查看对象的前6行（默认）
3 安装并了解hgu95av2.db包
- 官网：http://www.bioconductor.org/packages/release/data/annotation/html/hgu95av2.db.html
- 安装
biocLite("hgu95av2.db") -
hgu95av2.db包
这个数据库中共有36个包，每个包都可以当成一个列表操作，可以用as.list 函数展示数据，
library(hgu95av2.db)
ls("package:hgu95av2.db")
as.list(hgu95av2SYMBOL[1]) #变成列表进行查看第一组元素
#结果得到的就是
#$`1000_at`
#[1] "MAPK3"
#不懂hgu95av2SYMBOL是什么意思，就查找
?hgu95av2SYMBOL
#结果得到：hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations【就是基因id与探针号的关系，这里显示的就是MAPK3这个基因对应的探针id是1000_at】
tmp=toTable(hgu95av2SYMBOL) #toTable把一对多的list压缩成表格，存储在tmp中
colnames(tmp) #看一下tmp的列名
tmp[grep("^TP53$",tmp[,2]),] #查找TP53基因对应的探针id
-
探针与基因的对应关系
不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量
length(unique(tmp$symbol)) #查看总共基因数（注意unique的使用，为了避免重复基因出现） tail(sort(table(tmp$symbol))) #table输出基因id与探针的对应关系，sort从小到大排序，tail取最大6个 table(sort(table(tmp$symbol))) # 排序后，将探针数量与基因数对应，比如有8个探针的基因存在5个 -
找到sCLLex表达矩阵(e)在hgu95av2.db包中没有交叉的探针
#找探针id的包含关系就看hgu95av2.db包中的hgu95av2SYMBOL（上面的tmp） table(rownames(e)%in%tmp[,1]) #%%是管道符号，相当于linux的“|”，%in%表示两者求交集 #结果FALSE TRUE 1165 11460 -
过滤掉那些没有被hgu95av2.db收录的探针
e=e[rownames(e)%in%tmp[,1]]
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com
Welcome to our bioinfoplanet!
网友评论