12-genomicsdataanalysis.Rmd

# 基因组学数据分析

- [主页](http://genomicsclass.github.io/book/)

## microarrays

### 原理

- 生成互补DNA探针
- 标记样品中的DNA单链 可以对不同样品标记不同颜色
- 特异性互补反应
- 测定标记物光信号

### 应用

- 测定基因表达

     - 已知序列
     - 3’端（降解从5’端开始）选取11个片段作为探针
     - 样品对11个片段都是高表达则基因高表达
     
- 寻找SNP
    
     - SNP 单核苷酸多态性 用来探索基因型
     - 合成SNP探针
     - 测定对不同探针的响应判断AA AG GG类型
     
- 寻找转录因子结合位点
    
     - 样品处理为含蛋白与不含蛋白两份 去除蛋白后扩增
     - 探针是基因组感兴趣的片段
     - 瓦片分析可知探针与含转录因子DNA结合位点
     - 总DNA作为对照
     
## NGS

### 原理

- DNA打成50~70片段 一个样品片段上亿
- 加上adaptor 固定在板上后原位扩增成束
- 使用标记过的单核苷酸逐碱基对测光强 同时测序大量片段
- 得到测序结果与强度

### 应用

- 寻找SNP
- RNA-seq 测定RNA表达量
- 寻找结合位点 表达量

## 数据分析应用背景

### DNA 甲基化

- CpG 5’端到3‘端 CG
- C上甲基化 复制时该特性会保留
- 临近CpG位点的基因不会被表达
- CpG 成簇存在 称为CpG islands
- bisulfite treatment 可以用来测定CpG是否被甲基化 通过将未甲基化的CpG中的C改为T
- 测序中测定改变率就可知CpG位点甲基化程度与位置

### CHIP-SEQ

- 蛋白结合后固定，洗掉其余片段，然后洗掉蛋白，对序列片段测序得到结合位点

### RNA 测序

- RNA反转录为cDNA测序
- 只有外显子
- 同一基因多种RNA片段
- 均值与方差有相关性 需要进行log变换后分析

## Bioconductor

- [官方说明](http://bioconductor.org/install/)
- 使用`biocLite()`安装，安装后仍需要`library()`才能使用 

```
source("http://bioconductor.org/biocLite.R")
biocLite()
```

### 数据结构

#### 分析数据

- F 行 S 列 F代表芯片特征数，S代表样本数

#### 表型数据

- S 行 V 列 V代表样本特征，为分类或连续变量
- 如果表型数据解释不清，可以建立一个解释样本特征的labelDescription数据框，通过`phenoData <- new("AnnotatedDataFrame",data=pData, varMetadata=metadata)` 建立AnnotatedDataFrame类型数据

#### 实验描述

- MIAME 类型对象
- 描述实验参数

#### 组装数据

- 将分析数据，表型数据，实验描述组装为一个ExpressionSet类型的对象
- `exampleSet <- ExpressionSet(assayData=exprs,phenoData=phenoData,experimentData=experimentData,annotation="hgu95av2")`
- annotation代表了一组相似实验设计的芯片数据的代号，通过相关代号可以索引到芯片特征信息并将其与其他数据如基因型，染色体位置等连接以便于分析
- 从ExpressionSet里可以按照表型数据提取子集，也就是对 S 截取 V 中特定子集 `exampleSet[ , exampleSet$gender == "Male"]`
- `esApply` 用来针对ExpressionSet应用函数

#### 数据集应用

- `library(Biobase)`
- `library(GEOquery)`
- `geoq <- getGEO("GSE9514")` 从基因表达精选集（GEO）上得到数据表达集
- `names(geoq)` 得到文件名
- `e <- geoq[[1]]` 得到数据集
- `dim(e)` 查看表达集维度 给出样本数与特征值，也就是测定序列数
- `dim(exprs(e))` 与上面等同，给出基因分析数据
- `dim(pData(e))` 给出8个样本的信息，信息头用`names(pData(e))`给出
- `dim(fData(e))` 给出特征与信息头列表
- exprs为特征数×样本数矩阵 pdata为样本数×信息头 fdata为特征数×信息
- experimentData(e) 给出实验信息
- annotation(e) 特征注释
- exptData(se)$MIAME 给出实验相关关键信息
- Y <- log2(exprs(bottomly.eset) + 0.5) 对NGS数据加0.5取2为底的对数（防0）得-1，排除掉0后可得MAplot观察数值分布，一般为均值小差异大，均值大相对稳定
- formula 用来定义公式
- model.matrix 用定义的公式生成矩阵
- rowttests(y[, smallset], group[smallset]) 定义分组，设定模型可进行t-test，用火山图来表示

##### Iranges

- `library(IRanges)` 序列范围
- `ir <- IRanges(start = c(3, 5, 17), end = c(10, 8, 20))` 定义序列
- `IRanges(5, 10)` 表示5到10这6个碱基对，可以shift
- `range(ir)` 表示存在ir中序列的起止范围
- `gaps(ir)` 表示寻找ir中间隔片段
- `disjoin(ir)` 表示将ir中序列碎片化后互不重叠的片段

##### GRanges and GRangesList

- `library(GenomicRanges)` 基因范围
- `gr <- GRanges("chrZ", IRanges(start = c(5, 10), end = c(35, 45)), strand = "+", seqlengths = c(chrZ = 100L))` 定义位于染色体chrZ上几个序列范围，认为这些范围共同定义一个基因
- 可以shift，可以定义长度后trim
- `mcols(gr)$value <- c(-1, 4)` 定义该基因类型中的列并赋值
- `grl <- GRangesList(gr, gr2)` 多个Granges定义一个基因库
- `length(grl)` 给出基因库里基因个数
- `mcols(grl)$value <- c(5, 7)` 定义该基因库类型中的列并赋值

##### findOverlaps 

- gr1 gr2 为两个基因范围对象
- `fo <- findOverlaps(gr1, gr2)` 寻找两个基因重叠序列
- `queryHits(fo)` 与 `subjectHits(fo)` 提取两个基因重叠序号 成对出现
- `gr1[gr1 %over% gr2]` 提取对应序列范围

##### Rle

- `Rle(c(1, 1, 1, 0, 0, -2, -2, -2, rep(-1, 20)))` 表示4组处理，每组各有 3 2 3 20 个重复
- Rle是一种压缩存储实验设计的方式，可以用`as.numeric()`提取原始数据
- `Views(r, start = c(4, 2), end = c(7, 6)` 提取对应实验组

### 数据读取

- microarray 或 NGS 数据由芯片厂商提供，常见读取原始信息的包有affyPLM、affy、 oligo、limma
- 在Bioconductor里，这些原始数据要转为ExpressionSet格式

#### Affymterix CEL files

- `library(affy)`
- `tab <- read.delim("sampleinfo.txt", check.names = FALSE, as.is = TRUE)` 读取样本信息
- `ab <- ReadAffy(phenoData = tab)` 读取样本数据，探针层次
- `ejust <- justRMA(filenames = tab[, 1], phenoData = tab)` 直接读取为基因层数据
- `e <- rma(ab)` 对样本进行背景校正与正则化，从探针层转化为基因层数据

#### 背景干扰

- spikein方法 梯度加入已知浓度的基因片段 阵列上进行shift 类似拉丁方设计
- 可以看到同一基因不同片段大致符合先平后增模式 开始阶段是噪声主导 后面是浓度主导
- 使用类似基因模拟噪声主导 相减后得到去干扰浓度效应 但低值部分会导致方差过大
- 也可以使用统计建模方法模拟背景值与响应 得到还原度更高的信号

#### 正则化

- 基因组数据大多数为0 加标样品变化 正则化是为了还原这一结果
- 分位数正则化
- 局部回归正则化
- 稳方差正则化
- 当重复实验时 直接用分位数正则会掩盖样品差异 可以考虑只对加标基因正则化 然后推广到全局

#### 探索分析作图

##### MA-plot

- x轴为两组基因组的均值，y轴为两组基因组的均值差
- 用来表示两组平行间的差异

##### Volcano plot

- 横坐标为处理间基因表达差异，纵坐标为差异的`-log10(p.value)`
- 一般为火山喷发状，差异越大，p值越小

## 示例：甲基化数据分析

### 读取数据


```r
devtools::install_github("coloncancermeth","genomicsclass")

library(coloncancermeth)
data(coloncancermeth)
```

该数据集为结肠癌病人与对照的DNA甲基化数据集。

### 数据说明


```r
dim(meth)
dim(pd)
length(gr)
```

meth为测序数据，pd为样本信息，gr测序片段信息。


```r
colnames(pd)
table(pd$Status)
X = model.matrix(~pd$Status)
```

查看病患与正常人的分组并构建模型。

```r
chr = as.factor(seqnames(gr))
pos = start(gr)

library(bumphunter)
cl = clusterMaker(chr,pos,maxGap=500)
res = bumphunter(meth,X,chr=chr,pos=pos,cluster=cl,cutoff=0.1,B=0)
```

按染色体生成因子变量，找出基因起始位点，然后利用bumphunter包寻找甲基化数据中某个阈值（0.1）下甲基化基因聚类的后出现的位置，聚类号，聚类相关性等信息寻找问题基因，可从中提取相关信息

```r
cols=ifelse(pd$Status=="normal",1,2)
Index=(res$table[6,7]-3):(res$table[6,8]+3)
matplot(pos[Index],meth[Index,,drop=TRUE],col=cols,pch=1,xlab="genomic location",ylab="Methylation",ylim=c(0,1))

Index=(res$table[6,7]):(res$table[6,8])

test <- meth[Index,,drop=T]
colnames(test) <- pd$bcr_patient_barcode
test1 <- test[,cols==1]
test2 <- test[,cols==2]

test3 <- apply(test2, 2, mean)
apply(matrix, 1, rank)
```

从上面可以得到有差异的甲基化数据所在的基因位置并提取相关样本数据信息。可根据差异作图，得到两组数据甲基化水平差异所在的基因位置。可对差异进行平滑操作，得到位置。这样就可以知道甲基化发生的序列位置与水平差异的信息了。

-----

下面的例子是用人类基因组数据探索潜在的CpG岛。


```r
library(BSgenome.Hsapiens.UCSC.hg19)

Hsapiens[["chr1"]]

# 计算某染色体上潜在位点个数

countPattern('CG',Hsapiens[["chr1"]]) 

# 计算某染色体上特定序列比例 观察与期望出现的比例

CG <- countPattern('CG',Hsapiens[["chr1"]])/length(Hsapiens[["chr1"]])
GC <- countPattern('GC',Hsapiens[["chr1"]])/length(Hsapiens[["chr1"]])

table <- alphabetFrequency(Hsapiens[["chr1"]])
expect <- table['C']%*%table['G']/(length(Hsapiens[["chr1"]]))^2

CG/expect
```