chap二代测序数据分析实用实用教案

《chap二代测序数据分析实用实用教案》由会员分享，可在线阅读，更多相关《chap二代测序数据分析实用实用教案（51页珍藏版）》请在装配图网上搜索。

1、BLAST第1页/共50页第一页，共51页。BLAST第2页/共50页第二页，共51页。二代(r di)测序数据分析2012/9/25第3页/共50页第三页，共51页。测序数据分析流程(lichng) 短片(dun pin)段序列 短片(dun pin)段Mapping 基因型估计 关联分析短序列片段短序列片段(fastq)Mapping (BWA, Bowtie)基因型估计基因型估计 (Bayes)关联分析关联分析第4页/共50页第四页，共51页。短序列(xli)片段 Fastq格式 文本文件 包含(bohn)短片段序列（reads）和碱基测序质量 示例 HWI-EAS209_0006_FC

2、706VJ:5:58:5894:21141#ATCACG/1 TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT +HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1 efcfffffcfeefffcffffffddffeed_Ba_YBBBBBBBBBBRTTddddddddddaddBBBBBBBBBBBBBBBBBBBBBBBB第5页/共50页第五页，共51页。短序列(xli)片段 长度(chn

3、gd)：35 - 120bp 碱基测序质量：Phred Quality Score ASCII表示 P是base-calling错误率1010logQP 第6页/共50页第六页，共51页。ASCII码表第7页/共50页第七页，共51页。内容(nirng) 短片段序列 短片段Mapping 基因型估计 关联(gunlin)分析第8页/共50页第八页，共51页。短片(dun pin)段Mapping 输入 一个(y )参考基因组 大量(10-1000M)的25100bp的reads 输出 成功map到参考基因组上的每一个(y )位点信息 未成功map比例第9页/共50页第九页，共51页。短片(du

4、n pin)段Mapping 问题 不唯一位置 不确切位置 方法 哈希表：迅速，需要完美(wnmi)匹配 阵列扫描：无法处理空隙 动态规划：Indels，最优，但速度慢 Burrows-Wheeler Transform (BW Transform): 快速，但对于空隙和错配，缺乏敏感度第10页/共50页第十页，共51页。BW Transform X Bacaacg$acaacgaacg$acacaacg$acg$acacaacg$acg$acaag$acaacgc$aaacBurrows-Wheeler Matrix (BWM)BWT第11页/共50页第十一页，共51页。BW Transfo

5、rm 循环(xnhun)转换0 0a ac c a a a a c c g g $ $6 6 $ $ a a c c a a a a c c g g1 1 c c a a a a c c g g $ $ a a2 2 a a a a c c g g $ $ a a c c2 2 a a a a c c g g $ $ a a c c0 0 a a c c a a a a c c g g $ $3 3 a a c c g g $ $ a a c c a a3 3 a a c c g g $ $ a a c c a a4 4 c c g g $ $ a a c c a a a a1 1 c c a

6、 a a a c c g g $ $ a a5 5 g g $ $ a a c c a a a a c c4 4 c c g g $ $ a a c c a a a a6 6 $ $ a a c c a a a a c c g g5 5 g g $ $ a a c c a a a a c cn字典(zdin)方式排序第12页/共50页第十二页，共51页。LF mapping LF (last-first) mapping: 字符在最后(zuhu)一列与第一列的出现次数位置保持不变SASA seqseq0 06 6$ $a ac ca aa ac cg g0 01 12 2a aa ac cg

7、g$ $a ac c1 12 20 0a ac ca aa ac cg g$ $2 23 33 3a ac cg g$ $a ac ca a3 34 41 1c ca aa ac cg g$ $a a4 45 54 4c cg g$ $a ac ca aa a5 56 65 5g g$ $a ac ca aa ac c6 6$ $ 0 0a a 1 1c c 4 4g g 6 6第13页/共50页第十三页，共51页。LF mapping LFi = CLi + ri，CLi是Li在BWT中的起始(q sh)位置，ri是Li出现的次数 LF0 = CL0 + 1 = 6 + 1 -1 = 6 L

8、F1 = CL1 + 1 = 4 + 1 -1= 4 LF2 = CL2 + 1 = 0 + 1 -1= 0 LF3 = CL3 + 1 = 1 + 1 -1= 1 LF4 = CL4 + 2 = 1 + 2 -1= 2 LF5 = CL5 + 2 = 1 + 3 -1= 3 LF6 = CL6 + 2 = 4 + 2 -1= 5 LF = 6 4 0 1 2 3 5第14页/共50页第十四页，共51页。BW Transform 重构 恢复原序列(xli): X = L LF(i) + X; i = LF(i)g g$ $ a a c c a a a a c c g ga a a a c c

9、g g $ $ a a c ca a c c a a a a c c g g $ $a a c c g g $ $ a a c c a ac c a a a a c c g g $ $ a ac c g g $ $ a a c c a a a ag g $ $ a a c c a a a a c cc c g g$ $ a a c c a a a a c c g ga a a a c c g g $ $ a a c ca a c c a a a a c c g g $ $a a c c g g $ $ a a c c a ac c a a a a c c g g $ $ a ac c g g

10、$ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c ca a c c g g$ $ a a c c a a a a c c g ga a a a c c g g $ $ a a c ca a c c a a a a c c g g $ $a a c c g g $ $ a a c c a ac c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c cLF0=6, L6=ci=6LF6=5, L5=ai=5X=L0i=0第15页/共5

11、0页第十五页，共51页。BW Tranform 重构 恢复原序列(xli)a a a a c c g g$ $ a a c c a a a a c c g ga a a a c c g g $ $ a a c ca a c c a a a a c c g g $ $3 3 a a c c g g $ $ a a c c a ac c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c cc c a a a a c c g g$ $ a a c c a a a a c c g g1

12、 1 a a a a c c g g $ $ a a c ca a c c a a a a C C g g $ $3 3 a a c c g g $ $ a a c c a ac c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c ca a c c a a a a c c g g$ $ a a c c a a a a c c g g1 1 a a a a c c g g $ $ a a c ca a c c a a a a c c g g $ $3 3 a a c c g g

13、 $ $ a a c c a a4 4 c c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c cLF5=3, L3=ai=3LF3=1, L1=ci=1LF1=4, L4=ai=4第16页/共50页第十六页，共51页。BW Transform BW Transform是可逆的 时间(shjin)复杂度为O(n) 重构时间(shjin)复杂度为O(n)，空间复杂度为O(n)第17页/共50页第十七页，共51页。BWT mapping 查询(chxn)序列aacSASAa a

14、a a c c0 0 $ $ a a c c a a a a c c g g1 1 a a a a c c g g $ $ a a c c2 2 a a c c a a a a c c g g $ $3 3 a a c c g g $ $ a a c c a a4 4 c c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c ca a a a c c0 0 $ $ a a c c a a a a c c g g1 1 a a a a c c g g $ $ a a c c2 2

15、 a a c c a a a a c c g g $ $3 3 a a c c g g $ $ a a c c a a4 4 c c a a a a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c ca a a a c c0 0 $ $ a a c c a a a a c c g g1 1 a a a a c c g g $ $ a a c c2 2 a a c c a a a a c c g g $ $3 3 a a c c g g $ $ a a c c a a4 4 c c a a a

16、 a c c g g $ $ a a5 5 c c g g $ $ a a c c a a a a6 6 g g $ $ a a c c a a a a c c第18页/共50页第十八页，共51页。BWT mapping 对于参照序列X和查询序列W，计算( j sun)W的SA区间 SA区间可以通过从末端循环计算( j sun)得到 R(aW)为SA区间下限 R(aW)为SA区间上限 C(a)：X0,n-2排序后a前面的字符个数 O(a,i)：a在B0,i中出现的个数()( )( , () 1)()( )( , () 1R aWC aO a R WR aWC aO a R W第19页/共50页

17、第十九页，共51页。BWT mapping cSASAa a a a c c BWTBWT a a c c g g0 0 6 6 $ $ a a c c a a a a c c g gg g 0 0 0 0 1 11 1 2 2 a a a a c c g g $ $ a a c cc c 0 0 1 1 1 12 2 0 0 a a c c a a a a c c g g $ $ $ 0 0 1 1 1 13 3 3 3 a a c c g g $ $ a a c c a aa a 1 1 1 1 1 14 4 1 1 c c a a a a c c g g $ $ a aa a 2 2 1

18、 1 1 15 5 4 4 c c g g $ $ a a c c a a a aa a 3 3 1 1 1 16 6 5 5 g g $ $ a a c c a a a a c cc c 3 3 2 2 1 1$ $ a ac c g grankrank 0 0 1 14 4 6 6cfrom 4+0=4to 4+2-1=5第20页/共50页第二十页，共51页。BWT mapping acSASAa a a a c c BWTBWT a a c c g g0 0 6 6 $ $ a a c c a a a a c c g gg g 0 0 0 0 1 11 1 2 2 a a a a c c

19、 g g $ $ a a c cc c 0 0 1 1 1 12 2 0 0 a a c c a a a a c c g g $ $ $ 0 0 1 1 1 13 3 3 3 a a c c g g $ $ a a c c a aa a 1 1 1 1 1 14 4 1 1 c c a a a a c c g g $ $ a aa a 2 2 1 1 1 15 5 4 4 c c g g $ $ a a c c a a a aa a 3 3 1 1 1 16 6 5 5 g g $ $ a a c c a a a a c cc c 3 3 2 2 1 1a ac cg grankrank1 1

20、4 46 6afrom 1+1=2to1+3-1=3第21页/共50页第二十一页，共51页。BWT mapping aac所在参照序列(xli)中起始位置：2 （0起始）SASAa a a a c c BWTBWT a a c c g g0 0 6 6 $ $ a a c c a a a a c c g gg g 0 0 0 0 1 11 1 2 2 a a a a c c g g $ $ a a c cc c 0 0 1 1 1 12 2 0 0 a a c c a a a a c c g g $ $ $ 0 0 1 1 1 13 3 3 3 a a c c g g $ $ a a c c

21、a aa a 1 1 1 1 1 14 4 1 1 c c a a a a c c g g $ $ a aa a 2 2 1 1 1 15 5 4 4 c c g g $ $ a a c c a a a aa a 3 3 1 1 1 16 6 5 5 g g $ $ a a c c a a a a c cc c 3 3 2 2 1 1a ac cg grankrank14 46 6afrom 1+0=1to1+1-1=1第22页/共50页第二十二页，共51页。BWT mappingSASABWTBWTa ag gc ct t0 01111$ $c ct tg ga aa aa ac ct tg

22、 gg gt t0 00 00 01 11 13 3a aa aa ac ct tg gg gt t$ $c ct tg g0 01 10 01 12 24 4a aa ac ct tg gg gt t$ $c ct tg ga a1 11 10 01 13 35 5a ac ct tg gg gt t$ $c ct tg ga aa a2 21 10 01 14 40 0c ct tg ga aa aa ac ct tg gg gt t$ $2 21 10 01 15 56 6c ct tg gg gt t$ $c ct tg ga aa aa a3 31 10 01 16 62 2g ga

23、 aa aa ac ct tg gg gt t$ $c ct t3 31 10 02 27 78 8g gg gt t$ $c ct tg ga aa aa ac ct t3 31 10 03 38 89 9g gt t$ $c ct tg ga aa aa ac ct tg g3 32 20 03 39 91010t t$ $c ct tg ga aa aa ac ct tg gg g3 33 30 03 310101 1t tg ga aa aa ac ct tg gg gt t$ $c c3 33 31 13 311117 7t tg gg gt t$ $c ct tg ga aa aa

24、 ac c3 33 32 23 3n示例(shl)2第23页/共50页第二十三页，共51页。BWT mapping ctgg所在起始(q sh)位置：6 （0起始(q sh)）g gfrom from 6+0=66+0=6toto6+3-1=86+3-1=8g gfrom from 6+1=76+1=7toto6+2-1=76+2-1=7t tfrom from 9+2=119+2=11toto9+3-9+3-1=111=11c cfrom from 4+1=54+1=5toto4+2-1=54+2-1=5a ag gc ct trankrank1 14 46 69 9BWTBWTa ag g

25、c ct tSASA1111$ $c ct tg gt t0 00 00 01 11 13 3a aa aa ac cg g0 01 10 01 12 24 4a aa ac ct ta a1 11 10 01 13 35 5a ac ct tg ga a2 21 10 01 14 40 0c ct tg ga a$ $2 21 10 01 15 56 6c ct tg gg ga a3 31 10 01 16 62 2g ga aa aa at t3 31 10 02 27 78 8g gg gt t$ $t t3 31 10 03 38 89 9g gt t$ $c cg g3 32 20

26、 03 39 91010t t$ $c ct tg g3 33 30 03 310101 1t tg ga aa ac c3 33 31 13 311117 7t tg gg gt tc c3 33 32 23 3第24页/共50页第二十四页，共51页。BWT mapping 优势 需要较少的存储空间，3GB的参照序列需要46空间 快速(kui s)，线性时间复杂度 相关软件 BWA Bowtie第25页/共50页第二十五页，共51页。BWT based mapping 保存整个参照序列 一个(y )碱基一个(y )碱基的进行联配 如果没有发现匹配，把碱基替换成其他碱基重新匹配第26页/共50

27、页第二十六页，共51页。Mapping结果(ji gu)第27页/共50页第二十七页，共51页。内容(nirng) 短片段序列(xli) 短片段Mapping 基因型估计 关联分析第28页/共50页第二十八页，共51页。基因组多态 寡核苷酸多态：SNP (single nucleotide polymorphism) 插入(ch r)删除：Indel (Insertion/Deletion) 短片段重复：STR (short tandem repeat) 拷贝数目多态：CNV (copy number variation)第29页/共50页第二十九页，共51页。多态位点信息(xnx) 测序深度

28、(shnd)：Reads个数 （20X） 碱基：A、T、C、G、- 碱基reads：10、10、0、0、0 基因型：纯合子，杂合子第30页/共50页第三十页，共51页。基因型估计(gj) 简单统计： A(10 reads), T(10 reads) AT A(10 reads), T(2 reads) AA A(2 reads), T(10 reads) TT 贝叶斯估计(gj) 根据先验概率来估计(gj)三种基因型的概率，选取概率最大的作为位点基因型第31页/共50页第三十一页，共51页。基因型估计(gj) 贝叶斯公式是由英国学者贝叶斯发表的一篇论文论归纳推理的一种方法提出(t ch) 先验

29、概率：所有先验假设的信息 观察值：观察到的信息 后验概率：根据先验概率和观察值推断的结果(| )( )* (|)/ ( )P G RP GP R GP R第32页/共50页第三十二页，共51页。基因型估计(gj) SNP 1 : 6 reads (4 T, 2 A) 假设(jish)测序错误率为1%，人群中A出现的概率为75%，T出现的概率为25% |R|=6, |R(T)|=2, |R(A)|=4 Priorp(G=TA) = 2*0.75*0.25=0.375 Priorp(G=AA) = 0.75*0.75=0.5625 Priorp(G=TT) = 0.25*0.25=0.0625第3

30、3页/共50页第三十三页，共51页。基因型估计(gj) P(R|G=TA) = 0.56=0.015625 P(R|G=AA) = 0.012*0.994=0.00009605 P(R|G=TT) = 0.014*0.992=0.000000009801 P(R,G=TA) = P(R,G=TA) = 0.015625*0.375 P(R,G=AA) = P(R,G=AA) = 0.000096059601*0.5625 P(R,G=TT) = P(R,G=TT) = 0.000000009801*0.0625 P(R) = P(R,G=TA)+P(R,G=AA)+P(R,G=TT) = 0.

31、005913第34页/共50页第三十四页，共51页。基因型估计(gj) Postp(G=TA|R) = P(G=TA,R)/P(R) = 0.375*0.015625/0.005913 = 0.9909 Postp(G=AA|R) = P(G=AA,R)/P(R) = 0.5625*0.00009605/0.005913 = 0.009137 Postp(G=TT|R) = P(G=TT,R)/P(R) = 0.0625*0.000000009801/0.005913 =.0000001036 基因型：AT第35页/共50页第三十五页，共51页。内容(nirng) 短片段序列(xli) 短片段

32、Mapping 基因型估计 关联分析第36页/共50页第三十六页，共51页。关联(gunlin)分析 假设检验( jinyn) 零假设H0：SNP与疾病不相关 H1：SNP与疾病相关 检验( jinyn)统计量 质量性状：卡方检验( jinyn) (Pearson chi-square test)、Fisher确切检验( jinyn) (Fisher exact test) 数量性状：回归分析、方差分析第37页/共50页第三十七页，共51页。关联(gunlin)分析 根据统计量大小及相应自由度，求得p值 P值反映了零假设(jish)成立概率，p1 表明是危险因素 OR (odds ratio)

33、 1 表明是保护因素A AT T总计比例casecase400400800800120012000.460.46controlcontrol10001000400400140014000.540.54总计140014001200120026002600比例0.540.540.460.46第39页/共50页第三十九页，共51页。关联(gunlin)分析 卡方检验(jinyn)： X-squared = 377.3243, df = 1, p-value 1，危险因素；置信区间1，保护因素195%ln()ln() 1.96*0.5CI ofORORABCD第40页/共50页第四十页，共51页。基因

34、型关联(gunlin)分析 共显性模型 AA vs AT vs TT 考察不同(b tn)基因型对疾病的影响 分组：AA=3, AT=2, TT=1AAAAATATTTTTcasecase100100200200300300controlcontrol400400200200100100X-squared = 273.9286, df = 2, p-value 2.2e-16X-squared = 273.9286, df = 2, p-value 2.2e-16第41页/共50页第四十一页，共51页。基因型关联(gunlin)分析 显性模型(mxng) (AA+AT) vs TT 考察显性对

35、疾病的影响 分组：AA=2, AT=2, TT=1AA+ATAA+ATTTTTcasecase300300300300controlcontrol600600100100X-squared = 193.4524, df = 1, p-value 2.2e-X-squared = 193.4524, df = 1, p-value 2.2e-1616第42页/共50页第四十二页，共51页。基因型关联(gunlin)分析 隐性模型 AA vs (AT+TT) 考察隐形对疾病( jbng)的影响 分组：AA=2, AT=1, TT=1AAAAAT+TTAT+TTcasecase10010050050

36、0controlcontrol400400300300X-squared = 223.631, df = 1, p-value 2.2e-X-squared = 223.631, df = 1, p-value 2.2e-1616第43页/共50页第四十三页，共51页。多重检验(jinyn)校正 当多次进行假设检验( jinyn)时，通过调整p值来校正假阳性 多重检验( jinyn)校正方法第44页/共50页第四十四页，共51页。多重检验(jinyn)校正 Bonferroni correction 计算原有p值 Corrected P-value= p-value * n (number o

37、f tests)，小于0.05为显著 Bonferroni Step-down (Holm) correction 计算原始p值 对所有(suyu)p值从小到大排序 对于在i位置的test，Corrected P-value= p-value * (n-i+1) 直到位置k，k之后的corrected p-value大于0.05第45页/共50页第四十五页，共51页。多重比较(bjio)校正 Westfall and Young 置换 计算原有p值并排序(pi x) 把原始数据随机分成case组和control组，构成伪数据 在伪数据的基础上计算p值 最小的伪p值被保留下来，并与原始p值做比较

38、 重复这个过程，如果伪p值小于真是p值，则作为校正后的p值第46页/共50页第四十六页，共51页。多重检验(jinyn)校正 Benjamini and Hochberge FDR 计算原始p值 从小到大排序p值 对于排序i的检验( jinyn)，Corrected p-value = p-value*(n/n-i+1) ，corrected p-value0.05为显著第47页/共50页第四十七页，共51页。多重检验(jinyn)校正第48页/共50页第四十八页，共51页。 Thanks for your attention! Questions? Happy Holiday!第49页/共50页第四十九页，共51页。感谢您的欣赏(xnshng)！第50页/共50页第五十页，共51页。NoImage内容(nirng)总结BLAST。第1页/共50页。第2页/共50页。2012/9/25。第3页/共50页。第4页/共50页。HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1 TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT第五十一页，共51页。