GEO数据分析step2

最新推荐文章于 2024-09-15 12:50:10 发布

JiangCaifu

最新推荐文章于 2024-09-15 12:50:10 发布

阅读量1k

点赞数

分类专栏：生物信息

本文链接：https://blog.csdn.net/JiangCaifu/article/details/101304835

版权

生物信息专栏收录该内容

6 篇文章 6 订阅

订阅专栏

step2：检查数据

rm(list = ls())  ## 魔幻操作，一键清空~

options(stringsAsFactors = F)

load(file = 'step1-output.Rdata')

table(group_list)

# 每次都要检测数据

dat[1:4,1:4]

## 下面是画PCA的必须操作，需要看说明书。
#t表示转置
dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换

dat=as.data.frame(dat)#将matrix转换为data.frame
#cbind横向追加，即将分组信息追加到最后一列
dat=cbind(dat,group_list) 

library("FactoMineR")#画主成分分析图需要加载这两个包

library("factoextra") 

# The variable group_list (index = 54676) is removed

# before PCA analysis

dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)

fviz_pca_ind(dat.pca,
             
             geom.ind = "point", # show points only (nbut not "text")
             
             col.ind = dat$group_list, # color by groups
             
             # palette = c("#00AFBB", "#E7B800"),
             
             addEllipses = TRUE, # Concentration ellipses
             
             legend.title = "Groups"
             
)

ggsave('all_samples_PCA.png')

#删除当前工作环境中所用得对象
rm(list = ls())  ## 魔幻操作，一键清空~

load(file = 'step1-output.Rdata') #此步为一个小插曲，即计算一下从第一行开是计算每一行的sd值，知道最后一行所需要的时间

dat[1:4,1:4] 
#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
#tail函数：查看数据中的最后几行元素
#names对向量命名
#对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
cg=names(tail(sort(apply(dat,1,sd)),1000))
library(pheatmap)
#dat([cg,]).as.numeric([cg,])
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) 
# 'scale'可以对log-ratio数值进行归一化
n=t(scale(t(dat[cg,])))

n[n>2]=2 

n[n< -2]= -2

n[1:4,1:4]

pheatmap(n,show_colnames =F,show_rownames = F)
#class(group_list)
ac=data.frame(g=group_list)
#把ac的行名给到n的列名，即对每一个探针标记上分组信息
rownames(ac)=colnames(n) 

pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=ac,filename = 'heatmap_top1000_sd.png')