1. 引言
类群的适应性研究一直是系统发育研究的下游，同时也是将传统的分类学问题升华到物种多样性的重要桥梁。但是，传统的适应性演化大部分都停留在通过人为的将系统发育树或者多样化速率与一些气候因子数据或者海拔相互联系，很难进行量化。今天来介绍一个新的软件aPhyloGeo，这个软件主要的算法就是通过将基于气候因子构建的聚类树和基于遗传信息（基因）构建的基因树的拓扑结构进行比较，从而找到和特定气候因子相关联的基因，也就确定了那些基因驱动了类群的适应性演化。
2. 软件下载
个人看了半个小时的官方下载教程，并且在Gemini的帮助下，才搞明白这个软件的安装和使用逻辑，真的很反人类。逻辑大概就是要通过poetry这个软件，进入作者配置在安装路径中的一个环境，在这个环境完成后续的分析。
git clone https://github.com/tahiri-lab/aPhyloGeo.git
cd aPhyloGeo
pip install poetry
poetry install
poetry self add poetry-plugin-shell
# 进入环境
poetry shell
3 软件使用
在引言里面说了，这个软件其实和根据基因构建的遗传信息树和环境因子构建的聚类树进行结构相似性比较，从而得出哪个环境因子和那些基因或者那一段基因具有相关性，所以第一步就是准备每个样本（物种）采样地点的环境信息，例如平均气温，降水量等等。这时候可能会有出现问题了，谁采样的时候会记录什么气温，降水量啊？确实，如果不是专门做生态或者环境调查的不会记录这些。但是没关系，我们有采样点和经纬度啊，可以基于WorldClim的19个生物因子的数据，得到我们每个样本的采样点对应的一些环境因子的数据，这样就解决了我们分析的问题，因此第一步就是根据采样点的经纬度和19个气候因子，得到每个采样点的环境因子数据，至于怎么下载19个气候因子，请自行网上查找学习。
准备一个csv格式，样品的经纬度和海拔数据（海拔可以没有，是可选的）
| ID | Longitude | Latitude | Altitude |
|---|---|---|---|
| Sample_01 | 114.3 | 30.6 | 20 |
| Sample_02 | 116.4 | 39.9 | 50 |
| Sample_03 | 121.47 | 31.23 | 10 |
经纬度转化为环境因子，输出的文件就是aPhyloGeo的其中的一个输入文件
python ~/Gemini/Location2enviroment_1.py -c coordinates.csv -t wc_data -o test.csv
head test.csv
ID,Precipitation_of_Driest_Quarter,Mean_Temperature_of_Warmest_Quarter,Max_Temperature_of_Warmest_Month,Mean_Temperature_of_Wettest_Quarter,Precipitation_of_Wettest_Quarter,Temperature_Annual_Range,Mean_Temperature_of_Coldest_Quarter,Precipitation_of_Driest_Month,Min_Temperature_of_Coldest_Month,Isothermality,Precipitation_Warmest_Quarter,Annual_Mean_Temperature,Precipitation_Seasonality,Precipitation_of_Wettest_Month,Annual_Precipitation,Mean_Diurnal_Range,Temperature_Seasonality,Mean_Temperature_of_Driest_Quarter,Precipitation_of_Coldest_Quarter
Sample_01,129.0,28.108,32.996,25.945333,561.0,32.272,5.813,31.0,0.724,25.054485,510.0,17.257708,54.785477,203.0,1265.0,8.085583,904.89557,7.8065,129.0
Sample_02,14.0,25.4365,31.328,25.4365,440.0,40.255997,-1.7836666,4.0,-8.9279995,28.161226,440.0,12.493875,136.02074,206.0,578.0,11.336583,1107.1929,-1.7836666,14.0
Sample_03,144.0,26.649,31.302,26.614,442.0,29.947,5.8485003,48.0,1.355,23.248016,409.0,16.459126,48.190933,162.0,1061.0,6.9620833,853.7569,13.0485,153.0
接下来就是运行aPhyloGeo，我们已经有了其中一个输入数据，另外一个就是MAS。这里需要注意的是，这个软件的输入文件和输出文件以及一些参数是通过配置文件设置的，并且运行也比较奇怪，不需要指定配置文件的位置，软件会自己寻找，所以配置文件的位置是不能改的，配置文件的位置默认就在软件下载目录的/aphylogeo/params.yaml
file_name: './datasets/example/geo.csv' # 修改为刚生成的环境因子文件
specimen: 'id' #"Please enter the name of the column containing the specimens names: "
dist_threshold: 60
window_size: 200
step_size: 100
bootstrap_threshold: 100 #The number of bootstrap samples to generate.
reference_gene_dir: './datasets' # MSA文件路径
reference_gene_file: 'sequences.fasta' # MSA文件名
makeDebugFiles: True
alignment_method: '1' #Please select one ~ 0:no alignment, 1:pairwiseAligner, 2:MUSCLE, 3:CLUSTALW, 4:MAFFT
distance_method: '0' #Please select one ~ 0: All distance methods, 1: Least-Square distance, 2: Robinson-Foulds distance, 3: Euclidean distance (DendroPY)
fit_method: '1' #Please select one ~ 1:Wider Fit by elongating with Gap (starAlignment), 2:Narrow-fit prevent elongation with gap when possible
tree_type: '1' #Please select one ~ 1: BioPython consensus tree, 2: FastTree application
rate_similarity: 90
method_similarity: '1' #Please select one :
# 1: Hamming distance
# 2: Levenshtein distance
# 3: Damerau-Levenshtein distance
# 4: Jaro similarity
# 5: Jaro-Winkler similarity
# 6: Smith–Waterman similarity
# 7: Jaccard similarity
# 8: Sørensen-Dice similarity
preprocessing_genetic: 1
preprocessing_climatic: 1
preprocessing_threshold_genetic: 0.2 # proportion of gaps allowed
preprocessing_threshold_climatic: 0.7 #variance threshold for numeric columns
permutations_mantel_test: 999 #permutations for Mantel test
permutations_protest: 999 #permutations for PROTEST analysis
mantel_test_method: "pearson" #method for Mantel test
statistical_test: '1' #Please select one ~ 0: Both test, 1: Mantel test, 2: PROTEST (Procrustes) analysis
运行
make
# head /result/output.csv
Gene,Phylogeographic tree,Name of species,Position in ASM,Bootstrap mean,Least-Square Distance,Robinson-Foulds Distance,Normalised RF Dist,Euclidean Distance,window,r (Correlation coeffici>./datasets/example/sequences.fasta,ALLSKY_SFC_SW_DWN,OL989074,200_399,100.0,3.41,4.0,1.0,0.71,,,,
./datasets/example/sequences.fasta,T2M,OL989074,200_399,100.0,3.89,4.0,1.0,0.62,,,,
./datasets/example/sequences.fasta,QV2M,OL989074,200_399,100.0,4.54,4.0,1.0,0.67,,,,
./datasets/example/sequences.fasta,WS10M,OL989074,200_399,100.0,4.1,4.0,1.0,0.53,,,,
./datasets/example/sequences.fasta,ALLSKY_SFC_SW_DWN,OL989074,400_599,100.0,3.4,4.0,1.0,0.71,,,,
./datasets/example/sequences.fasta,T2M,OL989074,400_599,100.0,3.88,4.0,1.0,0.62,,,,
./datasets/example/sequences.fasta,QV2M,OL989074,400_599,100.0,4.53,2.0,0.5,0.67,,,,
./datasets/example/sequences.fasta,WS10M,OL989074,400_599,100.0,4.08,4.0,1.0,0.53,,,,
./datasets/example/sequences.fasta,ALLSKY_SFC_SW_DWN,OL989074,600_799,100.0,3.3,4.0,1.0,0.71,,,,
./datasets/example/sequences.fasta,T2M,OL989074,600_799,100.0,3.72,4.0,1.0,0.62,,,,
./datasets/example/sequences.fasta,QV2M,OL989074,600_799,100.0,4.35,4.0,1.0,0.65,,,,
./datasets/example/sequences.fasta,WS10M,OL989074,600_799,100.0,3.91,4.0,1.0,0.52,,,,
./datasets/example/sequences.fasta,ALLSKY_SFC_SW_DWN,OL989074,800_999,100.0,3.4,4.0,1.0,0.71,,,,
./datasets/example/sequences.fasta,T2M,OL989074,800_999,100.0,3.88,4.0,1.0,0.62,,,,
./datasets/example/sequences.fasta,QV2M,OL989074,800_999,100.0,4.53,4.0,1.0,0.67,,,,
./datasets/example/sequences.fasta,WS10M,OL989074,800_999,100.0,4.08,4.0,1.0,0.53,,,,
./datasets/example/sequences.fasta,ALLSKY_SFC_SW_DWN,OL989074,1000_1199,100.0,3.41,4.0,1.0,0.71,,,,
./datasets/example/sequences.fasta,T2M,OL989074,1000_1199,100.0,3.89,4.0,1.0,0.62,,,,
./datasets/example/sequences.fasta,QV2M,OL989074,1000_1199,100.0,4.54,4.0,1.0,0.67,,,,
./datasets/example/sequences.fasta,WS10M,OL989074,1000_1199,100.0,4.1,4.0,1.0,0.53,,,,
,,,,,,,,,Mantel test result:,-0.4671097202423033,0.27,5.0
随后会在当前文件夹下生成一个result文件夹，里面就是分析的接过，主要看output.csv就可以。个人认为比较重要的就是Normalised RF Dist（衡量基因树和环境因子树的拓扑结构相似性指标）、Mantel test result（衡量环境对该类群分化的一个作用）。
| Gene | 输入的序列文件路径（./datasets/example/sequences.fasta） |
| Phylogeographic tree | 环境/气候变量名称，用于构建地理树 |
| Name of species | 物种/序列标识符（如 OL989074，这是 GenBank 登录号） |
| Position in ASM | 序列滑动窗口位置，格式为 起始_终止（如 200_399 表示第200-399bp） |
| Bootstrap mean | Bootstrap 平均值，衡量系统发育树分支的置信度（100.0 表示完全支持） |
| Least-Square Distance | 最小二乘距离，衡量系统发育树与地理树之间的差异 |
| Robinson-Foulds Distance | Robinson-Foulds 距离，比较两棵树拓扑结构差异的经典指标 |
| Normalised RF Dist | 标准化 RF 距离（0-1之间，1.0 表示树结构完全不同） |
| Euclidean Distance | 欧氏距离，衡量两棵树之间的距离 |
| window | 窗口信息（此处为空） |
| r (Correlation coefficient) | Mantel 检验相关系数（仅最后一行有值：-0.467） |
| p (Significance level) | Mantel 检验显著性水平（0.27，>0.05 表示不显著） |
| n (Number of observations) | 观测数量（5个样本）
4. 结语
由于一些不可抗力的原因，这个博客写了三天才写完。我觉得这个软件在生物多样性研究领域是有很高潜力的，比如研究某个地理部分很有特点的类群，就可以知道那些环境因素是真正推动了物种的多样化，而不是在讨论里面意淫。另外，软件的本身其实有一些局限性，在分析的过程中，是利用单个基因采用滑动窗口的形式进行构建多个系统发育，从而衡量每个窗口的的基因树和环境因子树的拓扑结构相似性，从而进行判断相关性。但是在大多数的情况下，都是多基因建树，少则几十，多则几百上千，和软件作者的的设想还是不太一样。个人认为如果在真实的研究中，有两种策略（1）我们可以把多个基因串联起来，然后滑动窗口，最后筛选出和环境因子具有相关性的片段，并进一步找出这些片段来自于那些基因；（2）每个基因单独进行一次分析，最后对结果进行汇总。