分析基因芯片的数据，提取出差异表达的基因
这次试验的数据来源于文献：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")

到这里差异表达基因的提取和注释就完成了