043 判断Geo数据是否需要进行log处理

实验目的

实验代码

###加载R包
library(readxl)
library(tidyverse)
library(GEOquery)
library(tidyverse)
library(GEOquery)
library(limma) 
library(affy)
library(stringr)
library(FactoMineR)
library(factoextra)
library(sva)
###下载数据，如果文件夹中有会直接读入
gset = getGEO('GSE205185', destdir=".", AnnotGPL = T, getGPL = T)
class(gset)
gset1
gset[["GSE205185_series_matrix.txt.gz"]]@phenoData@data[["data_processing"]]
gset2 = getGEO('GSE29431', destdir=".", AnnotGPL = T, getGPL = T)
class(gset2)
gset21
gset2[["GSE29431_series_matrix.txt.gz"]]@phenoData@data[["data_processing"]]
gset3 = getGEO('GSE20711', destdir=".", AnnotGPL = T, getGPL = T)
class(gset3)
gset31
gset3[["GSE20711_series_matrix.txt.gz"]]@phenoData@data[["data_processing"]]
#提取子集
plf1<-gset1@annotation
plf2<-gset21@annotation
plf3<-gset31@annotation
#提取平台文件
GPL_data<- getGEO(filename ="GPL21185.soft.gz", AnnotGPL = T)
GPL_data_11 <- Table(GPL_data)
GPL_data1<- getGEO(filename ="GPL570.annot.gz", AnnotGPL = T)
GPL_data_22 <- Table(GPL_data1)
GPL_data2<- getGEO(filename ="GPL570.annot.gz", AnnotGPL = T)
GPL_data_33 <- Table(GPL_data2)
###提取表达量
exp <- exprs(gset1)
###使用boxplot函数查看
boxplot(exp)
###使用判别函数判断
ex <- exp
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <-(qx[5]>100)||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2]>0&&qx[2]<1&&qx[4]>1&&qx[4]<2)
if(LogC){
  ex[which(ex<=0)]<-NaN
  exprSet <- log2(ex)
  print("需要取log2")}else{print("无需取log2")
  }
probe_name<-rownames(exp)

###提取表达量
exp2 <- exprs(gset21)
###使用boxplot函数查看
boxplot(exp2)
###使用判别函数判断
ex <- exp2
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <-(qx[5]>100)||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2]>0&&qx[2]<1&&qx[4]>1&&qx[4]<2)
if(LogC){
  ex[which(ex<=0)]<-NaN
  exprSet <- log2(ex)
  print("需要取log2")}else{print("无需取log2")
  }
probe_name2<-rownames(exp2)

###提取表达量
exp3 <- exprs(gset31)
###使用boxplot函数查看
boxplot(exp3)
###使用判别函数判断
ex <- exp3
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <-(qx[5]>100)||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2]>0&&qx[2]<1&&qx[4]>1&&qx[4]<2)
if(LogC){
  ex[which(ex<=0)]<-NaN
  exprSet <- log2(ex)
  print("需要取log2")}else{print("无需取log2")
  }
probe_name3<-rownames(exp3)

#我改进的代码
setwd("F:/001微生物课题/公共数据库生信挖掘/02.Hp感染24小时后AGS的差异基因/001.判断GEO是否需要log2") ##💚💚💚💚💚设置路径



### 加载所需的R包
# readxl: 用于读取Excel文件
library(readxl)
# tidyverse: 数据处理和可视化的综合包
library(tidyverse)
# GEOquery: 用于从GEO数据库下载和处理数据
library(GEOquery)
library(tidyverse)
library(GEOquery)
# limma: 用于差异表达分析
library(limma) 
# affy: 用于处理Affymetrix芯片数据
library(affy)
# stringr: 用于字符串处理
library(stringr)
# FactoMineR和factoextra: 用于主成分分析和可视化
library(FactoMineR)
library(factoextra)
# sva: 用于批次效应处理
library(sva)

### 下载并读取GEO数据集💚💚💚💚💚设置自己的GEO号
# 使用getGEO函数下载GSE70394数据集
# destdir="."表示下载到当前目录
# AnnotGPL = T表示下载注释信息
# getGPL = T表示同时下载平台信息
gset = getGEO('GSE70394', destdir=".", AnnotGPL = T, getGPL = T)
class(gset)  # 查看数据类型
gset[[1]]    # 查看第一个数据集的信息
# 查看数据处理方法信息
gset[["GSE74577_series_matrix.txt.gz"]]@phenoData@data[["data_processing"]]


# 获取每个数据集的平台注释信息
plf1<-gset[[1]]@annotation



### 处理第一个数据集(GSE70394)
# 使用exprs函数提取表达矩阵
exp <- exprs(gset[[1]])
# 使用箱线图可视化数据分布情况
boxplot(exp)
# 使用分位数方法判断是否需要log2转换
ex <- exp
# 计算数据的多个分位数点
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
# 判断是否需要log2转换的条件：
# 1. 99%分位数大于100 或
# 2. 数据范围大于50且25%分位数大于0 或
# 3. 25%分位数在0-1之间且75%分位数在1-2之间
LogC <-(qx[5]>100)||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2]>0&&qx[2]<1&&qx[4]>1&&qx[4]<2)
if(LogC){
  ex[which(ex<=0)]<-NaN  # 将小于等于0的值设为NaN
  exprSet <- log2(ex)    # 进行log2转换
  print("需要取log2")}else{print("不需要取log2")
  }
probe_name<-rownames(exp)  # 保存探针名称

最终的结果

image.png