多序列比对

多序列比对 孟雪红 mengxuehong@genomics.cn Tel: +8600000000 January 2011

序列比对的意义 • 不同物种基因组共线性分析可以知道物种间亲缘关系，利于基因预测和功能注释(熊猫文章)

同一物种SD(片段复制)分析(蚂蚁文章)

主要内容 • 两物种基因组比对（lastz/chainnet） • 多物种基因组比对（multiz）

Target file Target sequence in put Repeat with reverse complement Scoring parameters Alignment output Lastz workflow Lastz/chainnet Indexing target seed words interpolation Back-end filtering Gapped extension HSP chaining Gap-free extension seeding Query file

Scoring Inference • HOXD70 http://genomewiki.ucsc.edu/index.php/Hg19_conservation_lastz_parameters

Indexing Target Seed Words ACGTGACATCACACATGGCGACGTCGCTTCAC target seed word position table |... | |10325-->12, 255, 28451, 36512 | |10326-->365, 5475, 47154, 225641 | |...

repeat • 如果知道repeat序列，将target和query序列在比对之前将repeats mark成小写字母。不参Indexing Target Seed Words步骤和seeding步骤。 • 如果repeat位点不知道，设置参数‑‑maxwordcount，在Indexing Target Seed Words步骤中将出现次数过多的seeds去掉。 • ‑‑masking比对过程中动态的mark掉比对多次的位点，只影响后续的query序列。

Seeding spaced seeds • seed=12of19（1110100110010101111） • seed=14of22（1110101100110010101111） target:ACGTGACATCACACATGGCGACGTCGCTTCACTGG query: GTAGCTTCAC GTAGCTTCAC pattern: 110 0 10111 1 110 010 1111

x-drop=910 X-drop extension | --> HSP ？ <-- | |-->seed<--| CACGAAACCAGCACGTATCCAAGGGACTATCCCC CATGGAACCAGCACGTATCCAAGGGCCTTTCCCC HSPs : high-scoring segment pairs Hsp-threshold=3000

Gapped Extension 1、仿射空位罚分 gap_open_penalty=400 gap_extend_penalty=30 公式：Wk=400+30k（k为gap长度） 2、y_drop=9400 3、gapped_threshold=3000 1 3 2

Needleman-Wunsch算法 1、DNA序列： S1 = GCCCTAGCG S2 = GCGCAATG 核苷酸替换打分矩阵S gap扣分d=-2 A T C G A 1 -1 -1 -1 T -1 1 -1 -1 C -1 -1 1 -1 G -1 -1 -1 1

2、算法规则 这个算法使用二维表格，一个序列沿顶部展开，一个序列沿左侧展开。通过以下三个途径到达每个单元格： 1）来自上面的单元格，代表将左侧的字符与空格比对。 2）来自左侧的单元格，代表将上面的字符与空格比对。 3）来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）

首先对第一行、第一列初始化

F(i-1,j-1)+s(xi,yj) F(i,j)=max F(i-1,j)-d F(i,j-1)-d

从右下角的单元格开始反向回溯，即可得到比对结果从右下角的单元格开始反向回溯，即可得到比对结果 S1' = GCCCTAGCG S2' = GCGC- AATG

Back-end Filtering • Identity • Continuity • Coverage • Match count

Interpolation

ChainNet • axtChain:将相邻的block连接起来，打分矩阵和blastz相同，gap打分改变。 • chainNet：对target序列，确定最优比对区域。 1、首先将所有的染色体或scaffold的碱基标记未用的。 2、按打分由高到低排序，形成列表。 3、迭代：每次从列表中取出一个chain，扔掉与已经存在的chain有overlap的区域，余下的部分添加上去，如果与之前的chain有gap，标记成子集，通过这种方式形成的层级称为net。记录overlap的区域，用于下一步识别重复 • netSyntenic：处理inersion、duplication。

运行 lastz_chainnet.py • step1_lastz_target_query.sh • step2_chain_target_query.sh • step3_net_target_query.sh 去除repeat，-M参数；切割文件（方式、数量）；切割脚本输出：maf格式

Multiz 提供信息： 1、物种的拓扑结构: ((t1 q1) q2) 2、两两物种lastz比对maf文件(以同一个物种为参考序列) 3、储存物种信息的list文件：

TBA(Threaded-Blockset Aligner) 将reference至于顶行，按照reference坐标对排列其余物种，按照系统发育树重新对行排列。 reference：h reference：m

打分：使用与lastz相同的核苷酸替换打分矩阵，每一列的打分为两两物种打分之和。打分：使用与lastz相同的核苷酸替换打分矩阵，每一列的打分为两两物种打分之和。 • Gap惩罚(quasi-natural)：400+30(L-1)

运行 python ../bin/run_multiz.py --pair_align pairwise_alignment1.list --tree "((t1 q1) q2)" --out `pwd`/output list文件： t1 q1 input/t1_q1.axt.chain.prenet.net.axt.maf t1 q2 input/t1_q2.axt.chain.prenet.net.axt.maf run_multiz.sh 输出maf文件

Thanks!

多序列比对

多序列比对

Presentation Transcript