生物信息学 第十一章转录组与转录调控分析 (2)
基因集富集分析 Gene Set Enrichment Analysis (GSEA) 考虑所有的基因而不仅是差异表达基因 根据差异表达水平将所有基因排序 发现基因集里的基因是否倾向于分布在所有基因列表的前部 ( 上调 ) 或后部 ( 下调 ) 计算流程 : 计算一个富集分值 (Enrichment score, ES) 根据 ES 分值估计统计显著性 针对多重检验进行校正
GSEA 的输入 基因表达数据 D, 包括 N 个基因和 k 个样本 排序产生基因列表 L 上调 - 不变 - 下调 差异表达基因的 p-value 用指数 p 来控制每一步的权重 p=0, 标准 Kolmogorov Smirnov 统计 p=1,gsea 软件 独立的基因集 S, 如 GO, KEGG 包括 N H 个基因
Enrichment Score ES(S) 基因集 D 排序为 L= {g 1,..., g N },r(g j ) = r j 对基因列表 L, 计算每一个位置 i 以上, 包括在基因集 S 中的 ( hits ) 部分, 以及不在 S 中的 ( misses ) 部分 ES = P hit -P miss 的最大值
ES 计算 : 实例 基因列表 Gene E-ratio 1 5 2 4 3 3 4 2 5 1.5 6 0.04 7 0.007 基因集 S= 1, 3, 4, 5 P hit (S, 1) = 5/5=1 P miss (S, 1) = 0 P hit (S, 2) = 5/5=1 P miss (S, 2) =0.33 P hit (S, 3) = 1 + 4/(5+4) =1.44 P miss (S, 3) =0.33 P hit (S, 4) = 1.44 + 3/(5+4+3)=1.69 P miss (S, 4) = 0.33 ES(1) = 1 ES(2) = 0.67 ES(3) = 1.11 ES(4) = 1.36
ES 计算 : 实例 基因列表 Gene E-ratio 1 5 2 4 3 3 4 2 5 1.5 6 0.04 7 0.007 基因集 S= 1, 3, 4, 5 P hit (S, 5) = 1.69+ 1.5/(5+4+3+2+1.5)=1.8 P miss (S, 5) = 0.33 ES(5) = 1.47 P hit (S, 6) = 1.8 ES(6) = 1.14 P miss (S, 6) = 0.66 P hit (S, 7) = 1.8 ES(7) = 0.8 P miss (S, 7) = 1
显著性估算 计算 ES 的显著性 : 与 ES NULL 比较 将原数据的 p-value 取出, 随机分给每个基因, 重新排序 重新计算 ES 值, 重复 1,000 次置换 (Permutation), 建立 ES NULL 分布的直方图 估算 p-value = 大于或等于 ES 的个数 / 总数 Null distribution of enrichment scores Actual ES
GSEA 算法总结
Kolmogorov Smirnov 检验 GSEA, p = 0, 不考虑 p-value 的权重 非参检验方法 (Nonparametric test) 检验分布是否符合正态分布 两样本的 Kolmogorov Smirnov 检验 差异最大的值算显著性
基因调控网络 早期观点 : 表达谱相似的基因可能存在功能上的关联, 可能相互作用 ( 直接作用 ) 当前的观点 : 表达谱相似的基因可能具有共同的调控元件 ( 基因 UTR 区域存在共同的启动子 ), 能够被同一个上游因子所调控
相关系数 : 基因共表达网络 与光合效率和气孔发育相关的基因 :ERL2 在 Wild-type 中与之显著相关, 但在 Mutant 中显著不相关的基因 Wild-type Mutant ERL2 SKP1 Unknown ChS1
相关系数 : 基因共表达网络
Microarray: 工具 & 数据库
GEO - NCBI
Array Express - EMBL
SMD - Stanford
RNA-Seq 分析流程 分离纯化所有的 mrna 利用反转录酶将 mrna 转变成 cdna cdna 测序 将序列映射到基因组上 序列检测到越多, 则表达越高 RNA-Seq 数据分析 : 挑战! 将读段映射到参考基因组上 定量已知的基因 发现新的转录本 定量可变剪切异构体
RNA-Seq 数据 : 序列读段 FASTQ: Illumina 测序读段文件 Line 1: Line 2: Line 3: Line 4: @EAS042_0001:1:1:1061:20798#0/1 TNTCTGTGTCCTGGGGCATCAATGATAGTCACATAGTACTTGCTGGTCTCAAATTTCCACAAGGAGATATCAATGG +EAS042_0001:1:1:1061:20798#0/1 ab\^^y]a^]cde`daayaaa_bc\\`b^y\a\aauqy\]a\`aa\w ]HVZ]VQF^[`UH]\J^F^T^\\I] Line 1 EAS042_0001 the unique instrument name 1 flowcell lane 1 tile number within the flowcell lane 1061 20798 #0 'x'-coordinate of the cluster within the tile 'y'-coordinate of the cluster within the tile index number for a multiplexed sample (0 for no indexing) Line 2: 原始序列 Line 3: +? Line 4: 序列的质量分值 -5 ~ 62 使用 ASCII 59 ~ 126 /1 the member of a pair, /1 or /2 (paired-end or mate-pair reads only)
RNA-Seq 分析流程 Poly(A) end-reads 样品准备 高通量测序 数据分析 : 读段映射 可视化 重头组装 定量
RNA-Seq vs. DNA 芯片 RNA-seq 可以发现新的转录本和可变剪切异构体,DNA 芯片只能检测已知转录本的表达 RNA-Seq 比全基因组铺瓦芯片的解析度更高 : 单个碱基 RNA-Seq 相同的实验流程可以做不同的分析 单核苷酸多态性 (SNP 芯片 ) 外显子连接点 ( 连接点芯片 ) 基因融合 ( 基因融合芯片 )
RNA-Seq 的优势 转录组分析方法 :RNA-Seq 比其他方法更有优势
RNA-Seq vs. DNA 芯片 营养丰富的培养基里的芽殖酵母 仅在中等表达时, 两者有较高一致性 基因表达低或者高的时候, 关联性差 DNA 芯片的对基因表达低或高的检测不够敏感
RNA-Seq 的问题 : 建库 mrna 与 cdna 碎裂的不同 利用 oligo-dt 引物得到的 cdna 更倾向于获得转录本的 3 端 RNA 碎裂在 5 和 3 端都较少
RNA-Seq 的准确性 HiSeq 2000 和 SOLiD 测序分析结果 HuGene2:Affymetrix 芯片数据 利用 qpcr (Taqman 或 PrimePCR) 验证表达水平 不同方法一致性较高
RNA-Seq 数据 : 基因模型 基因 C17orf45 Bioinformatics, (ENSG00000175061) 2018, HUST 的结构
RNA-Seq 数据的归一化 不同样本之间比较 : 测序深度不同 不同基因之间 : 表达相同时序列越长读段越多 RPKM (The number of reads per kilobase of transcript per million mapped reads) 每百万可映射读段每千 kb 的读段数 例如, 测序 10 7 读段,8*10 6 可映射 (N) 1kb 转录本 (L), 总共有 1000 个读段 (C) RPKM = 1000/(1 * 8) = 125 Multiread ( 多重定位读段 ) 计数时需要除以定位个数
归一化的 RNA-Seq 数据 技术重复的结果一致性较高 体外合成的转录本加入样本中检测 RPKM 与加入量线性相关
多重定位读段 RNA-Seq vs. DNA 芯片 不考虑多重定位读段造成部分基因 RPKM 值偏低 考虑多重定位读段提高与芯片结果的线性关系
归一化 : 不同样本之间的比较 不同样本相同基因的比较不需要考虑长度因素 Y gk : 基因 g 的读段数目,N k : 读段总数
差异表达基因分析
λ < 10: 有偏分布 λ < 10: 近似为正态分布 差异表达基因分析 利用 Fisher s exact test qpcr 与芯片结果拟合更好!? R=0.86 R=0.82
二次泊松分布? TSPM: two-stage Poisson model 测序可以看成是泊松过程 建库也可以近似看成是泊松过程 两者的 λ 不同, 需要分别估计 mrna 分子 读段 接头 连接扩增测序
负二项分布 :DESeq NB: Negative Binomial (Pascal) Distribution 每次试验的成功率 p, 失败率 1-p 经历 r 次失败后观察到 k 次成功事件, 则 RNA-Seq 数据的归一化 总共 m 个样本 基因 i 在样本 j 中检测到 K ij 个读段 样本 j 的采样深度 s j
负二项分布 :DESeq RNA-Seq 数据的归一化 总共 m 个样本 基因 i 在样本 j 中检测到 K ij 个读段 样本 j 的采样深度 s j 样本 j 的条件 p(j) q ip 基因 i 在条件 p 下的平均归一化表达量 用最大似然性估算 (MLE) 获得
DESeq: 差异表达基因分析 两样本的 p-value 计算 p(a, b) = Pr(K ia = a) Pr(K ib = b) K is = K ia + K ib
负二项分布 vs. 泊松分布 DESeq ( 橙 ),EdgeR ( 橙虚线 ), 泊松 ( 紫 ) NB 估算的方差更大
负二项分布 vs. 二次泊松分布 DESeq 和 EdgeR 的重合较高 :NB 分布 TSPM ( 二次泊松 ) 与其他工具重合较低
RNA-Seq 差异表达基因分析工具 目前最好 :DESeq 和 EdgeR DESeq: http://bioconductor.org/packages/release/bioc/html/deseq.html EdgeR: http://bioconductor.org/packages/release/bioc/html/edger.html R 语言 + Bioconductor 实现
可变剪切 AS, Alternatively Splicing 92 94% 的人类多外显子基因存在可变剪切现象 Major isoform vs. minor isoform 86% 包括至少有一个次要异构体 ( 比例 15%) SLC25A3 基因
果蝇 Dscam 基因的可变剪切 12*48*33*2=38,016 个不同的异构体 果蝇 : ~250,000 个神经元 轴突导向 (Axon guidance) 调控神经发育 神经元之间的连接
可变剪切的类型 8 种类型 仅统计外显子上的读段数 不全面
可变剪切的组织特异性 脑的可变剪切特异性最高 不同个体之间存在差异 细胞系有特别的剪切模式
可变剪切鉴定 :Cufflinks 用 Bowtie ( 领带 ) 将所有读段映射到基因组上 不能完美匹配的, 用 TopHat ( 礼帽 ) 发现剪切位点 http://cole-trapnell-lab.github.io/cufflinks/
可变剪切鉴定 :Cufflinks OLC 算法 (Overlap graph) 读段是图中的结点, 若可以匹配则置边 用 OLC 算法计算多条路径 用最大似然性方法估算每个异构体的表达水平
可变剪切鉴定 :Scripture DBG 算法 (Connectivity graph) 单个碱基是结点, 边是两个读段共有的序列 用 DBG 算法算出多个路径 Cuffdiff: 利用双端测序数据校正不同异构体的比例 RPKM 归一化表达 http://www.broadinstitute.org/software/scripture/
常用的 RNA-Seq 分析软件
表观遗传学 (Epigenetics) Dolly Ian Wilmut 多莉 : 生于 1996 年 7 月 5 日, 死于 2003 年 2 月 14 日
Dr. Evil Dr. Evil s Clone Mini Me Dr. Evil
The calico cat? Rainbow Allie Copy Cat
PEV: 位置效应斑 Position Effect Variegation white 基因 vs. 转座子
表观遗传学 基因的 DNA 序列不发生改变的情况下, 基因的表达水平与功能发生改变, 并产生可遗传的表型 特征 : (1) 可遗传 ;(2) 可逆性 ;(3) DNA 不变 表观遗传学机制 : DNA 甲基化 组蛋白修饰 MicroRNA 染色质重塑
Waddington's epigenetics 基因型 表型
DNA 甲基化 SAH, S-adenosylhomocysteine; SAM, S-adenosylmethionine
DNA 甲基化
叶酸 : 很重要?
DNA 甲基化 vs. 转座子 叶酸 : 很重要 IAP: intracisternal A particle 非甲基化 : 异常表达 Agouti 甲基化 : 正常表达
组蛋白的共价修饰
Paramutation 一个沉默的基因常常导致等位基因也发生沉默 maize P-rr gene:myb 转录因子
Paramutation 分子机理 (RNA-mediated trans-induction of chromatin):rna 参与调控
染色质重塑 染色质 (Chromatin): 染色体上高度致密的部分, 通常不表达基因
ChIP-Seq 转录因子 组蛋白等与 DNA 结合 Chromatin immunoprecipitation (ChIP) 结合高通量测序
ChIP-Seq 的读段分布 不同蛋白质与 DNA 结合的模式不同 A. 转录因子 ;B. RNA 聚合酶和组蛋白质等
ChIP-Seq 校正 测序读段的双峰校正为单峰 例如,CTCF 结合双峰的校正
常见的 ChIP-Seq 数据分析工具
不同软件结果获得的峰不同 人类三个转录因子 NRSF (neuron-restrictive silencer factor) GABP (growthassociated binding protein) FoxA1 (hepatocyte nuclear factor 3α)
qpcr 验证 NRSF: 大致相当 GABP:SISSRS 和 CisGenome 较差
DNA 甲基化的检测 : 芯片
MethylC-seq: Bisulfite Modification 重亚硫酸盐测序 : 非甲基化的 C 将被测序成 T, 反链为 A
Methyl-Seq 技术 Bisulfite-based MethylC-seq Reduced representation bisulfite sequencing (RRBS) Enrichment-based Methylated DNA immunoprecipitation sequencing (MeDIP-seq) Methylated DNA binding domain sequencing (MBD-seq) MethylC-seq 全基因组覆盖度高, 测 CpG 岛其他方法有优势