#程序前准备

Sys.setenv(LANGUAGE = "en")
options(stringsAsFactors = FALSE)
rm(list=ls())

#调用R包

library(GEOquery)
library(dplyr)
library(tidyverse)

#下载并读取数据

gset = getGEO('GSE149507', destdir=".",getGPL =F)#F表示不下载探针文件

 

 

soft=getGEO(filename ="GSE149507_family.soft/GSE149507_family.soft")
gpl=soft@gpls[["GPL23270"]]@dataTable@table

gset=gset[[1]]
pdata=pData(gset)#获得临床信息
exprSet=exprs(gset)#获得表达矩阵

#把探针文件的ENTREZ_GENE_ID转换为GENE SYMBOL

library(org.Hs.eg.db)

keytypes(org.Hs.eg.db)#查看一下有哪些基因ID

x1=gpl$ENTREZ_GENE_ID
x1=as.character(x1)
x2=AnnotationDbi::select(org.Hs.eg.db, keys=x1, columns=c("ENTREZID", "SYMBOL"), keytype="ENTREZID")
gpl=gpl[,c(1,2)]
gpl$gene=x2$SYMBOL

#将平台文件和矩阵文件合并

exp=as.data.frame(exprSet)
exp.pl=merge(gpl,exp,by.x=1,by.y=0)

#将GENE SYMBOL变为行名

rownames(exp.pl)=exp.pl$gene#发现报错了，因为行名有重复

exp.pl.1=distinct(exp.pl,gene,.keep_all = T)#去除重复

rownames(exp.pl.1)=exp.pl.1$gene#发现报错了，因为有缺失值

exp.pl.2=na.omit(exp.pl.1)

rownames(exp.pl.2)=exp.pl.2$gene#转换成功

exp.pl.3=exp.pl.2[,-c(1:3)]#去除无用信息，得到基因表达矩阵

#查看临床分组信息

View(pdata)

a=c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35)
a2=a+1

exp.4=exp.pl.3[,c(a2,a)]#调整顺序，前18列全是normal，后面全是cancer
group_list=c(rep('control',18),rep('tumor',18))

group_list=factor(group_list)

group_list <- relevel(group_list, ref="control")#强制限定顺序

#构建差异分析的矩阵

design=model.matrix(~ group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(exp.4)

#limma包差异分析
library(limma)
fit=lmFit(exp.4,design)#构建线性拟合模型
fit=eBayes(fit)#eBayes()使用trend=TRUE对标准误差进行经验贝叶斯平滑，计算每个对比中每个基因的moderated t-statistic和log-odds。 
allDiff=topTable(fit,coef=2,adjust='fdr',number=Inf)#topTable()给出一个最有可能在给定对比下差异表达的基因列表。 
#topTable函数的coef参数，coef=2是指design的第2列，即tumour，即把tumour与normal进行对比
write.table(allDiff,file = "allDiff.txt",sep = "\t",col.names = NA)

#有配对信息的差异分析
pairinfo = factor(rep(1:18,2))
design=model.matrix(~ pairinfo+group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit) 
allDiff_pair=topTable(fit,adjust='BH',coef="T",number=Inf,p.value=0.05)
##adjust="BH"表示对p值的校正方法，包括了："none", "BH", "BY" and "holm"。