阅读 24 SEO

01_2_差异表达分析实操

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

作者:生信小白

原文链接:https://www.jianshu.com/p/dd9f42cecf2a

文章分类
后端
文章标签
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐