GEO数据挖掘
GEO TCGA:
数据下载:
首先在搜索栏搜索相应的癌症
一个课题Series里面有10个样本Samples,制作基因芯片的公司,检测平台,芯片编号等Platforms
点进自己需要的项目之后
点击
上面肿瘤组,下面正常组
分组,建组之后,数据很少,把相应的选中,再点组名放到组里面
根据p值由小到大一次排列,把具有统计学意义的基因筛查出来
P值,可以见到是由小到大的,p值越小,证明一个基因在正常组和肿瘤组间的差异越有统计学意义,t值是两组值t检验后的值
B:经贝叶斯调整后的得到的标准
LogFC:两组表达量间以2为底数的变化倍数(正常组比肿瘤组,再放上log2)
有些gene是空着的,这可能是还没研究清楚或者功能没确定
1.最好下载原始数据,不要下载矩阵数据,因为它有可能没有经过log2处理,探针表达量在几千几万以上,表达量在14以下,是经过log2处理的
2.原始数据里有探针信息没有基因信息,要利用平台的注释文件将探针转化为基因,通过基因表达量的差异筛选出不同组织间的差异基因
从上至下分别是平台文件,矩阵文件和原始文件
←关于探针
通过扫描得到DAT文件(荧光信号图像文件) 通过GCOS(图像处理软件)得到下面的文件,CEL文件通过RMA算法,背景矫正,标准化等处理得到基因表达数据,如果要处理数据,要选择CEL文件,不要选择其他类型的
将下载的文件解压,然后放到cel文件夹了,然后根据网页上的信息分类
然后运行代码,进行表格生成和合并
代码:
RMA法预处理normal样本
setwd("D:\\GEO\\cel\\normal")
library(affyPLM)
library(affy)
Data<-ReadAffy()
sampleNames(Data)
N=length(Data)
#用RMA预处理数据
eset.rma<-rma(Data)
#获取表达数据并输出到表格
normal_exprs<-exprs(eset.rma)
probeid<-rownames(normal_exprs)
normal_exprs<-cbind(probeid,normal_exprs)
write.table(normal_exprs,file="normal.expres.txt",sep='\t',quote=F,row.names=F)
RMA法预处理tumor样本
setwd("D:\\GEO\\cel\\tumor")
library(affyPLM)
library(affy)
Data<-ReadAffy()
sampleNames(Data)
N=length(Data)
#用RMA预处理数据
eset.rma<-rma(Data)
#获取表达数据并输出到表格
normal_exprs<-exprs(eset.rma)
probeid<-rownames(normal_exprs)
normal_exprs<-cbind(probeid,normal_exprs)
write.table(normal_exprs,file="tumor.expres.txt",sep='\t',quote=F,row.names=F)
合并N和T的数据
#setwd(" ")
normal_exprs<-read.table("normal.expres.txt",header=T,sep="\t")
tumor_exprs<-read.table("tumor.expres.txt",header=T,sep="\t")
#讲T和N合并
probe_exprs<-merge(normal_exprs,tumor_exprs,by="probeid")
write.table(probe_exprs,file="cancer.probeid.exprs.txt",sep='\t',quote=F,row.names=F)
芯片质量评估:
第一个是灰度图,合格的芯片灰度分布很均一,如果某一处有白色残块,说明有风险
第二个是权重图,绿色表示权重低,白色灰色权重高,如果权重随机分布,如果权重随机分布,绿色点就会分布的比较均匀
第三个是残差图,红色代表低残差,蓝色代表高残差,白色代表负残差,随机分布也会是均匀的
如果第二第三张图都均匀说明芯片质量过关
第四张,残差图符号图
一个条表示一个样本,在0附近表示质量好,不在就说明质量不好,无法利用,结论不靠谱(即使仍是一条线,比如说5,也不行,也代表质量不好)
相对对数表达(RLE)
一个探针组在某个样品的表达值除以该探针组在所有样品中表达之的中位数后取对数
反映平行实验的一致性
如果RLE仍无法确定可以使用下面的:
相对标准差(NUSE)
一个探针组在某个样品的PM值的标准差除以该探针组在各样品中的PM值标准差的中位数后取对数。反映平行实验的一致性
比RLE更为敏感
接近于纵坐标,然后虽然有两个样本略高,但是总体来说还是处在一个水平线上
RNA降解图,5低3高,rna降解从5端开始,如果说为直线或者低斜率,说明芯片不行,可能是这些样品的rna全部被降解了,斜率很大也不行
关于如何画这些图,按下面方法下载数据之后运行代码
代码:
setwd("D:\\GEO\\GSE79973_RAW")
#自己设定
#调用R包
library(affyPLM)
#载入数据
Data<-ReadAffy()
#对数据集进行回归计算
Pset <- affyPLM::fitPLM(Data)
质量控制:查看灰度图
image(Data[,1])
#根据计算结果,画权重图
image(Pset,type="weights",which=1,main="Weights")
#根据计算结果,画残差图
image(Pset,type="resids",which=1,main="Residuals")
#根据计算结果,画残差符号图
image(Pset,type="sign.resids",which=1,main="Residuals.sign")
质量控制:相对对数表达(RLE)
一个探针组在某个样品的表达值除以该探针组在所有样品中表达之的中位数后取对数
反映平行实验的一致性
#调用R包
library(affyPLM)
library(RColorBrewer)
#对数据集进行回归计算
Pset<-fitPLM (Data)
#载入颜色
colors<-brewer.pal(12,"Set3")
#绘制RLE箱线图
Mbox(Pset,col=colors,main="RLE",las=3)
质量控制:相对标准差(NUSE)
一个探针组在某个样品的PM值的标准差除以该探针组在各样品中的PM值标准差的中位数后取对数。反映平行实验的一致性
比RLE更为敏感。
#调用R包
library(affyPLM)
library(RColorBrewer)
#对数据集进行回归计算
Pset<-fitPLM (Data)
#载入颜色
colors<-brewer.pal(12,"Set3")
#绘制NUSE箱线图
boxplot(Pset,col=colors,main="NUSE",las=3)
质量控制:RNA降解图
原理:RNA降解从5’端开始,因为芯片结果5’端荧光强度要远低于3’端
#调用R包
library(affy)
#获取降解数据
data.deg<-AffyRNAdeg(Data)
#绘制RNA降解图
plotAffyRNAdeg(data.deg,col=colors)
#在左上部位添加图注
legend("topleft",sampleNames(Data),col=colors,lwd=1,inset=0.05,cex=0.2)
基因注释:
寻找差异基因:
热图,绿色低表达,红色高表达
假设文章不自带热图要怎么制作呢?
点进去后这样选择
上面一列代表一个样本,左边一列代表一个基因
这样子,点stack up可以扩大,看这括起来的基因怎么表达的
火山图,一个点代表一个样本
自己制作火山图和热图:
首先就是保证文件夹里有这两个文件,并且要处理GPL570-55999.txt这个文件,使它只有这三列,可以先转成excel处理,再复制到txt里
使用代码,最后生成这三个文件,三种矩阵文件分别是(探针,id,symbol):
代码:
Probe ID转换为Gene symbol(对平台文件进行整理)
setwd("D:\\GEO")
#读取基因表达文件
probe_exp<-read.table("cancer.probeid.exprs.txt",header=T,sep="\t",row.names=1)
#读取探针文件
probeid_geneid<-read.table("GPL570-55999.txt",header=T,sep="\t")
probe_name<-rownames(probe_exp)
#probe进行匹配
loc<-match(probeid_geneid[,1],probe_name)
#确定能匹配上的probe表达值
probe_exp<-probe_exp[loc,]
#每个probeid应对的geneid
raw_geneid<-as.numeric(as.matrix(probeid_geneid[,3]))
#找出有geneid的probeid并建立索引
index<-which(!is.na(raw_geneid))
#提取有geneid的probe
geneid<-raw_geneid[index]
#找到每个geneid的表达值
exp_matrix<-probe_exp[index,]
geneidfactor<-factor(geneid)
#多个探针对应1个基因的情况,取平均值
gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean))
#geneid作为行名
rownames(gene_exp_matrix)<-levels(geneidfactor)
geneid<-rownames(gene_exp_matrix)
gene_exp_matrix2<-cbind(geneid,gene_exp_matrix)
write.table(gene_exp_matrix2,file="Gastric.cancer.geneid.exprs.txt",sep='\t',quote=F,row.names=F)
#将gene id 转换为gene symbol
loc<-match(rownames(gene_exp_matrix),probeid_geneid[,3])
rownames(gene_exp_matrix)=probeid_geneid[loc,2]
genesymbol<-rownames(gene_exp_matrix)
gene_exp_matrix3<-cbind(genesymbol,gene_exp_matrix)
write.table(gene_exp_matrix3,file="Gastric.cancer.genesyb.exprs.txt",sep='\t',quote=F,row.names=F)
在运行代码的时候就会有警告说有NA缺失,使用KNN法把缺失值补充上:
功能分析 (GO、KEGG、GSEA):
GO:圆圈越大,更多的差异基因富集在圆圈上,颜色越深,富集程度越高
KEGG:反应寻找出来的差异基因,上游和下游蛋白直接的相互作用,红星代表筛选出来的
差异基因富集在通路上
共表达,每一个颜色代表一个模块,模块内基因功能相似