(19)一个affymetrix表达芯片实战-生信菜鸟团博客2周年精选文章集

这个实例上部分包括:

如何用R包下载GEO数据(只限单一平台,其余平台需要修改下面的代码)

如何对GEO的芯片数据归一化并且得到表达量矩阵,

如何用limma包做差异分析,

对找到的差异基因如何做GO和KEGG注释

首先下载两个GEO数据:

平台是:Affymetrix U133 gene chips

67 diseased triple negative breast cancer samples(GSE31519 )and 42 control samples (GSE20437)

都是表达量数据,同一种芯片。分成两组,正好做差异表达。

数据来源的文献是:

文章title:A clinically relevant gene signature in triple negative and basal-like breast cancer

结论(We describe a ratio of high B-cell presence and low IL-8 activity as a powerful new prognostic marker for TNBC. )

地址:http://www.breast-cancer-research.com/content/13/5/R97

GEO数据地址: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31519

Platform: GPL96 67 Samples

Download data: GEO (CEL, TXT)

SeriesAccession: GSE31519ID: 200031519

文章title:Histologically normal epithelium from breast cancer patients and cancer-free prophylactic mastectomy patient

GEO数据地址:  http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20437

Platform: GPL96 67 Samples

Download data: GEO (CEL, TXT)

SeriesAccession: GSE31519  ID: 200031519

Platform: GPL96 Series: GSE20437   42 Samples

Download data: GEO (CEL)

DataSetAccession: GDS3716  ID: 3716

我首先用R的GEOquery包来下载。(其实你完全可以直接去GEO网站下载数据,然后解压的)

suppressMessages(library(GEOquery))

setwd("D:\\test_analysis\\TNBC")

gse31519=getGEO("GSE31519",GSEMatrix = T,destdir = "./")

getGEOSuppFiles("GSE31519",baseDir = "./")

gse31519=getGEO("GSE20437",GSEMatrix = T,destdir = "./")

getGEOSuppFiles("GSE20437",baseDir = "./")

这样下载之后的数据都存在D:\\test_analysis\\TNBC下面

接下来我们就用affy包和limma来进行差异分析:

library(affy)

library(limma)

affy.data=ReadAffy(celfile.path="./cel_files")

请先搞清楚,ReadAffy 这个函数的用法!

当前工作目录下面有没有cel_files文件夹?

如果你搞不清楚什么是工作变量,什么是函数,什么是参数,如何赋值,请不要再看下去了,我懒得跟你说!

cel_files文件夹下面有没有文件?

eset.rma=rma(affy.data)

exprSet=exprs(eset.rma)

write.table(exprSet,"expr_rma_matrix.txt",quote=F,sep="\t")

group=factor(c(rep("control",42),rep("case",67)))

design = model.matrix(~0+group)

colnames(design)=c("case","control")

rownames(design)=sampleNames(affy.data)

fit=lmFit(exprSet,design)

cont.matrix = makeContrasts(contrasts="case-control",levels=design)

fit2=contrasts.fit(fit,cont.matrix)

fit2=eBayes(fit2)

diff_dat=topTable(fit2,coef=1,n=Inf)

write.table(diff_dat,"diff_dat.txt",quote=F)

这样得到的diff_dat就是我们差异分析的结果啦

we choose the log fold cut off change to be “2” to get a manageable set of genes.
原文说:we were able to get a list of 2567 genes after removing the duplicates and the not available genes

我们仅仅根据一个标准来挑选差异基因, the log fold cut off change to be “2”,我只挑出来了782个探针

接下来对这些探针进行注释,得到基因名,我这里用biomart包来进行注释

我们的平台是:Affymetrix U133 gene chips,虽然有22283个探针,但是只有13908个基因

所以代码如下:

ensembl=useMart("ensembl",dataset="hsapiens_gene_ensembl")

gene_probe=getBM(attributes=c("hgnc_symbol","affy_hg_u133a"),filter="affy_hg_u133a",value=rownames(diff_dat),mart=ensembl)

diff_probe=rownames(diff_dat[abs(diff_dat[,1])>2,])

diff_gene=gene_probe[match(diff_probe,gene_probe[,2]),1]

diff_gene=na.omit(diff_gene)

diff_gene=unique(diff_gene)

length(diff_gene)

这样会得到604个差异基因

然后我做一下GO和KEGG的富集分析!

gene_entrez=getBM(attributes=c("hgnc_symbol","entrezgene"),filter="hgnc_symbol",value=diff_gene,mart=ensembl)

## 必须安利一下Y叔的包,希望Y叔看到了会表扬我!

## 如果需要知道Y叔这个包的详细用法,请直接联系他本人!

require(DOSE)

require(clusterProfiler)

gene_entrez=na.omit(gene_entrez)

gene=as.character(gene_entrez[,2])

ego <- enrichGO(gene=gene,organism="human",ont="CC",pvalueCutoff=0.01,readable=TRUE)

ekk <- enrichKEGG(gene=gene,organism="human",pvalueCutoff=0.01,readable=TRUE)

write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)

write.csv(summary(ego),GO-enrich.csv,row.names =F)

懒得上传图片了,大家可以用同样的代码自己实现所有的流程

(0)

相关推荐