6.探针与基因的关系
原题目描述：
理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？
length(unique(ids$symbol))
# [1] 8585
tail(sort(table(ids$symbol)))
#YME1L1 GAPDH INPP4A MYB PTGER3 STAT1
# 7 8 8 8 8 8
table(sort(table(ids$symbol)))
# 1 2 3 4 5 6 7 8
# 6555 1428 451 102 22 16 6 5
答：不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量。
7.提取目的基因
原题目描述：
第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。
#芯片表达的数据存放在data_expression这个变量中，它有22个样本的12625个探针数据，如下所示：
class(data_expression)
# [1] "matrix"
str(data_expression)
# num [1:12625, 1:22] 5.74 2.29 3.31 1.09 7.54 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...
# ..$ : chr [1:22] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" ...
现在要解决的一个问题是：在data_expression中有一些基因没有收录在hug95av2.db包中，现在要找到这些基因，把它们都去掉。
思路：找到hug95av2.db包中能够与基因映射起来的探针，再找到这些探针对应的基因SYMBOL，然后在data_expression中取出这些基因即可。此时要使用到hug95av2.db包中的hug95av2SYMBOL这个文件，如下所示：
?hug95av2SYMBOL
查阅其文档，了解到如下信息：hug95av2SYMBOL是一个R对象，它提供的是芯片生产厂家与基因缩写之间的映射信息。这个映射的信息主要依据Entrez Gene数据库。现在我们通过mappedkeys()这个函数，得到映射到基因上的探针信息，如下所示：
probe_map <- hgu95av2SYMBOL
length(probe_map)
#全部的探针数目
# [1] 12625
probe_info <- mappedkeys(probe_map)
length(probe_info)
#探针与基因产生映射的数目
# [1] 11471
gene_info <- as.list(probe_map[probe_info])
# 转化为数据表
length(gene_info)
# [1] 11471
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
###其中，gene_symbol就是有映射关系的基因，它里面的数据是探针的编码以及探针对应的基因，在data_expression中它们挑出来即可。这里使用到了mappedkeys()函数，这个函数用于处理映射文件(bimap)，它的参数就是一个Bimap对象，例如上面的hgu95av2SYMBOL这个对象，关于Bimap对象，可以看这篇笔记《Bimap对象笔记》。
8. 过滤表达矩阵
原问题描述：过滤表达矩阵，删除那1165个没有对应基因名字的探针。这时再看一下原始数据，如下所示：
head(data_expression)
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL CLL16.CEL CLL17.CEL
# 1000_at 5.743132 6.219412 5.523328 5.340477 5.229904 4.920686 5.325348
# 1001_at 2.285143 2.291229 2.287986 2.295313 2.662170 2.278040 2.350796
# 1002_f_at 3.309294 3.318466 3.354423 3.327130 3.365113 3.568353 3.502440
# 1003_s_at 1.085264 1.117288 1.084010 1.103217 1.074243 1.073097 1.091264
# 1004_at 7.544884 7.671801 7.474025 7.152482 6.902932 7.368660 6.456285
# 1005_at 5.083793 7.610593 7.631311 6.518594 5.059087 4.855161 5.176975
思路：在这个原始数据中，有12625个探针，下面通过isNA函数可以得到一个逻辑向量，其中TRUE表示没有映射关系的探针，FALSE表示有映射关系的探针。然后使用seq生成一个整数向量，并提取那些有映射关系的探针所在的列，提取原始数据即可，如下所示：data_filter <- data_expression[rownames(data_expression)%in%probe_info]
data_filter <- data_expression[rownames(data_expression)%in%probe_info,]
# rownames(data_expression)是原始数据表达值的名称，这里提探针的编号
# probe_info是有映射的探针的编号
# A%in%B，表示A中有哪些成员在B中，有的返回TRUE，无的返回FALSE
# 这行命令表示，提取那些在probe_info中的探针的表达值，赋值给变量data_filter
table(rownames(data_expression) %in% probe_info)
# FALSE TRUE
# 1154 11471
# 从上面的数据可以看出来，原始表达值中有11471个数据在probe_info中
table(rownames(data_expression) %in% rownames(data_filter))
# FALSE TRUE
# 1154 11471
现在就得了过滤后的表达矩阵，即data_filter。
9.整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。
A. 提示，理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
B. 然后根据得到探针去过滤原始表达矩阵
ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
View(head(exprSet))
10.把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# MAPK3 5.743132 6.219412 5.523328 5.340477 5.229904
# TIE1 2.285143 2.291229 2.287986 2.295313 2.662170
# CYP2C19 3.309294 3.318466 3.354423 3.327130 3.365113
# CXCR5 1.085264 1.117288 1.084010 1.103217 1.074243
# CXCR5 7.544884 7.671801 7.474025 7.152482 6.902932
library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
# probe sample value group
#1 MAPK3 CLL11.CEL 5.743132 progres.
#2 TIE1 CLL11.CEL 2.285143 progres.
#3 CYP2C19 CLL11.CEL 3.309294 progres.
#4 CXCR5 CLL11.CEL 1.085264 progres.
#5 CXCR5 CLL11.CEL 7.544884 progres.
#6 DUSP1 CLL11.CEL 5.083793 progres.
View(head(exprSet))
