【生信技能树】GEO数据库挖掘 P6 5了解矩阵

【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibilijimmy开始招学徒啦~这个课程就是周末的线下课的录屏哦~戳它了解→ https://mp.weixin.qq.com/s/gzyCRNnfgYkSsnjPSr-MGw 专栏：https://www.bilibili.com/read/cv719181这个系列视频包学包会，学不会……那就多看几遍……或者……欢迎报名jimmy的线下课呀~o(*￣▽￣*)ブhttps://www.bilibili.com/video/BV1is411H7Hq?p=6&spm_id_from=333.1007.top_right_bar_window_history.content.click本章为该section的笔记

【生信技能树】GEO数据库挖掘 P6 5 了解你的表达矩阵

# 检查matrix表达量是否合适
# 载入内参基因的表达量，看是否与相应的基因名相统一
# GAPDH 以及 ACTB（b-actin的基因名）是常用的内参
# 转换exprset的rownames为基因名
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]

# 检查matrix表达量是否合适
# 载入内参基因的表达量，看是否与相应的基因名相统一
# GAPDH 以及 ACTB（b-actin的基因名）是常用的内参
# 转换exprset的rownames为基因名
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]

后续检查样本是否正常表达一些基因。

【方法1】用GAPDH和ACTB基因作为阳性参照，看是否表达量正常。（应保守的高表达）

【方法2】用melt将样本变成纵列后，重命名以及予以分组信息，使用ggplot2看基因表达量的分布。

# 检查1，常见基因表达量，确定管家基因是否弄错
exprSet['GAPDH',]
exprSet['ACTB',]
# 检查2 检验基因分布图
if(T){
  library(reshape2)
  exprSet_L=melt(exprSet)
  colnames(exprSet_L)=c('probe','sample','value')
  # group_list获取每个样本的分组信息
  group_list=c('Control','Control','Control','Vemurafenib','Vemurafenib','Vemurafenib')
  exprSet_L$group=rep(group_list,each=nrow(exprSet))
  head(exprSet_L)
  ### ggplot2 
  library(ggplot2)
  p=ggplot(exprSet_L,
           aes(x=sample,y=value,fill=group))+geom_boxplot()
  print(p)
  p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
  print(p)
  p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
  print(p)
  p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
  print(p)
  p=ggplot(exprSet_L,aes(value,col=group))+geom_density() 
  print(p)
  p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
  p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
  p=p+theme_set(theme_set(theme_bw(base_size=20)))
  p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
  print(p)
  # 如果来自不同样本或基线不匹配，需使用sv包中的combine函数来矫正
  # 检查样本是否有问题，是否有异常值或不一致的分布
}

【方法3】用hclust包看以下样本间的关系，是否相同处理聚类到了一起。

grouplist可以根据网址上的样本信息予以确定，也可以根据table确定。

# 检查3 利用hclust包进行聚类，简单看一下进化关系，是否每个样本都类似
if(T){
  # pdata=pData(data_information)
  group_list=as.character(pdata[,2])
  group_list
  dim(exprSet)
  exprSet[1:5,1:5]
  ## hclust 
  colnames(exprSet)=paste(group_list,1:6,sep='')
  # Define nodePar
  nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                  cex = 0.7, col = "blue")
  hc=hclust(dist(t(exprSet)))
  par(mar=c(5,5,5,10)) 
  plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)
}

【方法4】PCA主成分分析看一下样本情况。

# 检查4 利用PCA来看样本分布是否一致
if(F){
  # BiocManager::install('ggfortify')
  library(ggfortify)
  df=as.data.frame(t(exprSet))
  df$group=group_list 
  autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
  
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra") 
  df=as.data.frame(t(exprSet))
  dat.pca <- PCA(df, graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind = group_list, # color by groups
               # palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
}