• 这次我来详细的讲讲数据预处理，数据预处理这个过程呢，按我的理解就是先把高通量荧光信号转化成基因表达数据，然后把各种各样的的数据(有缺失值，有异常值，有值很小的，表达水平为负的等等数据)转化为正常的，完整的，有可比性的数据。(刚接触生物信息，说的可能不对，希望各位老师指正)
  我是用R语言来做的，现在网络资源这么发达，一开始我就做个有感情的代码搬运工，所以运用了来源于医学方的代码，再加上我自己对代码的理解。

第一步，先提取数据，把我们需要的矩阵读入

这是一个有12列，(样本)54676行的表达谱矩阵(探针id)要特别注意的是，在GEO下载的数据要经过处理一下，删除一些注释文件(开头和结尾都有)

probe_exp<-read.table("GSE62533_series_matrix.txt",header=T,sep="\t",row.names=1)#读入基因表达文件

• 这个文件是有三列数据，第一列是探针id，列名是ID,第二列是基因名，列名是Gene_symbol，第三列是基因id，列名为ENTREZ_GENE_ID

probeid_geneid<-read.table("GPL570_probe_geneid.txt",header=T,sep="\t")#读入探针与基因的文件

**如图所示

第二步，探针id对应基因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对应的geneid
index<-which(!is.na(raw_geneid))#找出有geneid的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))#多个探针对应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转化为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)#读出经过缺失值处理的数据

第四步，数据标准化(核心)

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)#贝叶斯检验得到差异结果
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

写这些东西单纯是为了记录我的学习过程，内容述说的可能很粗糙，但是有一个东西我觉得对学的是这个方向的人是很重要的(对于所有学习知识的人都适用的)，就是在各种时候报错(遇到困难的时候)，一定要先从自己代码本身或者是数据格式或者是其他找原因，(迎难而上！)把问题弄明白，而不是说就这样放弃了，相信勤奋的人运气都不会太差