分析基因芯片的數(shù)據(jù)，提取出差異表達的基因
這次試驗的數(shù)據(jù)來源于文獻：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20986
文獻的目的是比較 非傳代，增殖的HUVEC（人臍靜脈血管內(nèi)皮細胞）和iris（虹膜），retina（視網(wǎng)膜），choroid(脈絡膜)的基因表達譜，研究眼部疾病的病理生理機制，看看能否用HUVEC 細胞替代其他細胞作為研究眼科疾病的樣本
1.下載芯片數(shù)據(jù)并解壓
#設置工作路徑
> setwd("C:\\Users\\18019\\Desktop\\xinpian")
#安裝GEOquery包并下載原始數(shù)據(jù)
>source("http://www.bioconductor.org/biocLite.R")
>biocLite("GEOquery")
> library(GEOquery)
> getGEOSuppFiles("GSE20986")#下載數(shù)據(jù)
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讀取數(shù)據(jù)并標準化
使用read.affy函數(shù)讀取cel文件需要兩個輸入
第一個輸入是一個文本文檔，能夠被讀取成數(shù)據(jù)框，第一列無列名，為cel文件名稱，第二列為樣品標簽，可以在文獻上面找到
第二個為文件所在的路徑
第一個輸入我們直接在data文件夾里編輯文本，命名為phenodata.txt。列與列之間用分隔符隔開 如圖
然后就可以載入數(shù)據(jù)并標準化
>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#查看一下標準化后的數(shù)據(jù)，這里可以明顯看到數(shù)據(jù)集已經(jīng)是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.數(shù)據(jù)芯片可視化
在進行下一步的差異分析前，我們要對標準化后的數(shù)據(jù)進行質(zhì)量檢查
#做箱線圖對比標準前后數(shù)據(jù)
#設置顏色
> library(RColorBrewer)
> cols<-brewer.pal(8,"Set1")
#做箱線圖
> boxplot(celfiles,col=cols)
> library(affyPLM)#做標準化后的圖要載入affyPLM包，以便解析標準化后的數(shù)據(jù)
載入需要的程輯包：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.確認完數(shù)據(jù)質(zhì)量后，就可以開始過濾探針了
探針與目標之間一定存在著非特異性結(jié)合，所以所有的探針均會產(chǎn)生信號。
如果不加以過濾，認為這些探針對應的基因都表達，即不符合事實，也會對后續(xù)的分析產(chǎn)生影響。因此，過濾掉表達量低的探針是十分必要的。
rnsFilter() 函數(shù)提供一站式的探針過濾，能夠一步對Expression Set的探針進行過濾，返回一個list。list中的eset為過濾后的Expression Set, filter.log為每一步過濾到多少探針的記錄。nsFilster也可以用于oligo包讀取的對象
常用參數(shù)：
（1）require.entrez 若為TRUE，過濾掉沒有ENTREZ ID的探針。默認為TRUE，一般設為FALSE
（2）remove.dupEntrez 若為TRUE，當幾個探針對應統(tǒng)同一ENTREZ ID的時候，留下 var.func 值最大的探針，其余過濾。默認為TRUE,一般設為FALSE
（3）var.func 用于過濾的統(tǒng)計參數(shù)。默認為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 個探針位點因為是內(nèi)參而被過濾掉。
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.構(gòu)建差異分析的比對矩陣
根據(jù)文獻意思，我們將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.開始差異分析
#將線性模型擬合到過濾后的表達數(shù)據(jù)集上
> fit<-lmFit(celfiles.filtered$eset,design)
#將比對矩陣與線性模型擬合，比較不同細胞系表達數(shù)據(jù)
> huvec_fit<-contrasts.fit(fit,contrasts.matrix)
# 使用經(jīng)驗貝葉斯算法計算差異表達基因的顯著性
> 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")
#將差異基因與值是結(jié)合并輸出
> 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")
到這里差異表達基因的提取和注釋就完成了
