R語言_Affymetrix芯片數(shù)據(jù)處理
采用數(shù)據(jù)集GSE66360
##設(shè)置工作路徑(事先放好原始數(shù)據(jù)，分組信息，注釋文件)
setwd("/Users/apple/Desktop/生信學(xué)習(xí)材料/生信入門之路/代碼數(shù)據(jù)挖掘/實(shí)用數(shù)據(jù)挖掘/學(xué)習(xí)記錄 /01_差異基因分析/實(shí)操模仿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 樣本信息讀取(手動(dòng)整理分組信息)
phenoData = read.AnnotatedDataFrame('target.txt')
pheno = pData(phenoData)
View(pheno)
##import Annotaion 注釋信息讀?。ㄗ⑨屘结樉仃嚕?/h1>#僅留下ID列和gene symbol列
anno = read.csv("Annotation.csv",head=T)
View(head(anno))
##Affy包，RMA函數(shù)用來標(biāo)準(zhǔn)化。
文章來源：https://mp.weixin.qq.com/s/hvKxGm9nxWr5CThAZ4PK_g
#01_eset.rma = RMA(Data) 設(shè)置文件路徑，讀入原始CEL文件
eset.rma <- justRMA(filenames=paste(rownames(pheno),'.CEL',sep=''), celfile.path='./GSE66360_RAW/')
#獲得矯正以后具體的表達(dá)值，提取探針?biāo)奖磉_(dá)矩陣
datExpr = exprs(eset.rma)
#02_安裝impute包計(jì)算缺失值（使用KNN法計(jì)算缺失值）
#來源生信技能樹帖：https://mp.weixin.qq.com/s/6bR21_LvOelwZN5qen_CRw
#impute參數(shù)默認(rèn)的k = 10, 選擇K個(gè)鄰居的值平均或者加權(quán)后填充
默認(rèn)的rowmax = 0.5, 就是說該行的缺失值比例超過50%就使用平均值而不是K個(gè)鄰居
默認(rèn)的colmax = 0.8，意思是該列缺失值超過80%就報(bào)錯(cuò)
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為計(jì)算缺失值的矩陣
#輸出未均值處理的表達(dá)矩陣datExpr2，手動(dòng)添加ID列名
write.table(datExpr2,file="Unproccesed_datExpr2.txt",sep="\t")
#讀入手動(dòng)修改的表達(dá)矩陣datExpr3
datExpr3 <- read.csv('./proccesed_datExpr2.csv',header = T)
View(datExpr3)
#將(未處理的表達(dá)矩陣datExpr3)和(注釋探針矩陣anno)合并為一個(gè)矩陣datExpr4
datExpr4 <- merge(anno,datExpr3,by='ID')
View(datExpr4)
#03_多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況取平均值得到datExpr5表達(dá)矩陣
datExpr5 <- aggregate(x = datExpr4[,3:ncol(datExpr4)],
? ? ? ? ? ? ? ? ? ? ? by = list(datExpr4$Gene.Symbol), #根據(jù)gene名
? ? ? ? ? ? ? ? ? ? ? FUN = mean) #重復(fù)基因表達(dá)量取平均值
View(datExpr5)
##04導(dǎo)出最終表達(dá)矩陣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) ?#導(dǎo)出csv文件
#######差異表達(dá)分析
#樣本分組（樣本表達(dá)矩陣）
Group = factor(pheno$group,levels=c('MI','control')) #樣本進(jìn)行分組
design = model.matrix(~0+Group) #對(duì)樣本的實(shí)驗(yàn)設(shè)計(jì)進(jìn)行計(jì)算
colnames(design) <- c('MI','control')
design
###Limma包差異分析三步：
#手動(dòng)刪除datExpr5第一列數(shù)字列，得到datExpr6，操作使得基因名作為行名
datExpr6 <- read.csv('./Final_Expdata.csv',header = T) ?
rownames(datExpr6)=datExpr6[,1] ? ?#取出第一列，基因名作為行名，為后續(xù)lmfit()函數(shù)構(gòu)建
datExpr6=datExpr6[,-1]? ? ? ? ? #將第一列刪除
head(datExpr6)
##01_第一步：線性模型擬合
fit <- lmFit(datExpr6, design) #構(gòu)建比對(duì)模型，比較兩個(gè)條件下的表達(dá)數(shù)據(jù)。注意基因名作為行名，否則可能會(huì)報(bào)錯(cuò)。文章來源：https://shengxin.ren/question/1446
contrast.matrix <- makeContrasts(MI-control,levels=design)? #正值代表R對(duì)比NR高表達(dá),負(fù)代表低表達(dá)（通常是實(shí)驗(yàn)組減對(duì)照組）
fit2 <- contrasts.fit(fit, contrast.matrix)? #比對(duì)模型進(jìn)行差值計(jì)算
##02_第二步：貝葉斯檢驗(yàn)
fit2 <- eBayes(fit2)
##03_第三步：找出差異基因檢驗(yàn)結(jié)果，并且輸出符合條件的結(jié)果;
#以下選擇一種差異基因篩選模式運(yùn)行即可
##01_All_Degs輸出
diff <-topTable(fit2, coef=1, n=Inf) %>% na.omit() ??
#02_P值和LogFC調(diào)整閾值輸出
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 輸出差異表達(dá)基因
write.xlsx(diff,'DEG_MI vs control.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=T,append=T)
記錄：2021.4.19
