Trinity 软件运行原理

2018-06-04

转录组 Trinity

通过阅读Trinity的源码，学习Trinity的运行过程以及使用的算法，写这篇博客的时候已经阅读Trinity的源码有一段时间了，但是觉得很多细节的地方还是理解不透彻，所以把能理解的部分先写出来，方便感兴趣的人一起学习

Trinity作为现阶段主流的转录组组装软件，主要有三部分组成Inchworm、Chrysalis、Butterfly。可以分为三个部分进行运行，而每一个部分又可以根据算法的实现拆分成不同的小块进行运行。主要的流程图如下:

Inchworm过程：这一步会读入所有的reads的fq文件，然后将fq的reads文件转换成fa格式的reads文件，并将3端的reads和5端的reads进行合并得到both.fa。然后利用jellyfish将reads打断成kmer，kmer默认是25bp的长度。然后统计kmer的种类和每一类kmer的数目，根据kmer的频数按照从高到底进行排序。选取频数最高的kmer作为seed开始向3端延伸，每次延伸一个碱基，统计延伸之后每种kmer出现的次数，选择kmer频数最高的那一条kmer作为延伸路径。如果在某一支路出现了两种kmer频数一样的情况，就分别延伸这两种kmer，看下一条路径的kmer频数，取 kmer频数较大的那一条路径，直到kmer不能再延伸。再从5端开始向5端继续进行这样的操作，直到不能延伸，这样就得到了一条线性的contig，我们将所使用的kmer从我们存储的容器中去掉。继续生成下一条contig，如此重复我们会得到各种长短不一的线性contig。
Chrysalis过程：这一步可以是对得到的线性contig的聚类。
1. 我们会对得到的contig根据长度进行过滤，去掉小于100bp的contig
inchworm.K25.L25.fa通过脚本过滤得到了inchworm.K25.L25.fa.min100
这里的覆盖率是id上的值，程序取的阈值是10 例如: >a12;15 15就是覆盖率

利用bowtie2将both.fa(所有的reads)比对到inchworm.K25.L25.fa.min100得到比对的bam文件
根据reads的比对结果和inchworm.K25.L25.fa得到初步聚类的contig文件iworm_scaffolds.txt

overlap的判断条件条件有两种情况：
- 1 先判断成对的reads的比对情况。如果成对的reads1和read2分别比对到了两条不同的contig1和contig2，那么这两条contig也是属于同一个类的，会输出这对contig和支持的reads数。
- 2 无论是read1还是read2，如果这条reads同时比对到contig1和contig2，除开reads的比对到不同contig的起始位置和终止位置不是包含关系这种情况，都会认为两条contig存在overlap，我们会输出这两条contig，并输出支持的reads数目。
- 需要注意的是第二种情况比较的是reads的起始位置和终止位置。通过阅读源码，发现没有对出现包含情况的contig进行输出，可能是因为下一步聚类时是通过contig的overlap情况来进行聚类，避免这里重复进行计算
通过GraphFromFasta将iworm_scaffolds.txt的结果进一步聚类得到iworm_cluster_welds_graph.txt

这里是根据overlap继续进行判断：这一步的输入文件为iworm_scaffolds.txt，线性的contig结果，both.fa文件。首先会读入iworm_scaffolds.txt存一个map，map的结构是key是读入的id，value是一个Pool类
利用BubbleUpClustering根据iworm_cluster_welds_graph.txt的结果得到聚类的contig文件

聚类本质是将有联系的contig放在一个pool也就是一个类里面进行输出，没有聚类的contig单独进行输出，得到一个fasta格式的文件。

如果我来做这一步可能会用比较笨的思路：就是从第一行开始与每一行分别求取交集，有交集就合并，没有就下一个循环，遍历完一遍之后再用新的集合再与剩下的元素进行一次遍历，如此重复一个集合与其他的集合都没有交集，就输出这个集合，然后从第二个集合开始重复进行相同的操作，直到最后一个集合。但是这样应该会耗费很多时间来进行遍历

利用CreateIwormFastaBundle将属于同一类的contig合并得到bundled_iworm_contigs.fasta

这一步比较简单，将属于一个类(pool)的序列用X连接在一起输出成一条序列
通过ReadsToTranscripts将both.fa比对到bundled_iworm_contigs.fasta得到readsToComponents.out同时进行排序
将分类的reads分别再运行一次Trinity，这一步是作者写了一个多线程的工具Parafly来运行，同时在对于每一类reads重新运行Trinity时会调用了Butterfly这个jar包，Butterfly这个jar包是利用德布鲁因图的算法来得到拼接完成的contig信息。最后将每一类的trinity结果进行合并得到最终的Trinity.fasta
- 这里在源码中作者使用了一个递归的方式，对每一类聚类的后的reads分别重新运行Trinity同时会有一个TRINITY_COMPLETE_FLAG的变量来控制是否要运行butterfly这一步

线性contig的格式如下：
GraphFromFasta过程的描述：
iworm_scaffolds.txt 文件格式:

读入iworm_scaffolds.txt 文件
- 读入第一行的时，存一个map类型的数据结构，map的键的类型为整型，值得类型为Pool类。以上图中文件第一行为例，map的键为228，值是一个Pool类，Pool类中含有479这个元素；然后再存一个同样类型的map 键为479 值为Pool类，Pool类中含有228这个元素。数据结构示意图如：228=>[479];479=>[228]
- 读入第二行，数据结构示意图如：0=>[480];480=>[0]
- 读入第三行,由于 228这个map已经存在，因此向228对应的值得Pool中添加了505这个元素，228的数据结构改变为如下：228=>[479,505]同时会得到一个键为505的map，示意图如：505=>[228]
- 遍历完文件得到一个map容器
读入contig序列文件

contig文件格式：

读入contig文件的时候，一边读入文件，一边将序列打断程kmer，kmer的初始长度为12。
kmer在保存的过程中会转化成一个int类型的数字，每一种kmer用一个int类型的数字表示，这样做可能是为了减少程序运行过程中内存的消耗
例如现在有序列：

>a1

GTGCTGTTTGGTT
这条序列就会保存为下面的数据格式 a1为id的编号，0为相对于序列的起始位置，实际运行过程中GTGCTGTTTGGT应该是用一个数字表示，数据结构类型也是map类型

GTGCTGTTTGGT => [(a1,0)]

TGCTGTTTGGTT => [(a1,1)]
读入第二条的时候

>a2

GTGCTGTTTGGTA
由于GTGCTGTTTGGT已经存在因此

GTGCTGTTTGGT => [(a1,0),(a2,0)]

TGCTGTTTGGTA => [(a2,2)]

这样读完所有的contig序列，就会得到不同的kmer，并且kmer可以代表有联系的contig

然后再遍历一遍contig序列，将contig打断成24kmer的长度
- 对24kmer分成两部分进行判断
  
  ------------————
- 假设这一条contig来源于contig1，长度是24kmer，红色的12kmer在前面存的数据中是属于contig1，contig2，contig4；黑色的12kmer在前面的结果中属于contig1,contig3,contig2;对于这种情况我们会输出 24kmer类型和对应的contig分类，格式和12kmer的存储格式一致。
- 如果没有设置reads支持的参数，那么我们就会将属于同一个24kmer的元素分别添加到各自的map中，如果不存在相应的map，但是又有聚类在一起，那么就建立一两个新的map，相互存到对应的Pool元素中。例如 4-> 6 6 -> 4;
- 如果设置有reads支持，将reads打断成48kmer长度，那么就将有overlap的contig 成对左右各延伸12bp变成一个48bp的序列，统计reads的支持情况，将reads数大于某一个阈值的成对的contig添加其对应的map中。

48bp序列的延伸方式如下

例如有三条contig属于同一个24kmer contig1，contig2，contig3


contig1：----------------`---------------`

                                   `--------------------`--------------- contig2

contig1：             `------------------`---------

                -----------`------------------`                contig2

contig1：------------`-------------------`

                              `--------------------`---------------   contig3

contig1：             `------------------`---------

                -----------`------------------`                contig3

contig2：------------`-------------------`

                              `--------------------`---------------   contig3

contig2：             `------------------`--------

                ----------`-------------------`                contig3

BubbleUpClustering过程：
对聚类的结果进行简化：

1 -> 6

2 -> 6

3 -> 6

4 -> 7

5 -> 7

6 -> 1 2 3 8

8 -> 6 7 9

7 -> 8 4 5

15 -> 16

16 -> 15

9 -> 8

首先读入文件，将文件存为一个向量，向量的元素为Pool类型，Pool的m_id为每一行的第一个元素

m_index的元素为 ->之后的元素。然后对这个向量进行排序，Pool的大小，由小到大进行排序。

排序后的向量如下：

1 -> 6

2 -> 6

3 -> 6

4 -> 7

5 -> 7

9 -> 8

15 -> 16

16 -> 15

8 -> 6 7 9

7 -> 8 4 5

6 -> 1 2 3 8

开始遍历这个向量：
- 建立一个map pool_idx_to_containment，map的key是1 value是一个Pool Pool的id是1，里面暂时没有元素
- 建立一个map pool_idx_to_vec_idx，map的key是1 value是对应的向量的下标值
- 遍历向量时，p是对元素的引用，第一个元素是Pool类，id为1，p是对Pool的引用，如果p的元素个数大于0那么遍历这个Pool类，如果p的元素个数不大于0就进行下一个循环。
- 这个Pool类的元素赋值给other_id = 6
- 在属于6的Pool中添加1
- pool_idx_to_containment[other_id].add(id)
- 在属于6的Pool中添加1对应的Pool中的元素
- pool_idx_to_containment[other_id].add(pool_idx_to_containment[id])
- 在属于6的Pool中去掉6
- pool_idx_to_containment[other_id].exclude(other_id)

在原始的6属于的向量Pool中去掉1
pool_vec[ pool_idx_to_vec_idx[ other_id ] ].exclude( id )
执行完这一步之后会跳出循环
重新对1的Pool中的元素进行遍历
adjacent_id = 6
adjacent_pool这个Pool是原始向量中的6对应的Pool
if (adjacent_id != other_id && adjacent_id != id)
如果7不等于6 也不等于8那么,此时adjacent_pool这个Pool是原始向量中的7对应的Pool,进行下面的判断
if (! adjacent_pool.contains( id )) 退出程序
adjacent_pool.exclude( id ) 在adjacent_pool这个Pool是原始向量的pool中去掉id元素
if (! adjacent_pool.contains( other_id )) adjacent_pool.add( other_id )
other_id_pool是other_id对应的原始Pool
if (! other_id_pool.contains( adjacent_id )) other_id_pool.add( adjacent_id )
遍历完成之后会清空目前的p属于的Pool

pool_idx_to_containment[id].clear()

然后进去遍历向量元素的下一个循环，直到向量元素遍历完成

正常情况下会对向量元素进行两边遍历，保证聚类的结果是独立的

然后对聚类后的结果进行输出，属于一个类别的就输出，保证同一个类别的序列长度不小于200，对于没有聚类的序列单独作为一类进行输出，序列长度也不小于200

reads进行分类 * 首先将contig打断成kmer，k的长度为25，这个时候会储存kmer属于哪一条序列的角标，然后

将reads打断成kmer，判断这条kmer属于那一条序列，并将这条序列的角标放入一个向量,统计向量里面出现次数最多的哪一个角标为最优的角标，出现的次数减1除以所有的kmer数加上0.5就是reads占序列的百分比

上一篇：Trinity的介绍

下一篇：hexo 主题模板的设置