(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)
懒得上传图片了,大家可以用同样的代码自己实现所有的流程