GEO下载的数据在R中读取和注释

#  读取GEO下载的矩阵和注释
####  加载包 ----
library(readxl)
library(tidyverse)
library(GEOquery)# 读取 series_matrix文件 ----
GSE_data <- getGEO(filename = "data\\GSE14520-GPL571_series_matrix.txt.gz", getGPL = F)pd_GSE_data <- pData(GSE_data) #观察临床信息中的data processing,Microarray suite,MAS 5.0，即标准化方法，如果已经经过了MAS，一般情况下不需要再标准化了，避免矫枉过正
pd_GSE_data$data_processing[1]  #查看数据处理的方式
GSE_data_targets <-  pd_GSE_data %>%  # 提取有用的信息以便后续分析dplyr::select(sample_id = geo_accession,sample_name = title,tissue_type = `Tissue:ch1`)  GSE_data_expr <- exprs(GSE_data) %>% as.data.frame()  # 提取表达矩阵
# 数据集的标准化结果查看------
boxplot(exprs(GSE_data))  ## 箱线图查看数据是否标准化
## PCA  查看不同样本间的区别
PCA_new <- function(expr, ntop = 500, group, show_name = F){library(ggplot2)library(ggrepel)rv <- genefilter::rowVars(expr)select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]pca <- prcomp(t(expr[select, ]))#最核心的代码percentVar <- pca$sdev^2/sum(pca$sdev^2)d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], group = group, name = colnames(expr))attr(d, "percentVar") <- percentVar[1:2]if (show_name) {ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group")) + geom_point(size = 2) +xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) +geom_text_repel(aes(label = name),size = 3,segment.color = "black",show.legend = FALSE )} else {ggplot(data = d, aes_string(x = "PC1", y = "PC2",color = "group")) + geom_point(size = 2) +xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance"))}
}
PCA_new(exprs(GSE_data), nrow(GSE_data),group = GSE_data_targets$tissue_type)  # 组别以颜色显示## getGEO 读取注释信息 ----
##　AnnotGPL = T，用的是.gz的文件
GPL_data<- getGEO(filename ="data//GPL571.annot.gz", AnnotGPL = T)
GPL_data_11 <- Table(GPL_data)## AnnotGPL = F，用的是.soft的文件
GPL_data_2 <- getGEO(filename = "data//GPL10687_family.soft.gz", AnnotGPL = F)
GPL_data_21 <- Table(GPL_data_2)# 进行注释------
GSE_data_expr <- GSE_data_expr %>% rownames_to_column() %>% reshape::rename(c(rowname = "prode_id"))GPL <- GPL_data_11 %>% dplyr::select(ID,`Gene symbol`) %>% reshape::rename(c(ID = "prode_id",`Gene symbol` = "GeneSymbol")) %>% add_count(GeneSymbol,name = "n_gene") %>% add_count(prode_id,name = "n_prode_id") %>% dplyr::filter(n_gene < 1123) %>% # 选择能有对应基因的探针dplyr::select(c(-3,-4))GSE_expr <- GSE_data_expr %>% inner_join(GPL, ., by = "prode_id") %>% #p2s和上一步得到的结果再取交际，p2s放在右边是以它为准dplyr::select(-1) %>%  #去除第一列probe_idas.data.frame() %>% #因为aggregate需要数据框格式aggregate(.~ GeneSymbol, data = ., mean) %>%  #以symbol作为因子水平，把相似的数据放在一起取均值，最大值max，中位值mediancolumn_to_rownames(var = "GeneSymbol")## 可视化注释后的矩阵，因为注释过程丢失了部分探针
boxplot(GSE_expr)
PCA_new(GSE_expr, nrow(GSE_expr),group = GSE_data_targets$tissue_type)  # 组别以颜色显示
# 以RData 格式保存信息 ----
GSE14520_GPL571 <- GSE_data
GSE14520_GPL571_expr <- GSE_expr
GPL571 <- GPL_data_11
GSE14520_GPL571_targets <- GSE_data_targetspicDir <- "dealed_data"
dir.create(picDir )
save(GSE14520_GPL571,  # 矩阵GSE14520_GPL571_targets, # metadataGPL571, #注释信息file = paste('.',picDir,"GEOqury读取的GSE14520-GPL571_data.RData",sep="/"))
save(GSE14520_GPL571_expr,file = paste('.',picDir,"GSE14520-GPL571_expr.RData",sep="/"))

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！