分析基因芯片的数据，提取出差异表达的基因
这次试验的数据来源于文献：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20986
文献的目的是比较 非传代，增殖的HUVEC（人脐静脉血管内皮细胞）和iris（虹膜），retina（视网膜），choroid(脉络膜)的基因表达谱，研究眼部疾病的病理生理机制，看看能否用HUVEC 细胞替代其他细胞作为研究眼科疾病的样本
1.下载芯片数据并解压
#设置工作路径
> setwd("C:\\Users\\18019\\Desktop\\xinpian")
#安装GEOquery包并下载原始数据
>source("http://www.bioconductor.org/biocLite.R")
>biocLite("GEOquery")
> library(GEOquery)
> getGEOSuppFiles("GSE20986")#下载数据
trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE20nnn/GSE20986/suppl//GSE20986_RAW.tar?tool=geoquery'
Content type 'application/x-tar' length 56360960 bytes (53.8 MB)
> list.files()#显示文件
[1]"GSE20986"
#解压文件到data文件夹
> setwd("GSE20986")
> list.files()
[1] "data" "GSE20986_RAW.tar"
> untar("GSE20986_RAW.tar",exdir = "data")
> setwd("data")
> list.files()
[1] "GSM524662.CEL.gz" "GSM524663.CEL.gz" "GSM524664.CEL.gz"
[4] "GSM524665.CEL.gz" "GSM524666.CEL.gz" "GSM524667.CEL.gz"
[7] "GSM524668.CEL.gz" "GSM524669.CEL.gz" "GSM524670.CEL.gz"
[10] "GSM524671.CEL.gz" "GSM524672.CEL.gz" "GSM524673.CEL.gz"
> cels<-list.files()
> sapply(cels, gunzip)
GSM524662.CEL.gz GSM524663.CEL.gz GSM524664.CEL.gz GSM524665.CEL.gz
13555726 13555055 13555639 13560122
GSM524666.CEL.gz GSM524667.CEL.gz GSM524668.CEL.gz GSM524669.CEL.gz
13555663 13557614 13556090 13560054
GSM524670.CEL.gz GSM524671.CEL.gz GSM524672.CEL.gz GSM524673.CEL.gz
13555971 13554926 13555042 13555290
> list.files()
[1] "GSM524662.CEL" "GSM524663.CEL" "GSM524664.CEL" "GSM524665.CEL"
[5] "GSM524666.CEL" "GSM524667.CEL" "GSM524668.CEL" "GSM524669.CEL"
[9] "GSM524670.CEL" "GSM524671.CEL" "GSM524672.CEL" "GSM524673.CEL"
2.用simpleaffy读取数据并标准化
使用read.affy函数读取cel文件需要两个输入
第一个输入是一个文本文档，能够被读取成数据框，第一列无列名，为cel文件名称，第二列为样品标签，可以在文献上面找到
第二个为文件所在的路径
第一个输入我们直接在data文件夹里编辑文本，命名为phenodata.txt。列与列之间用分隔符隔开 如图
然后就可以载入数据并标准化
>library(simpleaffy)
> celfiles<-read.affy(covdesc = "phenodata.txt",path=".")
> celfiles#查看一下是否读取成功
AffyBatch object
size of arrays=1164x1164 features (22 kb)
cdf=HG-U133_Plus_2 (54675 affyids)
number of samples=12
number of genes=54675
annotation=hgu133plus2
notes=
Warning messages:
1: replacing previous import ‘AnnotationDbi::tail’ by ‘utils::tail’ when loading ‘hgu133plus2cdf’
2: replacing previous import ‘AnnotationDbi::head’ by ‘utils::head’ when loading ‘hgu133plus2cdf’
#用GC-RMA算法进行标准化
> celfiles.gcrma<-gcrma(celfiles)
Adjusting for optical effect............Done.
Computing affinitiesLoading required package: AnnotationDbi
Loading required package: stats4
Loading required package: IRanges
Loading required package: S4Vectors
Attaching package: ‘S4Vectors’
The following object is masked from ‘package:base’:
expand.grid
Attaching package: ‘IRanges’
The following object is masked from ‘package:simpleaffy’:
members
The following object is masked from ‘package:grDevices’:
windows
.Done.
Adjusting for non-specific binding............Done.
Normalizing
Calculating Expression
> celfiles.gcrma#查看一下标准化后的数据，这里可以明显看到数据集已经是ExpressionESet了，不是未标准化的AffyBatch
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 12 samples
element names: exprs
protocolData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL
(12 total)
varLabels: ScanDate
varMetadata: labelDescription
phenoData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL
(12 total)
varLabels: sample Target
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu133plus2
3.数据芯片可视化
在进行下一步的差异分析前，我们要对标准化后的数据进行质量检查
#做箱线图对比标准前后数据
#设置颜色
> library(RColorBrewer)
> cols<-brewer.pal(8,"Set1")
#做箱线图
> boxplot(celfiles,col=cols)
> library(affyPLM)#做标准化后的图要载入affyPLM包，以便解析标准化后的数据
载入需要的程辑包：gcrma
载入需要的程辑包：preprocessCore
> boxplot(celfiles.gcrma,col=cols)
如图 标准化之前的箱线图
标准化之后的箱线图
#做密度曲线图
> boxplot(celfiles.gcrma,col=cols)
> hist(celfiles,col=cols)
> hist(celfiles.gcrma,col=cols)
如图，标准化之前的密度曲线图
标准化之后的密度曲线图
#对单个芯片进行可视化
> celfiles.qc<-fitPLM(celfiles)
> image(celfiles.qc,which=1,add.lenge=TRUE)
> image(celfiles.qc,which=6,add.lenge=TRUE)
4.确认完数据质量后，就可以开始过滤探针了
探针与目标之间一定存在着非特异性结合，所以所有的探针均会产生信号。
如果不加以过滤，认为这些探针对应的基因都表达，即不符合事实，也会对后续的分析产生影响。因此，过滤掉表达量低的探针是十分必要的。
rnsFilter() 函数提供一站式的探针过滤，能够一步对Expression Set的探针进行过滤，返回一个list。list中的eset为过滤后的Expression Set, filter.log为每一步过滤到多少探针的记录。nsFilster也可以用于oligo包读取的对象
常用参数：
（1）require.entrez 若为TRUE，过滤掉没有ENTREZ ID的探针。默认为TRUE，一般设为FALSE
（2）remove.dupEntrez 若为TRUE，当几个探针对应统同一ENTREZ ID的时候，留下 var.func 值最大的探针，其余过滤。默认为TRUE,一般设为FALSE
（3）var.func 用于过滤的统计参数。默认为IQR
（4）var.cutoff 截断值。默认为0.5，即过滤到50%的基因。
> celfiles.filtered<-nsFilter(celfiles.gcrma,require.entrez = FALSE,remove.dupEntrez = FALSE)
> celfiles.filtered$filter.log
$`numLowVar`
[1] 27307
$feature.exclude
[1] 62
看出有 27307 个探针位点因为无明显表达差异（LowVar）被过滤掉，有 62 个探针位点因为是内参而被过滤掉。
5.用limma包开始差异分析
1.提取样品信息
>library(limma)
>sample<-celfiles.gcrma$Target
> design<-model.matrix(~0+sample)
>colnames(design) <- c("choroid", "huvec", "iris", "retina")
> design
choroid huvec iris retina
1 0 0 1 0
2 0 0 0 1
3 0 0 0 1
4 0 0 1 0
5 0 0 0 1
6 0 0 1 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
10 0 1 0 0
11 0 1 0 0
12 0 1 0 0
2.构建差异分析的比对矩阵
根据文献意思，我们将HUVEC分别与iris，retrina,choroid分别进行比对
> contrasts.matrix<-makeContrasts(huvec_choroid = huvec - choroid, huvec_retina = huvec - retina, huvec_iris <- huvec - iris, levels=design)
> contrasts.matrix
Contrasts
Levels huvec_choroid huvec_retina huvec_iris <- huvec - iris
choroid -1 0 0
huvec 1 1 1
iris 0 0 -1
retina 0 -1 0
3.开始差异分析
#将线性模型拟合到过滤后的表达数据集上
> fit<-lmFit(celfiles.filtered$eset,design)
#将比对矩阵与线性模型拟合，比较不同细胞系表达数据
> huvec_fit<-contrasts.fit(fit,contrasts.matrix)
# 使用经验贝叶斯算法计算差异表达基因的显著性
> huvec_ebFit <- eBayes(huvec_fit)
# coef=1 是 huvec_choroid 比对矩阵, coef=2 是 huvec_retina 比对矩阵，依此类推
> topTable(huvec_ebFit,number = 10,coef = 1)
logFC AveExpr t P.Value adj.P.Val B
204779_s_at 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
227377_at 3.672710 3.209739 36.08942 4.132286e-12 1.880603e-08 17.03075
204882_at -5.353547 6.513690 -34.70646 6.141288e-12 2.395629e-08 16.77396
38149_at -5.054267 6.483784 -31.45817 1.661577e-11 5.282891e-08 16.09358
205576_at 6.586393 4.139624 31.31294 1.741230e-11 5.282891e-08 16.06036
205453_at 3.624587 3.210563 30.74481 2.095602e-11 5.722251e-08 15.92787
> huvec_choroid<-topTable(huvec_ebFit,coef = 1)
> huvec_retina<-topTable(huvec_ebFit,coef = 2)
> huvec_iris<-topTable(huvec_ebFit,coef = 3)
6.获取差异基因的注释并输出
#安装探针注释的包，可以在文献中找到
> biocLite("hgu133plus2.db")
> library(hgu133plus2.db)
#载入注功能的包
> library(annotate)
#分别对三组差异表达的探针进行注释
> gene.symbols1 <- getSYMBOL(rownames(huvec_choroid), "hgu133plus2")
> gene.symbols2 <- getSYMBOL(rownames(huvec_retina), "hgu133plus2")
> gene.symbols3 <- getSYMBOL(rownames(huvec_iris), "hgu133plus2")
#将差异基因与值是结合并输出
> results<-cbind(gene.symbols1,huvec_choroid)
> head(results)
gene.symbols1 logFC AveExpr t P.Value adj.P.Val B
204779_s_at HOXB7 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at <NA> 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at <NA> 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at IL1RL1 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at NLGN1 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
227377_at IGF2BP1 3.672710 3.209739 36.08942 4.132286e-12 1.880603e-08 17.03075
> write.csv(results,file="huvec_choroid.csv")
> results<-cbind(gene.symbols2,huvec_retina)
> write.csv(results,file="huvec_retina.csv")
> results<-cbind(gene.symbols3,huvec_iris)
> write.csv(results,file = "huvec_iris.csv")
到这里差异表达基因的提取和注释就完成了
网友评论