R语言_Affymetrix芯片数据处理

采用数据集GSE66360

##设置工作路径(事先放好原始数据，分组信息，注释文件)

setwd("/Users/apple/Desktop/生信学习材料/生信入门之路/代码数据挖掘/实用数据挖掘/学习记录 /01_差异基因分析/实操模仿66360/")

##安装affy和limma包

if (!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install(version = "3.12")

BiocManager::install("affy")

BiocManager::install("limma")

BiocManager::install("impute")

install.packages("openxlsx")

install.packages("dplyr")

##加载R包

library("Biobase")

library("affy")

library("limma")

library(openxlsx)

library(futile.logger)

library(dplyr)

##import phenotype data 样本信息读取(手动整理分组信息)

phenoData = read.AnnotatedDataFrame('target.txt')

pheno = pData(phenoData)

View(pheno)

##import Annotaion 注释信息读取（注释探针矩阵）

#仅留下ID列和gene symbol列

anno = read.csv("Annotation.csv",head=T)

View(head(anno))

##Affy包，RMA函数用来标准化。

文章来源：https://mp.weixin.qq.com/s/hvKxGm9nxWr5CThAZ4PK_g

#01_eset.rma = RMA(Data) 设置文件路径，读入原始CEL文件

eset.rma <- justRMA(filenames=paste(rownames(pheno),'.CEL',sep=''), celfile.path='./GSE66360_RAW/')

#获得矫正以后具体的表达值，提取探针水平表达矩阵

datExpr = exprs(eset.rma)

#02_安装impute包计算缺失值（使用KNN法计算缺失值）

#来源生信技能树帖：https://mp.weixin.qq.com/s/6bR21_LvOelwZN5qen_CRw

#impute参数默认的k = 10, 选择K个邻居的值平均或者加权后填充

默认的rowmax = 0.5, 就是说该行的缺失值比例超过50%就使用平均值而不是K个邻居

默认的colmax = 0.8，意思是该列缺失值超过80%就报错

library(impute)

imputed_gene_exp = impute.knn(datExpr,k=10,rowmax = 0.5,

                              colmax=0.8,maxp =3000, rng.seed=362436069) 

datExpr2 = imputed_gene_exp$data

View(datExpr2)     #datExpr2为计算缺失值的矩阵

#输出未均值处理的表达矩阵datExpr2，手动添加ID列名

write.table(datExpr2,file="Unproccesed_datExpr2.txt",sep="\t")

#读入手动修改的表达矩阵datExpr3

datExpr3 <- read.csv('./proccesed_datExpr2.csv',header = T)

View(datExpr3)

#将(未处理的表达矩阵datExpr3)和(注释探针矩阵anno)合并为一个矩阵datExpr4

datExpr4 <- merge(anno,datExpr3,by='ID')

View(datExpr4)

#03_多个探针对应一个基因的情况取平均值得到datExpr5表达矩阵

datExpr5 <- aggregate(x = datExpr4[,3:ncol(datExpr4)],

                      by = list(datExpr4$Gene.Symbol), #根据gene名

                      FUN = mean) #重复基因表达量取平均值

View(datExpr5)

##04导出最终表达矩阵Final_Expdata

write.table(datExpr5,file="Final_Expdata.txt",sep="\t")  #输出txt文件

write.csv(datExpr5, file ="Final_Expdata.csv",sep ="", row.names =TRUE, col.names =TRUE, quote =TRUE)  #导出csv文件

#######差异表达分析

#样本分组（样本表达矩阵）

Group = factor(pheno$group,levels=c('MI','control')) #样本进行分组

design = model.matrix(~0+Group) #对样本的实验设计进行计算

colnames(design) <- c('MI','control')

design

###Limma包差异分析三步：

#手动删除datExpr5第一列数字列，得到datExpr6，操作使得基因名作为行名

datExpr6 <- read.csv('./Final_Expdata.csv',header = T)  

rownames(datExpr6)=datExpr6[,1]    #取出第一列，基因名作为行名，为后续lmfit()函数构建

datExpr6=datExpr6[,-1]          #将第一列删除

head(datExpr6)

##01_第一步：线性模型拟合

fit <- lmFit(datExpr6, design) #构建比对模型，比较两个条件下的表达数据。注意基因名作为行名，否则可能会报错。文章来源：https://shengxin.ren/question/1446

contrast.matrix <- makeContrasts(MI-control,levels=design)  #正值代表R对比NR高表达,负代表低表达（通常是实验组减对照组）

fit2 <- contrasts.fit(fit, contrast.matrix)  #比对模型进行差值计算

##02_第二步：贝叶斯检验

fit2 <- eBayes(fit2)

##03_第三步：找出差异基因检验结果，并且输出符合条件的结果;

#以下选择一种差异基因筛选模式运行即可

##01_All_Degs输出

diff <-topTable(fit2, coef=1, n=Inf) %>% na.omit()   

#02_P值和LogFC调整阈值输出

diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05,lfc=log(2,2),number=5000,sort.by = 'logFC')   

View(diff)

dim(diff)

#output 输出差异表达基因

write.xlsx(diff,'DEG_MI vs control.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=T,append=T)

记录：2021.4.19