GEO是当今最大、最全的公共基因数据资源库，包括基因的表达、突变、修饰等信息，涵盖几乎所有的疾病，且单个实验检测样品数目较多，是我们分析、学习的很好资源。
实验设计
原始文章对14个溃疡性结肠炎病人 (Ulcerative colitis, UC)和15个克罗恩病病人 (Crohn’s disease, CD)的发炎组织和未发炎组织活检采样，用Affy芯片检测基因表达谱，研究发炎组织和未发炎组织的基因表达差异。(Genome-wide Pathway Analysis Using Gene Expression Data of Colonic Mucosa in Patients with Inflammatory Bowel Disease. Inflamm Bowel Dis. 杂志影响因子不高，但是领域的专业杂志)
检测并安装依赖的包
a = rownames(installed.packages())
install_bioc <- c("Biobase","oligoClasses","ArrayExpress", "pd.hugene.1.0.st.v1",
"hugene10sttranscriptcluster.db", "oligo", "arrayQualityMetrics",
"limma", "topGO", "ReactomePA", "clusterProfiler", "gplots", "ggplot2",
"geneplotter", "RColorBrewer", "pheatmap", "dplyr","stringr","genefilter")
for(i in install_bioc) {if(! i %in% a) BiocManager::install(i, update=F)}
# General Bioconductor packages
library(Biobase)
library(oligoClasses)
# Annotation and data import packages
library(ArrayExpress)
library(pd.hugene.1.0.st.v1)
library(hugene10sttranscriptcluster.db)
# Quality control and pre-processing packages
library(oligo)
library(arrayQualityMetrics)
# Analysis and statistics packages
library(limma)
library(topGO)
library(ReactomePA)
library(clusterProfiler)
# Plotting and color options packages
library(gplots)
library(ggplot2)
library(geneplotter)
library(RColorBrewer)
library(pheatmap)
# Formatting/documentation packages
#library(rmarkdown)
#library(BiocStyle)
library(dplyr)
#library(tidyr)
# Helpers:
library(stringr)
library(matrixStats)
library(genefilter)
#library(openxlsx)
#library(devtools)
从ArrayExpress下载原始数据
实验检测采用的是Affy芯片 (A-AFFY-141 - Affymetrix GeneChip Human Gene 1.0 ST Array [HuGene-1_0-st-v1])，原始数据为CEL格式，存储于ArrayExpress，索引号是E-MTAB-2967。
image.png
使用getAE函数下载原始数据和注释数据到当前工作目录 (Rmd所在目录)，返回一个包含所有下载文件名的列表。
# type: 'raw' to download and extract only the raw data, 'processed' to download and extract only the processed data or 'full' to have both raw and processed data.
# anno_AE <- getAE("E-MTAB-2967", type = "raw")
若自动下载不成功，手动点击下载图中File部分所有文件，放置到当前目录。(也可去后台回复affy数据 获取)
image.png
# type: 'raw' to download and extract only the raw data, 'processed' to download and extract only the processed data or 'full' to have both raw and processed data.
# 修改参数local=T，读入当前目录的数据
anno_AE <- getAE("E-MTAB-2967", type = "raw", local=T)
## Warning in getAE("E-MTAB-2967", type = "raw", local = T): No processed data
## files found in directory /disk1/train/GEO/f1000
anno_AE里面存储了所有的文件的路径和名字信息。并把原始数据解压缩，释放出CEL文件到当前目录，方便后续读取。
anno_AE
## $path
## [1] "/disk1/train/GEO/f1000"
##
## $rawFiles
## [1] "164_I_.CEL" "164_II.CEL" "183_I.CEL" "183_II.CEL" "2114_I.CEL"
## [6] "2114_II.CEL" "2209_A.CEL" "2209_B.CEL" "2255_I.CEL" "2255_II.CEL"
## [11] "2400_I.CEL" "2400_II.CEL" "2424_A.CEL" "2424_B.CEL" "255_I.CEL"
## [16] "255_II.CEL" "2826_I.CEL" "2826_II.CEL" "2853_I.CEL" "2853_II.CEL"
## [21] "2978_I.CEL" "2978_II.CEL" "2987_I.CEL" "2987_II.CEL" "2992_I.CEL"
## [26] "2992_II.CEL" "2995_I.CEL" "2995_II.CEL" "321_I.CEL" "321_II.CEL"
## [31] "3222_I.CEL" "3222_II.CEL" "3223_I.CEL" "3223_II.CEL" "3226_I.CEL"
## [36] "3226_II.CEL" "3233_I.CEL" "3233_II.CEL" "3258_I.CEL" "3258_II.CEL"
## [41] "3259_I.CEL" "3259_II.CEL" "3262_I.CEL" "3262_II.CEL" "3266_I.CEL"
## [46] "3266_II.CEL" "3269_I.CEL" "3269_II.CEL" "3271_I.CEL" "3271_II.CEL"
## [51] "3302_I.CEL" "3302_II.CEL" "3332_I.CEL" "3332_II.CEL" "848_A.CEL"
## [56] "848_B.CEL" "888_I.CEL" "888_II.CEL"
##
## $rawArchive
## [1] "E-MTAB-2967.raw.1.zip" "E-MTAB-2967.raw.2.zip"
##
## $processedFiles
## NULL
##
## $processedArchive
## character(0)
##
## $sdrf
## [1] "E-MTAB-2967.sdrf.txt"
##
## $idf
## [1] "E-MTAB-2967.idf.txt"
##
## $adf
## [1] "A-AFFY-141.adf.txt"
原始数据解释
ArrayExpress的每个数据根据MAGE-TAB (MicroArray Gene Expression Tabular)指南存储，包含5种不同类型的文件：
IDF (Investigation Description Format)：研究描述文件，包含实验信息如题目、描述、提交者联系方式、实验操作过程等。
ADF (Array Design Format): 芯片设计文件
SDRF (Sample and Data Relationship Format): 样本属性文件，如实验分组、处理方式、取样部位等。
原始数据文件 (raw)
加工后的数据文件 (processed data files)
ExpressionSets数据结构描述
组学数据通常比较复杂，包含很多不同的部分，如实验样品信息、基因组注释信息和实验数据，对应到芯片数据是样品属性和分组信息、不同来源的基因ID信息和功能注释信息、基因表达矩阵。
为了更好的组织这些数据，Bioconductor的Biobase包定义了一个标准化的数据结构 (ExpressionSet类)存储这些数据。ExpressionSet类包含下面几部分：
assayData: 芯片实验的表达数据，探针在行，样品名字在列，用函数
exprs获取metaData
phenoData: 样品描述信息，样品名字在行，实验分组、处理方式、取样部位在列。通常是
SDRF文件的信息的读入。用函数pData获取。featureData: 基因注释特征信息，行为探针名字，列为探针对应的基因、转录本名字和相应的注释信息等。用函数
fData获取。experimentData: 对实验描述的其它信息。
ExpressionSet类可以帮我们协调数据修改、提取过程中的所有信息的一致性。但是要注意phenoData的行名字必须与assayData的列名字一致，都是样品标识符；assayData的行名字必须与featureData的行名字一致，都是基因标识符，具体见下图：
image.png
导入数据，存储为”ExpressionSet”
读入SDRF数据。
sdrf_location <- file.path("E-MTAB-2967.sdrf.txt")
SDRF <- read.delim(sdrf_location)
我们查看下，读进来的数据长什么样子，有哪些信息？
SDRF[1:3,1:5]
总共有哪些样品相关的信息？取样个体标记、物种、疾病、取样部分、是否患病、raw data存储的位置和名字等。
t(SDRF[1, ])
## 1
## Source.Name "164_I"
## Characteristics.individual. "164"
## Characteristics.organism. "Homo sapiens"
## Characteristics.disease. "Crohn's disease"
## Characteristics.organism.part. "colon"
## Characteristics.phenotype. "non-inflamed colonic mucosa"
## Material.Type "organism part"
## Protocol.REF "P-MTAB-41361"
## Protocol.REF.1 "P-MTAB-41363"
## Extract.Name "164_I"
## Protocol.REF.2 "P-MTAB-41364"
## Labeled.Extract.Name "164_I:Biotin"
## Label "biotin "
## Protocol.REF.3 "P-MTAB-41366"
## Assay.Name "164_I_"
## Technology.Type "array assay"
## Array.Design.REF "A-AFFY-141"
## Term.Source.REF "ArrayExpress"
## Protocol.REF.4 "P-MTAB-41367"
## Array.Data.File "164_I_.CEL"
## Comment..ArrayExpress.FTP.file. "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip"
## Factor.Value.disease. "Crohn's disease"
## Factor.Value.phenotype. "non-inflamed colonic mucosa"
Array.Data.File存储了样品名字信息 (CEL文件名)，用作行名字。把SDRF数据表转为AnnotatedDataFrame格式用于下游构建ExpressionSet对象。
rownames(SDRF) <- SDRF$Array.Data.File
SDRF <- AnnotatedDataFrame(SDRF)
Affymetrix芯片的原始数据是CEL文件，里面包含检测到的探针杂交密度值，代表原始基因表达量。另外每个CEL还含有额外信息，如芯片类型、扫描时间等，常用于做批次校正。
image.png
oligo包的read.celfiles函数可以读取这些文件。(虽然不影响，但有几个文件命名不规律，如164_I_多写了个下划线，有几个样品用A,B代替了I,II。不规范的名字是不太利于批量分析和下游识别的，引以为戒。)
raw_data <- oligo::read.celfiles(filenames = as.character(SDRF$Array.Data.File),
verbose = FALSE, phenoData = SDRF)
stopifnot(validObject(raw_data))
获得的raw_data就是一个ExpressionSet对象。raw_data@phenoData@data(或者pData(raw_data))是对应的样品属性信息，与SDRF信息一致。raw_data@annotation获取注释平台信息pd.hugene.1.0.st.v1。
str(raw_data)
## Formal class 'GeneFeatureSet' [package "oligoClasses"] with 9 slots
## ..@ manufacturer : chr "Affymetrix"
## ..@ intensityFile : chr NA
## ..@ assayData :<environment: 0x815aba0>
## ..@ phenoData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 23 obs. of 2 variables:
## .. .. .. ..$ labelDescription: chr [1:23] NA NA NA NA ...
## .. .. .. ..$ channel : Factor w/ 2 levels "exprs","_ALL_": 2 2 2 2 2 2 2 2 2 2 ...
## .. .. ..@ data :'data.frame': 58 obs. of 23 variables:
## .. .. .. ..$ Source.Name : Factor w/ 58 levels "164_I","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Characteristics.individual. : int [1:58] 164 164 183 183 2114 2114 2209 2209 2255 2255 ...
## .. .. .. ..$ Characteristics.organism. : Factor w/ 1 level "Homo sapiens": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.disease. : Factor w/ 2 levels "Crohn's disease",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.organism.part. : Factor w/ 1 level "colon": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.phenotype. : Factor w/ 2 levels "inflamed colonic mucosa",..: 2 1 2 1 2 1 2 1 2 1 ...
## .. .. .. ..$ Material.Type : Factor w/ 1 level "organism part": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF : Factor w/ 1 level "P-MTAB-41361": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.1 : Factor w/ 1 level "P-MTAB-41363": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Extract.Name : Factor w/ 58 levels "164_I","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Protocol.REF.2 : Factor w/ 1 level "P-MTAB-41364": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Labeled.Extract.Name : Factor w/ 58 levels "164_I:Biotin",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Label : Factor w/ 1 level "biotin ": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.3 : Factor w/ 1 level "P-MTAB-41366": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Assay.Name : Factor w/ 58 levels "164_I_","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Technology.Type : Factor w/ 1 level "array assay": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Array.Design.REF : Factor w/ 1 level "A-AFFY-141": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Term.Source.REF : Factor w/ 1 level "ArrayExpress": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.4 : Factor w/ 1 level "P-MTAB-41367": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Array.Data.File : Factor w/ 58 levels "164_I_.CEL","164_II.CEL",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Comment..ArrayExpress.FTP.file.: Factor w/ 2 levels "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Factor.Value.disease. : Factor w/ 2 levels "Crohn's disease",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Factor.Value.phenotype. : Factor w/ 2 levels "inflamed colonic mucosa",..: 2 1 2 1 2 1 2 1 2 1 ...
## .. .. ..@ dimLabels : chr [1:2] "rowNames" "columnNames"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ featureData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 0 obs. of 1 variable:
## .. .. .. ..$ labelDescription: chr(0)
## .. .. ..@ data :'data.frame': 1102500 obs. of 0 variables
## .. .. ..@ dimLabels : chr [1:2] "featureNames" "featureColumns"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ experimentData :Formal class 'MIAME' [package "Biobase"] with 13 slots
## .. .. ..@ name : chr ""
## .. .. ..@ lab : chr ""
## .. .. ..@ contact : chr ""
## .. .. ..@ title : chr ""
## .. .. ..@ abstract : chr ""
## .. .. ..@ url : chr ""
## .. .. ..@ pubMedIds : chr ""
## .. .. ..@ samples : list()
## .. .. ..@ hybridizations : list()
## .. .. ..@ normControls : list()
## .. .. ..@ preprocessing : list()
## .. .. ..@ other : list()
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 2
## .. .. .. .. .. ..$ : int [1:3] 1 0 0
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ annotation : chr "pd.hugene.1.0.st.v1"
## ..@ protocolData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 2 obs. of 2 variables:
## .. .. .. ..$ labelDescription: chr [1:2] "Names of files used in 'exprs'" "Run dates for files used in 'exprs'"
## .. .. .. ..$ channel : Factor w/ 2 levels "exprs","_ALL_": 2 2
## .. .. ..@ data :'data.frame': 58 obs. of 2 variables:
## .. .. .. ..$ exprs: Factor w/ 58 levels "164_I_.CEL","164_II.CEL",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ dates: Factor w/ 58 levels "2010-06-09T11:39:13Z",..: 56 55 9 10 6 8 58 57 5 7 ...
## .. .. ..@ dimLabels : chr [1:2] "rowNames" "columnNames"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. ..@ .Data:List of 4
## .. .. .. ..$ : int [1:3] 3 5 1
## .. .. .. ..$ : int [1:3] 2 42 0
## .. .. .. ..$ : int [1:3] 1 3 0
## .. .. .. ..$ : int [1:3] 1 0 0
# 如果用Rstudio, View(raw_data)显示更好
# View(raw_data)
查看样本属性信息
head(Biobase::pData(raw_data))
## Source.Name Characteristics.individual. Characteristics.organism.
## 164_I_.CEL 164_I 164 Homo sapiens
## 164_II.CEL 164_II 164 Homo sapiens
## 183_I.CEL 183_I 183 Homo sapiens
## 183_II.CEL 183_II 183 Homo sapiens
## 2114_I.CEL 2114_I 2114 Homo sapiens
## 2114_II.CEL 2114_II 2114 Homo sapiens
## Characteristics.disease. Characteristics.organism.part.
## 164_I_.CEL Crohn's disease colon
## 164_II.CEL Crohn's disease colon
## 183_I.CEL Crohn's disease colon
## 183_II.CEL Crohn's disease colon
## 2114_I.CEL Crohn's disease colon
## 2114_II.CEL Crohn's disease colon
## Characteristics.phenotype. Material.Type Protocol.REF
## 164_I_.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 164_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## 183_I.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 183_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## 2114_I.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 2114_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## Protocol.REF.1 Extract.Name Protocol.REF.2 Labeled.Extract.Name
## 164_I_.CEL P-MTAB-41363 164_I P-MTAB-41364 164_I:Biotin
## 164_II.CEL P-MTAB-41363 164_II P-MTAB-41364 164_II:Biotin
## 183_I.CEL P-MTAB-41363 183_I P-MTAB-41364 183_I:Biotin
## 183_II.CEL P-MTAB-41363 183_II P-MTAB-41364 183_II:Biotin
## 2114_I.CEL P-MTAB-41363 2114_I P-MTAB-41364 2114_I:Biotin
## 2114_II.CEL P-MTAB-41363 2114_II P-MTAB-41364 2114_II:Biotin
## Label Protocol.REF.3 Assay.Name Technology.Type Array.Design.REF
## 164_I_.CEL biotin P-MTAB-41366 164_I_ array assay A-AFFY-141
## 164_II.CEL biotin P-MTAB-41366 164_II array assay A-AFFY-141
## 183_I.CEL biotin P-MTAB-41366 183_I array assay A-AFFY-141
## 183_II.CEL biotin P-MTAB-41366 183_II array assay A-AFFY-141
## 2114_I.CEL biotin P-MTAB-41366 2114_I array assay A-AFFY-141
## 2114_II.CEL biotin P-MTAB-41366 2114_II array assay A-AFFY-141
## Term.Source.REF Protocol.REF.4 Array.Data.File
## 164_I_.CEL ArrayExpress P-MTAB-41367 164_I_.CEL
## 164_II.CEL ArrayExpress P-MTAB-41367 164_II.CEL
## 183_I.CEL ArrayExpress P-MTAB-41367 183_I.CEL
## 183_II.CEL ArrayExpress P-MTAB-41367 183_II.CEL
## 2114_I.CEL ArrayExpress P-MTAB-41367 2114_I.CEL
## 2114_II.CEL ArrayExpress P-MTAB-41367 2114_II.CEL
## Comment..ArrayExpress.FTP.file.
## 164_I_.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 164_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 183_I.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 183_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 2114_I.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 2114_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## Factor.Value.disease. Factor.Value.phenotype.
## 164_I_.CEL Crohn's disease non-inflamed colonic mucosa
## 164_II.CEL Crohn's disease inflamed colonic mucosa
## 183_I.CEL Crohn's disease non-inflamed colonic mucosa
## 183_II.CEL Crohn's disease inflamed colonic mucosa
## 2114_I.CEL Crohn's disease non-inflamed colonic mucosa
## 2114_II.CEL Crohn's disease inflamed colonic mucosa
这么展示更利于查看
t(Biobase::pData(raw_data)[1,])
## 164_I_.CEL
## Source.Name "164_I"
## Characteristics.individual. "164"
## Characteristics.organism. "Homo sapiens"
## Characteristics.disease. "Crohn's disease"
## Characteristics.organism.part. "colon"
## Characteristics.phenotype. "non-inflamed colonic mucosa"
## Material.Type "organism part"
## Protocol.REF "P-MTAB-41361"
## Protocol.REF.1 "P-MTAB-41363"
## Extract.Name "164_I"
## Protocol.REF.2 "P-MTAB-41364"
## Labeled.Extract.Name "164_I:Biotin"
## Label "biotin "
## Protocol.REF.3 "P-MTAB-41366"
## Assay.Name "164_I_"
## Technology.Type "array assay"
## Array.Design.REF "A-AFFY-141"
## Term.Source.REF "ArrayExpress"
## Protocol.REF.4 "P-MTAB-41367"
## Array.Data.File "164_I_.CEL"
## Comment..ArrayExpress.FTP.file. "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip"
## Factor.Value.disease. "Crohn's disease"
## Factor.Value.phenotype. "non-inflamed colonic mucosa"
筛选并保留关心的样本属性信息，个体信息 (Source.Name,Characteristics.individual.)，疾病信息 (Factor.Value.disease.), 发炎与否 (Factor.Value.phenotype.)。
Biobase::pData(raw_data) <- Biobase::pData(raw_data)[, c("Source.Name", "Characteristics.individual.", "Factor.Value.disease.", "Factor.Value.phenotype.")]
原始数据质控
exprs(raw_data)可获取原始的表达信息，行代表芯片探针在芯片上的位置，列代表每个样品。
Biobase::exprs(raw_data)[1:5, 1:5]
## 164_I_.CEL 164_II.CEL 183_I.CEL 183_II.CEL 2114_I.CEL
## 1 4496 5310 4492 4511 2872
## 2 181 280 137 101 91
## 3 4556 5104 4379 4608 2972
## 4 167 217 99 79 82
## 5 89 110 69 58 47
对表达数据取对数，并进行主成分分析。每个点代表一个样品，颜色代表是否发炎，形状代表疾病类型。从图中可以看出，在第一主成分两种疾病区分比较明显，在第二主成分疾病的区别也略大于是否发炎。而我们的关注点是发炎，而不是不同类型的疾病，需要在下游分析中考虑移除疾病的影响。
exp_raw <- log2(Biobase::exprs(raw_data))
PCA_raw <- prcomp(t(exp_raw), scale. = FALSE)
percentVar <- round(100*PCA_raw$sdev^2/sum(PCA_raw$sdev^2),1)
sd_ratio <- sqrt(percentVar[2] / percentVar[1])
dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
Disease = pData(raw_data)$Factor.Value.disease.,
Phenotype = pData(raw_data)$Factor.Value.phenotype.,
Individual = pData(raw_data)$Characteristics.individual.)
ggplot(dataGG, aes(PC1, PC2)) +
geom_point(aes(shape = Disease, colour = Phenotype)) +
ggtitle("PCA plot of the log-transformed raw expression data") +
xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) +
ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) +
theme(plot.title = element_text(hjust = 0.5))+
coord_fixed(ratio = sd_ratio) +
scale_shape_manual(values = c(4,15)) +
scale_color_manual(values = c("darkorange2", "dodgerblue4"))
image
另外，需要看下原始数据的表达分布。从图中看出，不同的芯片原始表达分布不均一，需要做合适的标准化处理。
# oligo::boxplot(raw_data, target = "core",
# main = "Boxplot of log2-intensitites for the raw data")
dataBoxplot <- data.frame(t(exp_raw), Sample=colnames(exp_raw),
Disease = pData(raw_data)$Factor.Value.disease.,
Phenotype = pData(raw_data)$Factor.Value.phenotype.,
Individual = pData(raw_data)$Characteristics.individual.)
dataBoxplot <- reshape2::melt(dataBoxplot, id.vars=c("Sample","Disease","Phenotype","Individual"))
ggplot(dataBoxplot, aes(x=Sample, y=value)) +
geom_boxplot(aes(color=Disease)) + theme_bw() +
theme(axis.text.x=element_text(angle=45, hjust=1, vjust=1))
image.gif
arrayQualityMetrics可以对芯片进行更多层面的检测，并生成一个检测报告。
arrayQualityMetrics(expressionset = raw_data,
outdir = "ysx",
force = TRUE, do.logtransform = TRUE,
intgroup = c("Factor.Value.disease.", "Factor.Value.phenotype."))
背景校正
鉴于探针的荧光强度会受到非特异性杂交和光学检测系统噪音的影响，因此需要对信号强度进行背景校正，从而能更真实的反映和比较特异性杂交产生的信号，更好的反应基因表达情况。
芯片间标准化 (calibration)
不同的芯片杂交受很多因素影响，比如反转录效率，标记或杂交反应，芯片物理性质，试剂批次和实验室环境等，校正后，才可以在不同芯片之间可比。
Summarization
在Affymetrix平台，一个转录本通常用多个探针检测。对每个基因来说，需要根据所有探针的信息估算出其相对于总的转录本的表达值，也就是Summarization。随后就可以使用注释数据集获得基因的symbols或ENSEMBL ID。我们这个平台的注释信息存储在hugene10sttranscriptcluster.db包中。
head(ls("package:hugene10sttranscriptcluster.db"))
## [1] "hugene10sttranscriptcluster"
## [2] "hugene10sttranscriptcluster_dbconn"
## [3] "hugene10sttranscriptcluster_dbfile"
## [4] "hugene10sttranscriptcluster_dbInfo"
## [5] "hugene10sttranscriptcluster_dbschema"
## [6] "hugene10sttranscriptcluster.db"
Affymetrix芯片”probesets”的变更
在采用芯片检测样品前，需要先提取RNA，加上一个荧光标签，并进行至少一轮PCR扩增以增加检测量。传统上这一步是根据mRNA``3'端的poly-A序列设计引物完成反转。这样最终的扩增产物里面3’端的序列会更多些。而且如果起始RNA片段化了或降解了，5’端被检测到的几率更低了。因此设计探针时，更多的探针也是针对3’端设置的，这就是3' IVT Affymetrix芯片。它是基于探针集的，探针集的一组固定探针用于检测一个特定的基因或转录本。在一些情况下，如选择性Poly-A位点出现时，一个基因需要设计多个探针集检测。
而后来的Gene/ExonAffymetrix芯片是基于外显子的，探针集来源于外显子区。这时为了获得基因的表达，需要做两步summarization，第一步是从探针集获得外显子的表达，然后根据多个来源于同一个基因或转录本的外显子获取基因的表达。因此需要合适的基因注释包，如我们这用到的hugene10sttranscriptcluster.db。
Gene芯片相比于Exon芯片探针更少，只保留Exon芯片的一部分”good”探针。在Exon芯片上至少需要4个探针来代表一个Exon,而在Gene芯片上，很多代表一个基因的探针集只有3个或更少的探针。
如下图所示，用单一颜色代表一个探针集，一个基因的探针集用一个颜色集合表示，如所有黄色探针代表一个外显子，所有黄色、橙色、红色探针集都用于检测同一个基因。
image.png
左图中，每个外显子或探针集包含很多探针 (同样颜色的块出现多次)，因此做探针集或外显子水平的summarization是合适的。而对于gene芯片，每个探针集的探针数目很少，因此在探针集或外显子水平的summarization是不推荐的，虽然也可以通过包hugene10stprobeset.db来完成。
而且在Gene和Exon芯片上不再有错配探针 (mismatch probes)。错配探针本来是用做背景校正排除非特异性杂交和背景噪音的，但被更好的计算方式取代。
一步完成背景校正、标准化、表达整合
oligo包提供了一个函数rma允许一步完成基于去卷积的背景校正、分位数校正标准化 (quantile normalization)和RMA (robust multichip average) 表达整合 (summarization)。
相对log表达 (RLE) 质控分析
我们先利用rma函数进行背景校正和summarization获得基因的表达量，而略过标准化。rma函数输出的结果已经做过对数处理。
palmieri_eset <- oligo::rma(raw_data, target = "core", normalize = FALSE)
## Background correcting
## Calculating Expression
然后进行RLE (Relative Log Expression)分析，首先计算每个转录本在所有样品中的表达中位数，然后每个转录本的表达减去这个中位数，整理格式进行绘图。
row_medians_assayData <-
Biobase::rowMedians(as.matrix(Biobase::exprs(palmieri_eset)))
# 减去中值
RLE_data <- sweep(Biobase::exprs(palmieri_eset), 1, row_medians_assayData)
RLE_data <- as.data.frame(RLE_data)
# 转换为长格式
RLE_data_gathered <-
tidyr::gather(RLE_data, patient_array, log2_expression_deviation)
ggplot2::ggplot(RLE_data_gathered, aes(patient_array, log2_expression_deviation)) +
geom_boxplot(outlier.shape = NA) +
ylim(c(-2, 2)) +
theme(axis.text.x = element_text(colour = "aquamarine4",
angle = 60, size = 6.5, hjust = 1 ,
face = "bold"))
image.gif
图中，Y-轴代表每个芯片中转录本表达与该转录本在所有芯片表达中位数的偏差。箱体越延展表明芯片中转录本的表达量与中位表达差别越大，箱体在Y-轴的偏移表示大部分转录本的表达定量存在系统的增高或降低。这可能存在质量问题或是检测批次效应造成的。因此如果一个样品的boxplot的形状和中值与其它样品差别很大，需要进一步评估是否为异常样本，并考虑移除。
图中5个中位数在-0.7~-0.8之间的样品 2826_I, 2826_II, 3262_II,3302_II和3332_II可能是异常样品。需要结合后面的聚类分析，再确认是否移除。
校正三部曲
background-correct, normalize and summarize.
palmieri_eset_norm <- oligo::rma(raw_data, target = "core")
## Background correcting
## Normalizing
## Calculating Expression
参数target定义summarization的程度，默认是core表示获得基因水平的表达 (using transcript clusters containing “safely” annotated genes) (对gene芯片，只有这么一个选项)。如果是外显子芯片，还可以选择extended和full。如果想获得外显子水平的表达，也可以用probeset做参数，但是一般不推荐。
除了RMA还有其它的背景校正和标准化方法，但RMA通常是一个好的默认选择。RMA在所有芯片间使用quantile normalization标准化方法使得不同样品检测的表达值可比。
image
出校准化后的表达值
norm_exprs <- Biobase::exprs(palmieri_eset_norm)
norm_exprs <- data.frame(ID=rownames(norm_exprs), norm_exprs)
write.table(norm_exprs, "normalized_probe_expr.tsv", sep="\t", row.names=F, quote=F)
校准化后数据的质量评估
主成分分析
exp_palmieri <- Biobase::exprs(palmieri_eset_norm)
PCA <- prcomp(t(exp_palmieri), scale = FALSE)
percentVar <- round(100*PCA$sdev^2/sum(PCA$sdev^2),1)
sd_ratio <- sqrt(percentVar[2] / percentVar[1])
dataGG <- data.frame(PC1 = PCA$x[,1], PC2 = PCA$x[,2],
Disease = Biobase::pData(palmieri_eset_norm)$Factor.Value.disease.,
Phenotype = Biobase::pData(palmieri_eset_norm)$Factor.Value.phenotype.)
ggplot(dataGG, aes(PC1, PC2)) +
geom_point(aes(shape = Disease, colour = Phenotype)) +
ggtitle("PCA plot of the calibrated, summarized data") +
xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) +
ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) +
theme(plot.title = element_text(hjust = 0.5)) +
coord_fixed(ratio = sd_ratio) +
scale_shape_manual(values = c(4,15)) +
scale_color_manual(values = c("darkorange2", "dodgerblue4"))
image.gif
第一主坐标轴大体可以分开发炎与未发炎样品，第二主坐标轴分开两种不同的疾病。(从图上看，还是疾病的区分更明显)
热图聚类分析
通过热图展示样品之间的相似度关系，并且标记样品的属性如是否发炎或是否来自同一个疾病。
# 生成行注释
phenotype_names <- ifelse(str_detect(pData(palmieri_eset_norm)$Factor.Value.phenotype., "non"), "non_infl.", "infl.")
disease_names <- ifelse(str_detect(pData(palmieri_eset_norm)$Factor.Value.disease., "Crohn"), "CD", "UC")
annotation_for_heatmap <- data.frame(Phenotype = phenotype_names, Disease = disease_names)
row.names(annotation_for_heatmap) <- row.names(pData(palmieri_eset_norm))
计算Manhattan距离，并进行聚类分析
# 使用dist计算样本之间的距离，首选做下转置，因为dist是对列进行操作的
# dist默认是欧式距离，反应直线路径的平方距离
# 这里用Manhattan距离，计算的是直角路径的绝对距离，计算结果更稳定
dists <- as.matrix(dist(t(exp_palmieri), method = "manhattan"))
rownames(dists) <- row.names(pData(palmieri_eset_norm))
hmcol <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "YlOrRd"))(255))
colnames(dists) <- NULL
# 设置对角距离为NA，增强颜色的展示力
diag(dists) <- NA
# www.ehbio.com/ImageGP
ann_colors <- list(
Phenotype = c(non_infl. = "chartreuse4", infl. = "burlywood3"),
Disease = c(CD = "blue4", UC = "cadetblue2"))
pheatmap(dists, col = (hmcol),
annotation_row = annotation_for_heatmap,
annotation_colors = ann_colors,
legend = TRUE,
treeheight_row = 0,
legend_breaks = c(min(dists, na.rm = TRUE),
max(dists, na.rm = TRUE)),
legend_labels = (c("small distance", "large distance")),
main = "Clustering heatmap for the calibrated samples")
image.gif
热图中上部黄色的条带可能是异常样品 (2826_II, 3262_II, 3271_I, 2978_II 和 3332_II)，部分样品与之前RLE图推测可能是异常的样品重合，可以考虑移除。这里为了跟原文一致，所有样品都保留了下来。
另外arrayQualityMetrics有更精确的计算异常样品的方式，我们后续也会提供基于网络图距离的异常值剔除方法。
过滤低表达基因
芯片数据中有一大部分探针的信号值与背景值相差不大，而且在不同样品中差别不大，不能给我们的生物问题提供更多信息。(注：虽然有部分表达低的基因表达变化比较小时就可以有比较大的作用，但因为其表达低，检测噪音也大，结果可信度一般)
首选绘制下每个转录本在所有芯片中的表达中位数的分布 (共计33,297转录本)。可以看到左侧区域富集低表达的转录本，具体选择移除的阈值主观性比较强，这里选择表达中位数需要大于4，给后续分析多保留一些基因。也可以排序一下，保留前75%的基因。
palmieri_medians <- rowMedians(Biobase::exprs(palmieri_eset_norm))
hist_res <- hist(palmieri_medians, 100, col = "cornsilk1", freq = FALSE,
main = "Histogram of the median intensities",
border = "antiquewhite4",
xlab = "Median intensities")
image
man_threshold <- 4
hist_res <- hist(palmieri_medians, 100, col = "cornsilk", freq = FALSE,
main = "Histogram of the median intensities",
border = "antiquewhite4",
xlab = "Median intensities")
abline(v = man_threshold, col = "coral4", lwd = 2)
image
假如最小的实验组有5个样品，如果一个转录本检测到表达值大于之前设定的阈值 (man_threshold)的样本数小于5，则移除。
首先获得各个实验组的样品数
no_of_samples <-
table(paste0(pData(palmieri_eset_norm)$Factor.Value.disease., "_",
pData(palmieri_eset_norm)$Factor.Value.phenotype.))
no_of_samples
##
## Crohn's disease_inflamed colonic mucosa
## 15
## Crohn's disease_non-inflamed colonic mucosa
## 15
## ulcerative colitis_inflamed colonic mucosa
## 14
## ulcerative colitis_non-inflamed colonic mucosa
## 14
# 获得最小样品数
samples_cutoff <- min(no_of_samples)
# sum(x > man_threshold)：转录本在多少样品中表达值高于给定阈值
idx_man_threshold <- apply(Biobase::exprs(palmieri_eset_norm), 1,
function(x){sum(x > man_threshold) >= samples_cutoff})
table(idx_man_threshold)
## idx_man_threshold
## FALSE TRUE
## 10493 22804
# 提取过滤转录本
palmieri_manfiltered <- subset(palmieri_eset_norm, idx_man_threshold)
注释转录本簇
在我们构建线性模型做差异分析之前，先把探针映射到Gene symbol，并只保留这部分基因用于后续分析。
# 使用select函数，获得探针ID对应的Gene symbol和描述
anno_palmieri <- AnnotationDbi::select(hugene10sttranscriptcluster.db,
keys = (featureNames(palmieri_manfiltered)),
columns = c("SYMBOL", "GENENAME"),
keytype = "PROBEID")
# 过滤掉没有Symbol的探针
anno_palmieri <- subset(anno_palmieri, !is.na(SYMBOL))
head(anno_palmieri)
如果一个探针对应于多个Gene symbol，我们没有办法判断到底是哪个，只能移除。
# 按探针ID分组
anno_grouped <- group_by(anno_palmieri, PROBEID)
# 计算每组Symbol的个数
anno_summarized <-
dplyr::summarize(anno_grouped, no_of_matches = n_distinct(SYMBOL))
head(anno_summarized)
# 过滤Symbol个数大于1的
anno_filtered <- filter(anno_summarized, no_of_matches > 1)
head(anno_filtered)
probe_stats <- anno_filtered
nrow(probe_stats)
## [1] 1763
有1763个探针ID对应多个Gene Symbol，从Assya数据中移除这些探针ID。
# featureNames(palmieri_manfiltered): 获取所有探针名字
ids_to_exlude <- (featureNames(palmieri_manfiltered) %in% probe_stats$PROBEID)
table(ids_to_exlude)
## ids_to_exlude
## FALSE TRUE
## 21041 1763
palmieri_final <- subset(palmieri_manfiltered, !ids_to_exlude)
validObject(palmieri_final)
## [1] TRUE
前面是从assay数据中移除了这些探针ID，后面需要从feature data中移除:
head(anno_palmieri)
现在的feature data是空的，增加一列PROBID。
# 增加一列PROBID
fData(palmieri_final)$PROBEID <- rownames(fData(palmieri_final))
# 使用left_join根据PROBID列，合并，增加注释信息
fData(palmieri_final) <- left_join(fData(palmieri_final), anno_palmieri)
## Joining, by = "PROBEID"
# restore rownames after left_join
rownames(fData(palmieri_final)) <- fData(palmieri_final)$PROBEID
validObject(palmieri_final)
## [1] TRUE
