11 ID转换——表达矩阵的整理

cleaning 选择表达量最大的基因

### 这里可以看到同一个探针可以对应多个ID,所以要先对表达矩阵进行过滤，保留表达量最大的那个ID
load(file = './Rdata/step0.Rdata')
### 首先要有个ID2Seq的矩阵
ID2Seq<- data.frame(ID=gpl18109$ID,Seq=gpl18109$SEQUENCE)
GSE53625_expr<- GSE53625_expr[GSE53625_expr$ID_REF%in%ID2Seq$ID,]
### 保证ID2Seq和GSE53625_expr 顺序相同
ID2Seq<- ID2Seq[match(GSE53625_expr$ID_REF,ID2Seq$ID),]
match(GSE53625_expr$ID_REF,ID2Seq$ID)
row.names(GSE53625_expr)<- GSE53625_expr$ID_REF
GSE53625_expr<- GSE53625_expr[,-1]
ID2Seq$max<- apply(GSE53625_expr,1,max)  #### 每一行的最大值，也就是每个基因在每个病人中的最大值
ID2Seq<- ID2Seq[order(ID2Seq$Seq,ID2Seq$max,decreasing = T),]
ID2Seq<- ID2Seq[!duplicated(ID2Seq$Seq),]
##### 最后完成表达矩阵的筛选
GSE53625_expr<- GSE53625_expr[row.names(GSE53625_expr)%in%ID2Seq$ID,]

去除没有注释的数据集

{
  exprSet = exprSet[ rownames(exprSet) %in% ID2gene[ , 1 ], ]
  ID2gene = ID2gene[ match(rownames(exprSet), ID2gene[ , 1 ] ), ]
}

dim( exprSet )
dim( ID2gene )
tail( sort( table( ID2gene[ , 2 ] ) ), n = 12L )

相同基因的表达数据取最大值

{
  MAX = by( exprSet, ID2gene[ , 2 ], 
            function(x) rownames(x)[ which.max( rowMeans(x) ) ] )
  MAX = as.character(MAX)
  exprSet = exprSet[ rownames(exprSet) %in% MAX , ]
  rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]
}
exprSet = log(exprSet)
dim(exprSet)
exprSet[1:5,1:5]

save( exprSet, group_list, file = 'exprSet_for_WGCNA.Rdata')

选择表达量最大的探针which.max(rowMeans(x)) 计算表达量,by()函数进行循环这个很关键

ID_trans <- function(exprSet,ids){
  tmp = by(exprSet, ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))] )
  probes = as.character(tmp)
  print(dim(exprSet))
  exprSet=exprSet[rownames(exprSet) %in% probes ,]
  
  print(dim(exprSet))
  rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
  return(exprSet)
}


new_exprSet <- ID_trans(exprSet,ids)

分步法解释选择表达量最大的探针

主要是利用了by这个函数，by 函数可以进行分类总汇

下面的这个代码就是ID_trans这个函数的分解版

利用by函数相当于批量对ids$symbol进行分类总汇
下面这些话是理解这个by 函数具体做了什么
exprSet[ids[,2]=='A2M',]
x= exprSet[ids[,2]=='A2M',]
row.names(x)
rowMeans(x)
which.max(rowMeans(x))
按照ids 分类

image.png image.png

# 确保probe_id是一致的
match(rownames(exprSet),ids$probe_id)
### 下面这些话是理解这个by 函数具体做了什么
# exprSet[ids[,2]=='A2M',]
# x= exprSet[ids[,2]=='A2M',]
# row.names(x)
# rowMeans(x)
# which.max(rowMeans(x))
# 按照ids 分类
# by是分类汇总，这步求的是基因表达量的均值
tmp= by(exprSet,ids$symbol,function(x)rownames(x)[which.max(rowMeans(x))])
probes<- as.character(tmp)
dim(exprSet)
exprSet= exprSet[rownames(exprSet)%in%probes,]
dim(exprSet)
exprSet$probe_id<- row.names(exprSet)
new_exprSet<- merge(ids,exprSet,by='probe_id')
rownames(new_exprSet)<- new_exprSet[,2]
new_exprSet<- new_exprSet[,-c(1,2)]

继续加油