本文代碼在mac環(huán)境下運(yùn)行
常用的表達(dá)譜分析軟件：
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")) ##同時安裝兩個包
數(shù)據(jù)預(yù)處理：
不同波長的熒光特異性涉及到前景值和背景值，因此要進(jìn)行背景矯正；之后剔除不可信的點(diǎn)且對確實(shí)值進(jìn)行估計；之后將數(shù)據(jù)歸一化（樣本之間歸一化使得樣本之間可以對比，基因之間標(biāo)準(zhǔn)化使得基因之間可以對比）以后有利于后續(xù)下游分析
數(shù)據(jù)標(biāo)準(zhǔn)化是為了消除系統(tǒng)偏差引起的高相關(guān)性，保留真正的生物學(xué)引起的基因表達(dá)的高相關(guān)性
數(shù)據(jù)預(yù)處理（Affy包，limma包）：
1.Affymetrix芯片都屬于單通道芯片(Oligonucleotide Arrays)，使用affy包
library("Affy")
ReadAffy(); ?##輸入數(shù)據(jù)
expresson(); ##背景矯正；標(biāo)準(zhǔn)化；歸一化
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(); ##輸入數(shù)據(jù)
backgroundCorrect(); ##背景矯正
normalizeWithinArrays(); ##芯片標(biāo)準(zhǔn)化
normalizeBetweenArrays(): ##芯片之間標(biāo)準(zhǔn)化
exprs.MA(); ##抽取表達(dá)值
avereps(); ##歸一化
plotMA(); ##繪圖
差異表達(dá)基因篩選（主流）
library(limma)
lmFit(); ##芯片數(shù)據(jù)的線性模型
eBayes(); ##貝葉斯統(tǒng)計方法
? ? ? ? ? ? ? ? ?##差異表達(dá)基因
實(shí)操代碼（修改版）
##設(shè)置工作路徑
setwd("/Users/apple/Desktop/生信學(xué)習(xí)材料/生信入門之路/代碼數(shù)據(jù)挖掘_淘寶/實(shí)用數(shù)據(jù)挖掘/學(xué)習(xí)記錄 /01_差異基因分析/實(shí)操/")?
##安裝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讀取原始數(shù)據(jù)?
#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?
##得到表達(dá)矩陣輸出?
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)?
#構(gòu)建比對模型，比較兩個條件下的表達(dá)數(shù)據(jù)?
contrast.matrix <- makeContrasts(Responder-Nonresponder, levels=design) ?##正值代表R對比NR高表達(dá),負(fù)代表低表達(dá)?
#########################################?
library(openxlsx)
library(futile.logger)?
#比對模型進(jìn)行差值計算?
fit2 <- contrasts.fit(fit, contrast.matrix)
#貝葉斯檢驗(yàn)?
fit2 <- eBayes(fit2)?
#找出差異基因檢驗(yàn)結(jié)果，并且輸出符合條件的結(jié)果?
diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05, lfc=log(2,2),number=5000,sort.by = 'logFC') ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ????????????????????? ? ? ? ? ? ? ? ? ? ? ? ? ##P值和LogFC調(diào)整閾值?
#diff ID轉(zhuǎn)換，基因注釋?
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 輸出差異表達(dá)基因?
write.xlsx(diff,'DEG.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=F,append=T)?
記錄2021.4.17
