14-shengxin.Rmd

# 生物信息 {#bios}

## 数据结构

- 列代表特征 行代表条目
- 每个条目有一个唯一性特征
- 数据表可通过列链接成为关系数据库

## Pubmed 搜索

- PubMed search tags
    - [AD] – Affiliation (company or school)
    - [ALL] – All fields (eliminates defaults)
    - [AU] or [AUTH] – Author
    - [1AU] – First author
    - [ECNO] – Enzyme Commission Numbers
    - [EDAT] – Entry date (YYYY/MM/DD)
    - [ISS] - Issue # of journal
    - [JOUR] - Journal (Title, Abbreviation , ISSN)
    - [LA] – Language
    - [PDAT] – Publication date (YYYY/MM/DD)
    - [PT] – Publication type
    - [SUBS] – Substance name
    - [TIAB] – Title/Abstract
    - [TW] – Text words
    - [UID] – Unique identifiers (primary keys)
    - [VOL] or [VI] – Volume of journal

- MeSH terms [MH][MAJR][SH]
    - 被 MeSH 索引的关系数据库
    - 保守性检索 有层级关系
- 时间段搜索 冒号分割 YYYY/MM/DD:YYYY/MM/DD
- 序列长度搜索 [SLEN] 可以是蛋白 可以是核酸
- 蛋白分子量搜索 [MOLWT]
- 物种搜索 [ORGN]
- Nucleotide 序列蛋白数据库
- [MMDB](http://www.ncbi.nlm.nih.gov/structure/) 3D结构数据库
- [Genome](http://www.ncbi.nlm.nih.gov/genome/) 基因组数据库
- [OMIM](http://omim.org/) 人类孟德尔遗传数据库 用来探索等位基因问题
- [分类数据库](http://www.ncbi.nlm.nih.gov/taxonomy) 用来界定分类
- [GEO](http://www.ncbi.nlm.nih.gov/geo/) 基因芯片的实验数据
- [SNP](http://www.ncbi.nlm.nih.gov/snp/) 基因指纹数据库

## 动态规划

- 用于序列比对
- 对角线得分 按总分评价比对结果
- 可全局 可局部
- 序列比对指标是特异性与相似性
- 特异性指精确匹配比率
- 相似性指精确匹配加化学相似性比率 结构相近则相似
- FASTA 慢准 BLAST 快
- 三种情况 匹配 不匹配 间隔
- 间隔罚分

## 得分矩阵

- 考虑突变的比对
- 蛋白的自然突变率矩阵PM1
- 矩阵自相乘得到外推矩阵 PM10 PM250 取对数为打分矩阵
- 取不同矩阵源于研究目的对多样性的判断

## E 值

- 表示序列的同源性 比对得分的稀有性
- 两个参数 数据库大小(N) 比对得分(S) E = N/S 
- 数据库越大越可能随机碰到相同序列 得分越高越可能同源
- E值很小说明同源性很高 E值很大什么说明不了
- 一般阈值1e-04

## PSI-BLAST

- 先用BLAST在一定E值上建库
- 计算新库的氨基酸概率 再与全库比对得分 得到统计显著性
- 可以发现BLAST未发现的序列 建立蛋白家族

## 蛋白

- Profiles 定量描述 
- Patterns 定性描述
- Signature 蛋白保守序列
- motif 少于20个氨基酸 指示二级结构
- Domains 超过40个氨基酸 蛋白的球状区
- 共同点 保守
- 正则表达式表示保守区
    - E-X(2,4)-[FHM]-X(4)-{P}-L
    - E后随意两个，三个，四个然后FHM其中一个，然后随意四个，然后一个不是P，最后为L
    - 可以精确可以模糊
    - 没有E值

## 蛋白结构预测

- 分子量 道尔顿（Da）描述质量
- 等电点 蛋白不带电的pH值
    - 小于7 酸性 中性带负电
    - 大于7 碱性 中性带正点
- 网站[计算](http://web.expasy.org/compute_pi/)
- 蛋白定位 分泌 胞内 核内
    - MITOPRED 预测线粒体蛋白

## 细菌基因组

- 细菌是环形DNA 真核是线性染色体
- 细菌不加工mRNA
- 细菌一段mRNA上有多个顺反子 也就是多个编码DNA序列
- 操纵子在mRNA编码的上游或下游调控转录
- [GLIMMER](http://www.ncbi.nlm.nih.gov/genomes/MICROBES/glimmer_3.cgi)与[FGENESB](http://www.softberry.com/berry.phtml?topic=fgenesb&group=programs&subgroup=gfindb)用来预测一段序列的转录情况

## 病毒

- 三种 RNA DNA 逆转录病毒 突变快
- RNA病毒三种 双链 正链 负链
- 逆转录基因组简单 Gag Pol Env
- 凝集素等决定病毒亚型

## 单核苷酸多态性（SNP）

- 至少1%种群中存在的DNA单核苷酸变化
- 后果
    - 编码区改变影响表型
    - 不改变蛋白序列的编码区可能影响mRNA加工
    - 启动子或调控区可能影响表达
    - 其他区没有影响 可作为染色体标记- 类型
    - 不改变氨基酸
    - 改变氨基酸
    - 非编码区
- 数据库
    - [dbSNP](http://www.ncbi.nlm.nih.gov/SNP/)
    - [SNPEffect](http://snpeffect.switchlab.org/) SNPs对蛋白的影响
    - [SNPedia](http://www.snpedia.com/index.php/SNPedia) SNPs的临床效应
    - [1000 基因组外显子计划](http://www.ncbi.nlm.nih.gov/pubmed/23128226) 第二代测序的发展

## 真核基因预测

- CDS是mRNA的子集
- CDS可能比mRNA外显子少
- 基因预测只能发现编码区外显子
- 有些转录变化不改变蛋白序列：UTR区与同义密码子

## DNA指纹

- 重复 突变会影响限制性片段长度
- VNTR 用来排除嫌犯
- PCR 用来扩增相关片段
- [CODIS](http://www.fbi.gov/about-us/lab/biometric-analysis/codis) 区域在美国用来鉴定身份

## Ensembl

- 外显子基因组学[数据库](ensembl.org)
- 可选择人类 鼠 斑马鱼等常见物种

## 基因组学数据分析

- [主页](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
```

## 链接 {#biolink}

- [Michael Love的教案](https://biodatascience.github.io/compbio/)
- [生信前沿信息集](http://www.gettinggeneticsdone.com/2017/02/staying-current-in-bioinformatics-genomics-2017.html)