第14章 生物信息
14.2 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 3D结构数据库
- Genome 基因组数据库
- OMIM 人类孟德尔遗传数据库 用来探索等位基因问题
- 分类数据库 用来界定分类
- GEO 基因芯片的实验数据
- SNP 基因指纹数据库
14.3 动态规划
- 用于序列比对
- 对角线得分 按总分评价比对结果
- 可全局 可局部
- 序列比对指标是特异性与相似性
- 特异性指精确匹配比率
- 相似性指精确匹配加化学相似性比率 结构相近则相似
- FASTA 慢准 BLAST 快
- 三种情况 匹配 不匹配 间隔
- 间隔罚分
14.5 E 值
- 表示序列的同源性 比对得分的稀有性
- 两个参数 数据库大小(N) 比对得分(S) E = N/S
- 数据库越大越可能随机碰到相同序列 得分越高越可能同源
- E值很小说明同源性很高 E值很大什么说明不了
- 一般阈值1e-04
14.7 蛋白
- Profiles 定量描述
- Patterns 定性描述
- Signature 蛋白保守序列
- motif 少于20个氨基酸 指示二级结构
- Domains 超过40个氨基酸 蛋白的球状区
- 共同点 保守
- 正则表达式表示保守区
- E-X(2,4)-[FHM]-X(4)-{P}-L
- E后随意两个，三个，四个然后FHM其中一个，然后随意四个，然后一个不是P，最后为L
- 可以精确可以模糊
- 没有E值
14.8 蛋白结构预测
- 分子量 道尔顿（Da）描述质量
- 等电点 蛋白不带电的pH值
- 小于7 酸性 中性带负电
- 大于7 碱性 中性带正点
- 网站计算
- 蛋白定位 分泌 胞内 核内
- MITOPRED 预测线粒体蛋白
14.11 单核苷酸多态性（SNP）
- 至少1%种群中存在的DNA单核苷酸变化
- 后果
- 编码区改变影响表型
- 不改变蛋白序列的编码区可能影响mRNA加工
- 启动子或调控区可能影响表达
- 其他区没有影响 可作为染色体标记- 类型
- 不改变氨基酸
- 改变氨基酸
- 非编码区
- 数据库
- dbSNP
- SNPEffect SNPs对蛋白的影响
- SNPedia SNPs的临床效应
- 1000 基因组外显子计划 第二代测序的发展
14.13 DNA指纹
- 重复 突变会影响限制性片段长度
- VNTR 用来排除嫌犯
- PCR 用来扩增相关片段
- CODIS 区域在美国用来鉴定身份
14.14 Ensembl
- 外显子基因组学数据库
- 可选择人类 鼠 斑马鱼等常见物种
14.15 基因组学数据分析
14.15.3 数据分析应用背景
14.15.4 Bioconductor
- 官方说明
- 使用
biocLite()安装，安装后仍需要library()才能使用
source("http://bioconductor.org/biocLite.R")
biocLite()14.15.4.1 数据结构
14.15.4.1.2 表型数据
- S 行 V 列 V代表样本特征，为分类或连续变量
- 如果表型数据解释不清，可以建立一个解释样本特征的labelDescription数据框，通过
phenoData <- new("AnnotatedDataFrame",data=pData, varMetadata=metadata)建立AnnotatedDataFrame类型数据
14.15.4.1.4 组装数据
- 将分析数据，表型数据，实验描述组装为一个ExpressionSet类型的对象
exampleSet <- ExpressionSet(assayData=exprs,phenoData=phenoData,experimentData=experimentData,annotation="hgu95av2")- annotation代表了一组相似实验设计的芯片数据的代号，通过相关代号可以索引到芯片特征信息并将其与其他数据如基因型，染色体位置等连接以便于分析
- 从ExpressionSet里可以按照表型数据提取子集，也就是对 S 截取 V 中特定子集
exampleSet[ , exampleSet$gender == "Male"] esApply用来针对ExpressionSet应用函数
14.15.4.1.5 数据集应用
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，用火山图来表示
14.15.4.1.5.1 Iranges
library(IRanges)序列范围ir <- IRanges(start = c(3, 5, 17), end = c(10, 8, 20))定义序列IRanges(5, 10)表示5到10这6个碱基对，可以shiftrange(ir)表示存在ir中序列的起止范围gaps(ir)表示寻找ir中间隔片段disjoin(ir)表示将ir中序列碎片化后互不重叠的片段
14.15.4.1.5.2 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)定义该基因库类型中的列并赋值
14.15.4.2 数据读取
- microarray 或 NGS 数据由芯片厂商提供，常见读取原始信息的包有affyPLM、affy、 oligo、limma
- 在Bioconductor里，这些原始数据要转为ExpressionSet格式
14.15.4.2.1 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)对样本进行背景校正与正则化，从探针层转化为基因层数据
14.15.4.2.2 背景干扰
- spikein方法 梯度加入已知浓度的基因片段 阵列上进行shift 类似拉丁方设计
- 可以看到同一基因不同片段大致符合先平后增模式 开始阶段是噪声主导 后面是浓度主导
- 使用类似基因模拟噪声主导 相减后得到去干扰浓度效应 但低值部分会导致方差过大
- 也可以使用统计建模方法模拟背景值与响应 得到还原度更高的信号
14.15.5 示例：甲基化数据分析
14.15.5.1 读取数据
devtools::install_github("coloncancermeth","genomicsclass")
library(coloncancermeth)
data(coloncancermeth)该数据集为结肠癌病人与对照的DNA甲基化数据集。
14.15.5.2 数据说明
dim(meth)
dim(pd)
length(gr)meth为测序数据，pd为样本信息，gr测序片段信息。
colnames(pd)
table(pd$Status)
X = model.matrix(~pd$Status)查看病患与正常人的分组并构建模型。
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）下甲基化基因聚类的后出现的位置，聚类号，聚类相关性等信息寻找问题基因，可从中提取相关信息
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岛。
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