- 這次我來詳細的講講數(shù)據(jù)預(yù)處理佣蓉，數(shù)據(jù)預(yù)處理這個過程呢，按我的理解就是先把高通量熒光信號轉(zhuǎn)化成基因表達數(shù)據(jù)亲雪，然后把各種各樣的的數(shù)據(jù)(有缺失值勇凭，有異常值，有值很小的义辕，表達水平為負的等等數(shù)據(jù))轉(zhuǎn)化為正常的虾标，完整的，有可比性的數(shù)據(jù)终息。(剛接觸生物信息夺巩，說的可能不對贞让，希望各位老師指正)
我是用R語言來做的周崭，現(xiàn)在網(wǎng)絡(luò)資源這么發(fā)達，一開始我就做個有感情的代碼搬運工喳张，所以運用了來源于醫(yī)學方的代碼续镇，再加上我自己對代碼的理解。
第一步销部，先提取數(shù)據(jù)摸航，把我們需要的矩陣讀入
這是一個有12列，(樣本)54676行的表達譜矩陣(探針id)要特別注意的是舅桩，在GEO下載的數(shù)據(jù)要經(jīng)過處理一下酱虎，刪除一些注釋文件(開頭和結(jié)尾都有)
probe_exp<-read.table("GSE62533_series_matrix.txt",header=T,sep="\t",row.names=1)#讀入基因表達文件
- 這個文件是有三列數(shù)據(jù)，第一列是探針id擂涛，列名是ID,第二列是基因名读串，列名是Gene_symbol，第三列是基因id，列名為ENTREZ_GENE_ID
probeid_geneid<-read.table("GPL570_probe_geneid.txt",header=T,sep="\t")#讀入探針與基因的文件
**如圖所示
第二步恢暖，探針id對應(yīng)基因symbol
probe_name<rownames(probe_exp)#提取probeid
loc<match(probeid_geneid[,1],probe_name)#probeid進行匹配
probe_exp<-probe_exp[loc,]#確定能匹配上的probe表達值
raw_geneid<-as.numeric(as.matrix(probeid_geneid[,3]))#每個probeid對應(yīng)的geneid
index<-which(!is.na(raw_geneid))#找出有g(shù)eneid的probeid并建立索引
geneid<-raw_geneid[index]#提取與geneid匹配的probeid
exp_matrix<-probe_exp[index,]#找到每個geneid的表達值
geneidfactor<-factor(geneid)
gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean))#多個探針對應(yīng)1個基因的情況排监，取平均值
rownames(gene_exp_matrix)<-levels(geneidfactor)#geneid作為行名
gene_exp_matrix2<-cbind(geneid,gene_exp_matrix)
write.table(gene_exp_matrix2,file="geneid_exp.txt",sep="\t",row.names=F)#寫出geneid表達譜矩陣
loc<match(rownames(gene_exp_matrix),probeid_geneid[,3])#geneid表達譜矩陣和geneid匹配
row.names(gene_exp_matrix)<-probeid_geneid[loc,2]
genesymbol<-rownames(gene_exp_matrix)
gene_exp_matrix3<-cbind(genesymbol,gene_exp_matrix#Gene_symbol這列為表達譜的行名，并與表達譜合并
write.table(gene_exp_matrix3,file="genesymbol_exp.txt",sep="\t",row.names=F,quote=F)#寫出genesymbol表達譜矩陣
- 最后就完成了將probeid轉(zhuǎn)化為genesymbol了
第三步杰捂，補充缺失值
library(impute)
gene_exp_matrix<read.table("genesymbol_exp.txt",sep="\t", header=T,row.names=1)
gene_exp_matrix<as.matrix(gene_exp_matrix)
imputed_gene_matrix<impute.knn(gene_exp_matrix,k=10,rowmax=0.5,colmax=0.8,maxp=3000,rng.seed=362436069)#用的是k值臨近法
GeneExp<imputed_gene_matrix$data#寫入表格
genesymbol<rownames(gene_exp_matrix)
GeneExp<cbind(genesymbol,GeneExp)
write.table(GeneExp,file="gene_exprs.txt",sep="\t",quote=F)#讀出經(jīng)過缺失值處理的數(shù)據(jù)
第四步舆床，數(shù)據(jù)標準化(核心)
library(limma)
setwd("D:/qiheyeban")
rt<read.table("gene_exprs.txt",row.names="genesymbol",sep="\t",header=T)
pdf(file="rawbox.pdf")
boxplot(rt,col="purple",xaxt="n",outline=F)#未校正的箱線圖
dev.off()
rt<normalizeBetweenArrays(as.matrix(rt),method="scale")
pdf(file="normalbox.pdf") boxplot(rt, col ="orange",xaxt ="n",outline = F) #校正后的圖
dev.off()
第五步，differential嫁佳，基因篩選
class <-c(rep("con",6),rep( "treat",6)) #找出對照組實驗組
design<-model.matrix(~factor(class))
colnames(design)<-c("con","treat")
fit<-lmFit(rt,design)#算對照組實驗組平均值方差
fit2<-eBayes(fit)#貝葉斯檢驗得到差異結(jié)果
allDiff=toptable(fit2, adjust= "fdr", coef=2, number=200000) #得到所有基因的差異信息
write.table(allDiff,file="allDiff.txt", sep="\t",quote=F )
diffLab<-allDiff[with(allDiff,((logFC>1|logFC<(-1))&adj.P.Val< 0.05)),] #foldchange差異兩倍以上,校正后的P.value<0.05的篩選出來
write.table(diffLab,file="diffexp.txt", sep="\t",quote=F)
diffExpLeve<-rt[rownames(diffLab),]#差異基因表
qvalue=allDiff[rownames(diffLab), ]$adj.P.Val
diffExpQvalue=cbind(qvalue,diffExpLeve1)
write.table(diffExpQvalue,file="diffexplevel.txt",sep="\t",quote=F)
heatmap挨队，熱圖
library(gplots)
hmExp=log10(diffExpLeve1+0.00001)#log值會報錯
hmMat<-as.matrix(hmExp)
pdf(file="heatmap.pdf",height=150,width=30)
par(oma=c(3,3,3,5))
heatmap.2(hmMat,col='greenred',trace="none",cexCo1=1)
dev.off()
image.png
volcano，火山圖
xMax=max(log10(allDiff$adj.P.Val))#x軸范圍為P.value最值
yMax=max(abs(allDiff$logFC))#y軸范圍為foldchange最值
plot(-log10(allDiff$adj.P.Val),allDiff$logFC,xlab="adj.P.Val",ylab="logFC",main="volcano") #把所有的表達值打黑點縱軸為foldchange蒿往，橫軸為P.value
diffSub=subset(allDiff,allDiff$adj.P.Val<0.05&abs(allDiff$logFC)>1)#篩選出差異基因
points(-log10(diffSub$adj.P.Val),diffSub$logFC,pch=20,col="red",cex=0.4)#差異的打紅點
abline(h=0,lty=2,lwd=3)#畫中線
dev.off()
image.png
寫這些東西單純是為了記錄我的學習過程瞒瘸，內(nèi)容述說的可能很粗糙，但是有一個東西我覺得對學的是這個方向的人是很重要的(對于所有學習知識的人都適用的)熄浓，就是在各種時候報錯(遇到困難的時候)情臭，一定要先從自己代碼本身或者是數(shù)據(jù)格式或者是其他找原因，(迎難而上赌蔑！)把問題弄明白俯在，而不是說就這樣放棄了，相信勤奮的人運氣都不會太差
