用于 IQtree 的序列串联方法

---

前言

因为这一段时间很懒而且现实生活的事情又多，所以系统发育分析的正式文章教程可能要搁置一会，但是想着总得更新些东西出来，所以就以这篇文章作为正餐开始前的一叠小菜分享分享。

事实上序列串联对于刚涉及到这一块的新手来说确实也会有些迷糊，比如说串联是怎么个串联法，需要准备什么文件，文件的内容格式应该是怎样等等问题都会存在，这篇文里就针对这一点讲一讲。

更新日志

2023-04-07 简化了脚本，并完善了脚本的注释。

2023-07-04 完善了脚本，并添加了部分内容。

需要知道的知识

一些系统发育分析中的概念是最基础的，如果未作了解或了解不深可以看一下这篇文章：生命之树及其应用 。写的非常全面和具体，里面涉及到的各种术语也很专业，同时也介绍了很多常用的软件等。本文涉及到的地方就是这篇文章展望部分中第一点讲的 “超大树构建方法的革新” 中的超矩阵方法。

序列串联方法

需要准备的文件

首先按照惯例推出几个能够实现的软件：Phylosuite、FASconCAT-G。前者用起来非常方便但是我只下了 windows 的版本，虽然它有 linux 的版本但我也没去调查怎么用，索性就自己写了个脚本进行序列串联。

先说明需要的文件内容格式：

>tax1@seq1
ABCDEFG
>tax2@seq1
ABCDEEG
>tax3@seq1
ACCDEEG

这是一个由不同物种的同源基因组成的多序列比对，其中不同 tax 代表不同的物种，以 @ 作为分隔符，其后对应不同物种的同源序列名称。对应的序列是已经经过比对修剪的序列，允许 gap 存在。每个物种至多存在一条序列。例（我随便编的）：

这是一个同源序列（蛋白）：

>Human@human_1215
AAABKISPLM
>Mouse@mouse_98455
AAXBKISPLM
>Fish@fish_1561
AAABKISPNM

这是另一个同源序列（蛋白）：

>Human@human_5372
CJBSODNONVOIAA
>Mouse@mouse_12389
CJBSODNNNVOIAS
>Fish@fish_1001
DJBSODNBNVOIAS

串联的意思是把它们串成这样（一个比对中不一定包含所有物种，串联以后该比对不存在的物种会全部以 gap 代替）：

>Human
AAABKISPLMCJBSODNONVOIAA
>Mouse
AAXBKISPLMCJBSODNNNVOIAS
>Fish
AAABKISPNMDJBSODNBNVOIAS

串成这样以后，也要有对应的 partition 信息以供软件进行超矩阵算法。

对我而言，实现上述操作的方法有两种：

1. ①、简化 ID 成物种编号（就是把 @ 以及后面的序列编号全删了）。

   ②、将处理后的序列输入 Phylosuite 处理。

2. 编写脚本一步到位。

对于第一种方法如果序列比较少是没什么问题，因为下载到 windows 不用多久（当然如果会 phylosuite 的 Linux 端操作也会变得非常方便），但是如果涉及到的序列很多且很大那就只能自己动手丰衣足食了。

脚本

2023-07-04 引流补充：

如果是 Windows 系统：不嫌麻烦的话可以下载我基于 PyQt5 编写的 JuseKit，其中包含序列串联功能且输出信息更完善。

如果是 Linux 系统，以下是我编写的序列串联脚本：

# -*- coding: utf-8 -*-
# sequence_con
# Verson: 1.1
# Date: 2023.1.14
# This script concatenates sequences and generates an IQ-TREE2 partition file.

import os

IQpartition = '#nexus\nbegin sets;\n'
total_length = 0
total_spe_seq = {}


total_num = len(os.listdir('./'))
compl_num = 0
aligns = os.listdir('./')
aligns.sort()

#对当前文件夹里所有文件逐个读取
for align in aligns:

	with open(align, 'r') as ali:

		#初始化当前文件信息
		tmp_len = 0
		tmp_species = []
		for line in ali:
			if line.startswith('>') and line.strip() != '>':
				#如果物种名以下划线分割则替换为_
				spe_name = line.split('@')[0].strip()

				tmp_species.append(spe_name)

				#到当前比对才出现的物种之前的序列都以gap表示
				if spe_name not in total_spe_seq.keys():
					total_spe_seq[spe_name] = '-'*total_length

			#跳过空序列
			elif line.strip() == '>' or line.strip() == '':
				continue

			#将第一个序列的长度读取以作为该多序列比对的长度
			else:
				if len(tmp_species) == 1:
					tmp_len += len(line.strip())

				total_spe_seq[spe_name] += line.strip('>').strip()

		#对不在该文件出现但在之前出现过的物种以gap作为其序列进行补充
		for spe in total_spe_seq.keys():
			if spe not in tmp_species:
				total_spe_seq[spe] += '-'*tmp_len
		#添加分区信息
		charsetid = align.split('.')[0]
		IQpartition += f"\tcharset {charsetid}={total_length+1}-{total_length+tmp_len};\n"

		total_length += tmp_len
	compl_num += 1
	print(f"The progress is {compl_num}/{total_num}...")


os.mkdir("con_res")
#生成物种log
with open("con_res/sequence_con.log", 'w') as log:

	log.write("Species\tLength\n")
	for spe in total_spe_seq.keys():

		if len(total_spe_seq[spe]) == total_length:
			log.write(f"{spe[1:]}\t{len(total_spe_seq[spe])}AA\t+\n")
		else:
			log.write(f"{spe[1:]}\t{len(total_spe_seq[spe])}AA\t-\n")

	print("See running log in con_res/sequence_con.log")
#生成串联文件
with open("con_res/concatenation_ortho.fasta", 'w') as f:

	for spe in total_spe_seq:
		f.write(spe + '\n' + total_spe_seq[spe] + '\n')
#生成分区文件
with open("con_res/IQ_partition.txt", 'w') as f:

	f.write(IQpartition)
	f.write("end;")


print('Finished, see output at con_res.')

原理比较简单，不过相对于其他软件的处理而言，这个脚本还可以解决 “空序列” 的问题。

运行也很简单，在存放多序列比对文件的文件夹中输入 python /path/sequence_con.py 即可（脚本不要放里面，也不要存在其他类型的文件，因为我没有加尾缀识别的判断语句）。

运行结束后会多出一个 con_res 的文件夹，里面有合并后的序列 concatenation_ortho.fasta 以及用于 IQtree 建树的分区信息（IQ_partition.txt），此外也有一个运行日志 sequence_con.log，里面记录了每个物种的序列长度信息，最后一列为 + 则表示与第一个长度一致，理应上所有的物种都应该为 + ，如果出现 - 号则说明该物种的某些比对存在错误，长度与其他物种不一致。

后记

这一段时间完全放飞自我，正事没咋干摸鱼摸一堆，一方面是身不由己另一方面确实打不起什么劲，说明环境对于人的影响还是非常显著的。

相关的 python 脚本同以前一样放在了 github 里，感兴趣可以自取，若是不嫌麻烦事实上 phylosuite 也是个不错的选择（国人软件，能多引一个是一个）。