覆盖率
对已知的21条序列中挑选若干条进行搜索，基本都可以根据名称定位TaxId。
| Annotation | organism | TaxId |
|---|
| >BAA03443.2 1,3,4,6-tetrachloro-1,4-cyclohexadiene hydrolase [Sphingobium indicum] | Sphingobium indicum | 33205 | | >WP_003413363.1 MULTISPECIES: haloalkane dehalogenase [Mycobacterium tuberculosis complex] | Mycobacterium tuberculosis complex | 77643 | | >AAL17946.1 putative linB haloalkane dehalogenase, partial [Mycolicibacterium smegmatis MC2 155] | Mycolicibacterium smegmatis MC2 155 | 246196 | | >P0A3G2.1 RecName: Full=Haloalkane dehalogenase [Rhodococcus rhodochrous] | Rhodococcus rhodochrous | 1829 | | >WP_007349233.1 haloalkane dehalogenase [Marinobacter sp. ELB17] | Marinobacter sp. ELB17 | 270374 | | >3U1T_A Haloalkane Dehalogenase, DmmA, of marine microbial origin [unidentified] | unidentified | 32644 | | >sp|P22643.2|DHLA_XANAU RecName: Full=Haloalkane dehalogenase | | |
某些序列可能没有 [] 中定义的来源生物。
在某些序列中可能有多个 []，例如0_blast.txt中的三个例子：
PHQ80941.1 haloalkane dehalogenase [Coxiella sp.] [Coxiella sp. (in: Bacteria)]
SRX93098.1 haloalkane dehalogenase [Amycolatopsis mediterranei S699] [Mycobacterium shimoidei]
EAT11207.1 haloalkane dehalogenase [Oceanobacter sp. RED65] [Bermanella marisrubri]
因此，根据名称定位 TaxID 的任务基本可以实现。
广义来说，序列定位与专利搜索定位相同：均是层次分类，将叶节点挂到层次分类的数据上去。专利的IPC号提供了完整的层次分类路径，通过Excel函数提取各层次，辅助以数据透视表，可以得到各种统计图表。在这里，rankedlineage提供了各层次数据，因此，也可以最终借助Excel中的数据透视表，给出各种统计图表。
层次分类数据的Map，应该由专门的非关系型数据库来实现。可以去查阅python, Neo4j等对于此类任务的常规做法。
任务
前期探索又打开了三个支线任务：
- 搜索资料，编程实现微生物数据对序列的引入（预计1天）
- 搜索文献，例如祖先重构、微生物谱系与序列的关系等，了解序列引入微生物数据的作用。（预计2-3天）
- 继续执行序列筛选任务，主要是Jalview的使用。（预计2-3天）
利用Jalview实现各种序列筛选功能仍然是先决条件，因此今天先进行JalView的学习，之后任务2，再任务1。
JalView(预计2-3天)
希望实现序列的 1. 导入 2. 标记 3. 编辑 4. 导出满足某种筛选条件的序列 5. 编辑过程可追溯。
MSA
上传 0_out.txt 文件进行 MSA。该文件中有14000余条序列，大小为5165 KB，超过了Clustal Omega的要求。
This tool can align up to 4000 sequences or a maximum file size of 4 MB.
https://www.ebi.ac.uk/Tools/msa/
EBI 提供了 Clustal Omega, EMBOSS Cons, Kalign 等各种 MSA 工具，但是其网络提供的 MSA 大小都有限制，反而是Clustal Omega的4000条序列是最大的。
这不仅让人反思，是否对一万多条序列进行MSA是没有意义的。这启示我们在使用工具时考虑问题规模以及流程。
其次是其他的在线或软件工具能否实现这一点？
序列聚类与筛选
本质上希望实现的功能是基于已有的21条序列，通过聚类算法找到与它们属于同一类的序列。
目前有以下方案：
| 来源 | 相似性计算 | 聚类方法 | 算法 | 超参数 | 状态 |
|---|
| 论文 | USEARCH 计算 pairwise 相似性 | UMPGA | 层次聚类 | 簇的数量 | 聚类软件有bug，难以调试 | | 脑补 | USEARCH 计算 pairwise 相似性 | biopython | 层次聚类 | 簇的数量 | 没有可视化，算法效率未测试 | | 脑补 | Clustal Omega 计算 MSA | JalView | 层次聚类 | 簇的数量(阈值) | MSA计算有规模限制（4000条，4MB） | | 搜索 | USEARCH 计算 pairwise 相似性 | USEARCH | (KNN)聚类 | 簇的数量 | 有学习成本 | | 脑补 | USEARCH 计算 pairwise 相似性 | scikit-learn | kNN聚类，层次聚类，谱聚类 | 簇的数量 | 处于脑补阶段 |
下面就MSA的大规模计算和USEARCH直接进行聚类进行探究，二者均基于USEARCH。
MSA的大规模计算
http://www.clustal.org/omega/ 10698 Clustal Omega
https://www.drive5.com/muscle/ 35122 MUSCLE
MUSCLE的引用量相对多，先看看MUSCLE。
资料
Multiple Sequence Alignment -
Bioinformatics Algorithms: An Active Learning Approach
上述视频讲了MSA的算法，多个序列比对的曼哈顿距离算法是O(2knk)，而通过profile方法进行 alignment，可以有效降低复杂度，得到的局部最优解大部分情况下在现实中可用。
（也解了惑，profile即对每个位置的每种残基/核苷酸出现概率的统计。）
https://www.bioinformaticsalgorithms.org/
下午看了第五章，算是了解了序列比对算法是如何做的。对于序列比对感兴趣的推荐看一看视频，共10篇，总共约2小时。
- 序列比对具有重要生物学价值；
- 序列比对可以形式化为最长公共子序列的寻找；序列比对最重要的概念就是曼哈顿距离(一个有向图)，一定要有这个图象，之后的各种变体都是在该图上进行的；
- recursive算法效率过低，通过动态规划法（查表），时间复杂度与边的个数成正比，空间复杂度与点的个数成正比，对于pairwise，为序列长度的二次方；
- Back tracking 算法可以得到最优匹配；
- Global 和 Local 的区别在于后者在起点/终点到每一个中间节点间加上了两族免费的Taxi；
- 边的价值通过评估矩阵来打分，对于蛋白有特殊的评估矩阵（BLOSUM, PAM）；
- 加入delete和insert罚分可以通过3-level的Manhattan来实现；
- 为了满足长序列的空间复杂度要求，可以时间换空间，不得到最优匹配时以线性空间和二次方时间得到最优评分，需要最优匹配时以线性空间和二倍的二次方时间得到最优评分和最优匹配；
- 多序列匹配变为高维空间中的Manhattan，时间和空间复杂度爆炸，因此需要启发式的方法，通过greedy方法 + profile，可以在有意义的时间内完成比对。
Multiple Sequence Alignment using MUSCLE
BioPandit is an open platform to learn bioinformatics through video lectures and hands-on tutorial in Youtube.
解释了MUSCLE最简单的用法：
muscle.exe -in <inputfile> -out <outputfile> [-clw]
默认的输出格式是Fasta, -clw 可以输入clustalW 格式的结果。
用了 UGENE 来做序列的可视化，与JalView看起来类似。
MUSCLE
在本机(i5-8265U,4核)上进行测试：
| 序列数 | 最长序列 | 平均长度 | Iter | 内存(MB) | 时间 |
|---|
| 30 | 302 | 295 | 20 | 14 | 0:00 | | 295 | 386 | 295 | 20 | 63 | 0:17 | | 876 | 765 | 295 | 23 | 339 | 4:07 | | 2029 | 765 | 295 | - | - | - |
到2000条已经很吃力了。可以将序列的长度缩减一下，实现一下python里面对于序列的裁剪。
在台式机(i5-10400F,6核)上进行测试：
| 序列数 | 最长序列 | 平均长度 | Iter | 内存(MB) | 时间 |
|---|
| 30 | 302 | 295 | 20 | 14 | 0:00 | | 295 | 386 | 295 | 20 | 63 | 0:13 | | 876 | 765 | 295 | 23 | 339 | 3:28 | | 2029 | 765 | 295 | 24 | 1417 | 33:24 | | 741 | 398 | 295 | 22 | 250 | 1:20 | | 1974 | 398 | 295 | 24 | 1344 | 16:22 |
表明PC对于大规模的比对(上千万条)是无能为力的，还是要用profile-profile的方法。而且可能更有意义。
profile-profile alignment
目前有
- 序列截断
Each sequence identified was then shortened so that it did not exceed any query by more than 20 amino acids at either end.
是否有现成的程序可以做，还是需要自己实现一个朴素的想法？比如根据截断后的pairwise-alignment 的值最大值来选。比如确定query上核心的domain，然后向外延申一定长度截断。
- cut-off 确认
现在再看有感觉了，将已知的序列分到4个cluster中，然后确定阈值-数量最稳定的那个阈值，得到的各个簇的序列即为最终输出的序列。
The tree was cut at varying cut-offs, and the cut-offs at which all 22 known
HLD sequences (Table S1) were grouped into four clusters, i.e. subfamilies HLD-I, HLD-II and HLD-III and a cluster consisting of the putative HLD from A. niger (“HLD-IIIb”), were further analyzed. For each such cut-off c, the difference between the total number of sequences comprising HLD clusters at the cut-off c+1 and the number at the given cut-off c was calculated. The cut-off value of 74, i.e. the highest cut-off value among all those values giving the smallest difference in the overall sequence count, was selected as the final cut-off for cutting the hierarchical sequence tree.
- Multiple Sequence Alignment
首先对4个簇进行alignment，然后进行profile-profile alignment。
First, each of the four HLD clusters was aligned individually, then profile-profile alignment of the two closest clusters was constructed in a stepwise manner until all clusters were aligned.
考虑到最终只有不到1000条序列，这种MSA的方法同时考虑了簇内和簇间，可能能得到更准确的结果。
USEARCH
http://www.drive5.com/usearch/manual/cmds_by_category.html
有以下几个目录可能是相关的。
Sequence alignment commands
Sequence, tree and graph-based clustering
Distance matrices
|