本文代码在mac环境下运行
常用的表达谱分析软件：
1.ArrayTools ；
2.DNA-Chip Analyzer ；
3.SAM ；
4.R（“affy”，“Limma”，“marray”）；
5.Matlab（Bioinformatics Toolbox）
安装Bioconductor Packages：
source("http://bioconductor.org/biocLite.R")
bioLite() ##安装Bioconductor
安装Limma包
bioLite("limma") ##安装limma包
bioLite(c("GenomicFeatures","AnnotationDbi")) ##同时安装两个包
数据预处理：
不同波长的荧光特异性涉及到前景值和背景值，因此要进行背景矫正；之后剔除不可信的点且对确实值进行估计；之后将数据归一化（样本之间归一化使得样本之间可以对比，基因之间标准化使得基因之间可以对比）以后有利于后续下游分析
数据标准化是为了消除系统偏差引起的高相关性，保留真正的生物学引起的基因表达的高相关性
数据预处理（Affy包，limma包）：
1.Affymetrix芯片都属于单通道芯片(Oligonucleotide Arrays)，使用affy包
library("Affy")
ReadAffy(); ##输入数据
expresson(); ##背景矫正；标准化；归一化
justRMA(); ##more efficient
exprs();
library(simpleaffy)
ampli.eset <- call.exprs(cel,"mas5",sc=target)
qcs <- qc(cel.ampli.eset)
2.双通道芯片(Two-Color Spotted Arrays)，可以使用limma包
library(limma)
read.maimages(); ##输入数据
backgroundCorrect(); ##背景矫正
normalizeWithinArrays(); ##芯片标准化
normalizeBetweenArrays(): ##芯片之间标准化
exprs.MA(); ##抽取表达值
avereps(); ##归一化
plotMA(); ##绘图
差异表达基因筛选（主流）
library(limma)
lmFit(); ##芯片数据的线性模型
eBayes(); ##贝叶斯统计方法
##差异表达基因
实操代码（修改版）
##设置工作路径
setwd("/Users/apple/Desktop/生信学习材料/生信入门之路/代码数据挖掘_淘宝/实用数据挖掘/学习记录 /01_差异基因分析/实操/")
##安装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")
##加载R包
library("Biobase")
library("affy")
library("limma")
##import phenotype data 样本信息读取
phenoData = read.AnnotatedDataFrame('target.txt')
pheno = pData(phenoData)
View(pheno)
#import Annotaion 注释信息读取
anno = read.csv("Annotation.csv",head=T)
View(head(anno))
##RMA normalization用RMA读取原始数据
#eset.rma = RMA(Data)
eset.rma <- justRMA(filenames=paste(rownames(pheno),'.CEL',sep=''), celfile.path='./GSE12251') ##文件路径
datExpr = exprs(eset.rma)
#安装impute包计算缺失值
library(impute)
#KNN法计算缺失值
imputed_gene_exp = impute.knn(datExpr,k=10,rowmax = 0.5, colmax=0.8, maxp=3000, rng.seed=362436069)
datExpr2 = imputed_gene_exp$data
##得到表达矩阵输出
write.table(datExpr2,file="Expdata.txt",sep="\t")
#######
#样本分组
Group = factor(pheno$group,levels=c('Responder','Nonresponder'))
design = model.matrix(~0+Group)
colnames(design) <- c('Responder','Nonresponder')
design
#线性模型拟合
fit <- lmFit(datExpr2, design)
#构建比对模型，比较两个条件下的表达数据
contrast.matrix <- makeContrasts(Responder-Nonresponder, levels=design) ##正值代表R对比NR高表达,负代表低表达
#########################################
library(openxlsx)
library(futile.logger)
#比对模型进行差值计算
fit2 <- contrasts.fit(fit, contrast.matrix)
#贝叶斯检验
fit2 <- eBayes(fit2)
#找出差异基因检验结果，并且输出符合条件的结果
diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05, lfc=log(2,2),number=5000,sort.by = 'logFC') ##P值和LogFC调整阈值
#diff ID转换，基因注释
diff$Gene = anno$Gene.Symbol[match(rownames(diff),anno$ID_REF)]
diff$ID_REF = rownames(diff)
diff = diff[,c(8,7,1:6)]
diff = diff[diff$Gene != '---',]
#output 输出差异表达基因
write.xlsx(diff,'DEG.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=F,append=T)
记录2021.4.17
网友评论