实验目的

代码

###加载R包
library(readxl)
library(tidyverse)
library(GEOquery)
library(tidyverse)
library(GEOquery)
library(limma) 
library(affy)
library(stringr)
###下载数据，如果文件夹中有会直接读入
gset = getGEO('GSE205185', destdir=".", AnnotGPL = T, getGPL = T)
class(gset)
gset1

#读取平台文件
GPL_data<- getGEO(filename ="GPL21185.soft.gz", AnnotGPL = T)
GPL_data_11 <- Table(GPL_data)

#提取表达量
exp <- exprs(gset1)
probe_name<-rownames(exp)

#转换ID
loc<-match(GPL_data_11[,1],probe_name)
probe_exp<-exp[loc,]
raw_geneid<-(as.matrix(GPL_data_11[,"GENE_SYMBOL"]))
index<-which(!is.na(raw_geneid))
geneid<-raw_geneid[index]
exp_matrix<-probe_exp[index,]
geneidfactor<-factor(geneid)
gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean))
rownames(gene_exp_matrix)<-levels(geneidfactor)
gene_exp_matrix=na.omit(gene_exp_matrix)

#####读取分组信息#####
pdata <- pData(gset1)
group_list <- ifelse(str_detect(pdata$source_name_ch1,"primary breast tumour"), "T",
                     "N")
group_list
group_list = factor(group_list,
                    levels = c("T","N"))
group_list
pdata$group=group_list

#####进行数据矫正#####

boxplot(gene_exp_matrix,outline=T, notch=T,col=group_list, las=2)
dev.off()
gene_exp_matrix_noemal=normalizeBetweenArrays(gene_exp_matrix)
boxplot(gene_exp_matrix_noemal,outline=T, notch=T,col=group_list, las=2)
range(gene_exp_matrix_noemal)
gene_exp_matrix_noemal <- log2(gene_exp_matrix_noemal+1)
gene_exp_matrix_noemal <-as.data.frame(gene_exp_matrix_noemal)
gene_exp_matrix_noemal=na.omit(gene_exp_matrix_noemal)#去除NA列
write.csv(gene_exp_matrix_noemal,file = "geo_exp.csv")
range(gene_exp_matrix_noemal)
dev.off()

#####进行差异分析#####
design=model.matrix(~group_list)
fit=lmFit(gene_exp_matrix_noemal,design)#这里要注意，分组的样本与矩阵样本是否相符，不相符则去文件中调整然后读入
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
write.table(deg, file = "deg_all.txt",sep = "\t",row.names = T,col.names = NA,quote = F)
logFC=1
adj.P.Val = 0.05
k1 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC < -logFC)
k2 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC > logFC)
deg$change = ifelse(k1,"down",ifelse(k2,"up","stable"))
table(deg$change)
write.csv(deg,file="upanddown.csv")

#这个是我针对GEO的幽门感染GES-1的测序数据进行写的代码。

setwd("F:/GEO数据处理") ##💚💚💚💚💚设置路径




###加载R包
# readxl: 用于读取Excel文件
# tidyverse: 数据处理和可视化的综合工具包
# GEOquery: 用于从GEO数据库下载和处理数据
# limma: 用于差异表达分析
# affy: 用于处理Affymetrix微阵列数据
# stringr: 字符串处理工具包
library(readxl)
library(tidyverse)
library(GEOquery)
library(tidyverse)
library(GEOquery)
library(limma) 
library(affy)
library(stringr)

###下载数据：可在文件夹中或直接读取
# 从GEO数据库下载GSE74577数据集，并保存在当前目录
# AnnotGPL=T表示下载注释信息，getGPL=T表示下载平台信息
gset = getGEO('GSE74577', destdir=".", AnnotGPL = T, getGPL = T)
class(gset)  # 查看数据类型
gset[[1]]    # 查看第一个数据集的详细信息💚💚💚运行完到这里以后可以修改下面的GPL17586.soft.gz

#读取平台文件
# 获取GPL17586.soft.gz平台的注释信息
GPL_data<- getGEO(filename ="GPL17586.soft.gz", AnnotGPL = T)
# 提取平台注释表格数据
GPL_data_11 <- Table(GPL_data)

#读取表达矩阵
# 从数据集中提取表达矩阵
exp <- exprs(gset[[1]])
# 获取探针名称（行名）
probe_name<-rownames(exp)

#转换ID
# 匹配探针ID，确保表达矩阵和注释信息对应
loc<-match(GPL_data_11[,1],probe_name)
# 获取匹配后的表达矩阵
probe_exp<-exp[loc,]
# 从gene_assignment列获取基因信息并处理
raw_geneid <- as.matrix(GPL_data_11[,"gene_assignment"])
# 从gene_assignment列中提取基因名称（通常是第二个"//"之前的部分）
raw_geneid <- sapply(strsplit(raw_geneid, "//"), function(x) {
  if(length(x) >= 2) {
    # 去除空格并获取第一个有效的基因名
    gene <- trimws(x[2])
    if(gene == "---" || gene == "") return(NA)
    return(gene)
  } else {
    return(NA)
  }
})

# 继续后续处理
index <- which(!is.na(raw_geneid))
geneid <- raw_geneid[index]
exp_matrix <- probe_exp[index,]
geneidfactor<-factor(geneid)
gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean))
rownames(gene_exp_matrix)<-levels(geneidfactor)
gene_exp_matrix=na.omit(gene_exp_matrix)

#####读取分组信息#####
# 获取样本表型数据
pdata <- pData(gset[[1]])
# 根据样本名称确定分组信息：Hp感染组和未感染对照组
group_list <- ifelse(str_detect(pdata$title, "noninfected"), "Control",
                     "Hp")
print(group_list)  # 打印查看分组结果
# 将分组信息转换为因子类型，设定水平顺序（对照组在前，处理组在后）
group_list = factor(group_list,
                    levels = c("Control", "Hp"))
print(group_list)  # 打印查看分组结果
# 将分组信息添加到表型数据中
pdata$group = group_list

#####数据标准化处理#####
# 绘制标准化前的箱线图
boxplot(gene_exp_matrix,outline=T, notch=T,col=group_list, las=2)
dev.off()
# 进行数组间标准化处理
gene_exp_matrix_noemal=normalizeBetweenArrays(gene_exp_matrix)
# 绘制标准化后的箱线图
boxplot(gene_exp_matrix_noemal,outline=T, notch=T,col=group_list, las=2)
# 查看数据范围
range(gene_exp_matrix_noemal)
# 转换为数据框格式
gene_exp_matrix_noemal <-as.data.frame(gene_exp_matrix_noemal)
# 去除含NA的行
gene_exp_matrix_noemal=na.omit(gene_exp_matrix_noemal)
# 保存标准化后的表达矩阵
write.csv(gene_exp_matrix_noemal,file = "geo_exp.csv")
dev.off()

#####进行差异分析#####
# 构建设计矩阵
design=model.matrix(~group_list)
# 拟合线性模型
fit=lmFit(gene_exp_matrix_noemal,design)
# 贝叶斯检验，计算经验贝叶斯统计量
fit=eBayes(fit)
# 获取所有基因的差异分析结果
deg=topTable(fit,coef=2,number = Inf)
# 保存所有差异分析结果
write.table(deg, file = "deg_all.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

# 设置差异表达基因的筛选标准
logFC=1           # log2倍变化阈值
adj.P.Val = 0.05  # 调整后的p值阈值

# 标记上调和下调的基因
k1 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC < -logFC)    # 下调基因
k2 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC > logFC)     # 上调基因
# 添加变化方向标记
deg$change = ifelse(k1,"down",ifelse(k2,"up","stable"))
# 统计上调、下调和稳定的基因数量
table(deg$change)
# 保存带有上下调标记的差异分析结果
write.csv(deg,file="upanddown.csv")

#这个是只需要输入GSE号，就能自动获得GLP号
setwd("F:/001微生物课题/公共数据库生信挖掘/03.Hp感染不同时间以后GES-1差异lncRNA/002.GEO数据处理") ##💚💚💚💚💚设置路径

###加载R包
# readxl: 用于读取Excel文件
# tidyverse: 数据处理和可视化的综合工具包
# GEOquery: 用于从GEO数据库下载和处理数据
# limma: 用于差异表达分析
# affy: 用于处理Affymetrix微阵列数据
# stringr: 字符串处理工具包
library(readxl)
library(tidyverse)
library(GEOquery)
library(tidyverse)
library(GEOquery)
library(limma) 
library(affy)
library(stringr)

###下载数据：可在文件夹中或直接读取
# 从GEO数据库下载GSE111763数据集，并保存在当前目录
# AnnotGPL=T表示下载注释信息，getGPL=T表示下载平台信息
gset = getGEO('GSE111763', destdir=".", AnnotGPL = T, getGPL = T)
class(gset)  # 查看数据类型
gset[[1]]    # 查看第一个数据集的详细信息

# 自动获取GPL平台号
gpl_number <- gset[[1]]@annotation
cat("数据集使用的平台号为:", gpl_number, "\n")

#读取平台文件
# 获取平台的注释信息
GPL_data <- getGEO(filename = paste0(gpl_number, ".soft.gz"), AnnotGPL = T)
# 提取平台注释表格数据
GPL_data_11 <- Table(GPL_data)

#读取表达矩阵
# 从数据集中提取表达矩阵
exp <- exprs(gset[[1]])
# 获取探针名称（行名）
probe_name<-rownames(exp)

#转换ID
# 匹配探针ID，确保表达矩阵和注释信息对应
loc<-match(GPL_data_11[,1],probe_name)
# 获取匹配后的表达矩阵
probe_exp<-exp[loc,]
# 从gene_assignment列获取基因信息并处理
raw_geneid <- as.matrix(GPL_data_11[,"gene_assignment"])
# 从gene_assignment列中提取基因名称（通常是第二个"//"之前的部分）
raw_geneid <- sapply(strsplit(raw_geneid, "//"), function(x) {
    if(length(x) >= 2) {
        # 去除空格并获取第一个有效的基因名
        gene <- trimws(x[2])
        if(gene == "---" || gene == "") return(NA)
        return(gene)
    } else {
        return(NA)
    }
})

# 继续后续处理
index <- which(!is.na(raw_geneid))
geneid <- raw_geneid[index]
exp_matrix <- probe_exp[index,]
geneidfactor<-factor(geneid)
gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean))
rownames(gene_exp_matrix)<-levels(geneidfactor)
gene_exp_matrix=na.omit(gene_exp_matrix)

#####读取分组信息#####
# 获取样本表型数据
pdata <- pData(gset[[1]])
# 根据样本名称确定分组信息：Hp感染组和未感染对照组
group_list <- ifelse(str_detect(pdata$title, "noninfected"), "Control",
                     "Hp")
print(group_list)  # 打印查看分组结果
# 将分组信息转换为因子类型，设定水平顺序（对照组在前，处理组在后）
group_list = factor(group_list,
                    levels = c("Control", "Hp"))
print(group_list)  # 打印查看分组结果
# 将分组信息添加到表型数据中
pdata$group = group_list

#####数据标准化处理#####
# 绘制标准化前的箱线图
boxplot(gene_exp_matrix,outline=T, notch=T,col=group_list, las=2)
dev.off()
# 进行数组间标准化处理
gene_exp_matrix_noemal=normalizeBetweenArrays(gene_exp_matrix)
# 绘制标准化后的箱线图
boxplot(gene_exp_matrix_noemal,outline=T, notch=T,col=group_list, las=2)
# 查看数据范围
range(gene_exp_matrix_noemal)
# 转换为数据框格式
gene_exp_matrix_noemal <-as.data.frame(gene_exp_matrix_noemal)
# 去除含NA的行
gene_exp_matrix_noemal=na.omit(gene_exp_matrix_noemal)
# 保存标准化后的表达矩阵
write.csv(gene_exp_matrix_noemal,file = "geo_exp.csv")
dev.off()

#####进行差异分析#####
# 构建设计矩阵
design=model.matrix(~group_list)
# 拟合线性模型
fit=lmFit(gene_exp_matrix_noemal,design)
# 贝叶斯检验，计算经验贝叶斯统计量
fit=eBayes(fit)
# 获取所有基因的差异分析结果
deg=topTable(fit,coef=2,number = Inf)
# 保存所有差异分析结果
write.table(deg, file = "deg_all.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

# 设置差异表达基因的筛选标准
logFC=1           # log2倍变化阈值
adj.P.Val = 0.05  # 调整后的p值阈值

# 标记上调和下调的基因
k1 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC < -logFC)    # 下调基因
k2 = (deg$adj.P.Val < adj.P.Val)&(deg$logFC > logFC)     # 上调基因
# 添加变化方向标记
deg$change = ifelse(k1,"down",ifelse(k2,"up","stable"))
# 统计上调、下调和稳定的基因数量
table(deg$change)
# 保存带有上下调标记的差异分析结果
write.csv(deg,file="upanddown.csv")