生信教程：使用拓扑加权探索基因组进化（2）
动动发财的小手，点个赞吧！使用 Twisst 探索整个基因组的进化关系的拓扑加权[教程](https://github.com/simonhmartin/tutorials/blob/master/topology_weighting/README.md "Source")。
## 简介
拓扑加权是量化不一定是单系群之间关系的一种方法。它通过考虑更简单的“分类单元拓扑”并量化与每个分类单元拓扑匹配的子树的比例，提供了复杂谱系的摘要。我们用来计算权重的方法称为 Twisst：通过子树迭代采样进行拓扑权重。
在本次实践中，我们将使用模拟数据来探索拓扑权重如何提供谱系历史。然后，我们将尝试使用针对窄窗口推断的邻居连接树来推断整个模拟染色体的拓扑权重。
## 工作流程
在真实情况中，我们不知道真正的家谱历史，但我们有一组序列，我们希望从中推断家谱。我们将使用一种简单的方法来做到这一点：使用标准的系统发育工具为整个基因组的窗口进行系统发育。通过将我们推断的历史与 R 中的事实进行比较，我们将深入了解谱系推断中功率和分辨率之间的权衡。
## 从序列数据推断权重
上面我们使用了模拟的“真实”家谱。在大多数情况下，我们拥有的只是序列数据，并且必须推断其进化历史。事实上，有两件事我们不知道：
- 我们不知道所有个体之间的谱系关系
- 当我们沿着染色体移动时，我们不知道重组改变关系的“断点”
在这一部分中，我们将从序列数据开始。我们将使用一种相当简单的方法，在沿着基因组的窗口中推断家谱。然后我们将对这些进行扭曲，看看我们是否可以恢复一些接近基本事实的东西。
请注意，这一部分的教训之一是，在窗口中推断树是粗略的，并且可能存在缺陷，特别是如果我们弄错了树的大小。
### 下载代码和数据
这部分的脚本位于github上的genomics_general包中，我们需要下载：
```sh
#download package
wget https://github.com/simonhmartin/genomics_general/archive/v0.3.tar.gz
#extract files from zipped archive
tar -xzf v0.3.tar.gz
#delete zipped file
rm v0.3.tar.gz
```
为了确保Python可以识别这些库，请将“genomics_general”目录添加到Python路径中
```sh
export PYTHONPATH=$PYTHONPATH:genomics_general-0.3
```
我们将使用的序列文件随第 1 部分中下载的 twisst 包提供。该文件采用简单的 .geno 格式，其中包含每个个体的染色体、位置和基因型列：
```sh
zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.seqgen.SNP.geno.gz | head -n 5 | cut -f 1-8
```
### 推断树
我们有一个分布在 50 kb 基因组区域的 SNP 文件。我们将在定义数量的 SNP 的窗口中推断树，这样每个窗口都具有相似的信息量，但其在染色体上的绝对跨度可能不同，具体取决于 SNP 密度。
这种方法存在一个潜在的权衡。我们希望选择一个足够大的窗口大小，以便为树推理提供必要的能力，但又足够小，以实现足够的分辨率来捕获家谱历史在染色体上的变化。
我们将运行在 Windows 中读取 SNP 文件的脚本，然后使用 Phyml 为每个窗口推断一棵树。 Phyml 能够进行最大似然推断，但这里我们不会使用优化，因此树输出将是邻接树，使用 BIONJ 算法推断。通过模拟，我们发现邻接算法对于短序列的性能优于最大似然推理。
首先检查脚本的选项
```sh
python genomics_general-0.3/phylo/phyml_sliding_windows.py -h
```
主要选项是 -g 指定输入 .geno 文件和 --prefix 指定输出文件的前缀
还有用于设置窗口类型和大小的选项。我们将使用 20、50、100 和 500 个 SNP 等不同窗口大小运行该脚本四次。请注意，输入文件仅包含 SNP，因此通过设置 --windTypesites，每个窗口将被设置为包含固定数量的 SNP。通过这种方式划分窗口，它们在染色体上的绝对大小将随 SNP 密度而变化。每个窗口的开始和结束位置将记录在输出文件中。
最后，还有如何运行 Phyml 的选项。这里我们只需要设置 --optimise n 来关闭最大似然，并用 --model HYK85 定义替换模型，因为这是模拟序列的模型。
使用循环运行脚本，每次设置不同的窗口大小
```sh
for x in 20 50 100 500
do
echo "Inferring trees with window size $x"
python genomics_general-0.3/phylo/phyml_sliding_windows.py \
-g twisst-0.2/examples/msms_4of10_l50k_r500_sweep.seqgen.SNP.geno.gz \
--prefix msms_4of10_l50k_r500_sweep.seqgen.SNP.w$x.phyml_bionj \
--windType sites -w $x--model HKY85 --optimise n
done
```
每次运行时，脚本都会生成两个输出文件： .trees.gz 文件包含每个窗口的树，采用 Newick 格式。
我们现在可以使用树文件作为输入来计算染色体上的权重。
使用为每个不同窗口大小生成的树文件运行 twisst
```sh
for x in 20 50 100 500
do
echo "Running Twisst for window size $x"
python twisst-0.2/twisst.py \
-t msms_4of10_l50k_r500_sweep.seqgen.SNP.w$x.phyml_bionj.trees.gz \
-w msms_4of10_l50k_r500_sweep.seqgen.SNP.w$x.phyml_bionj.weights.tsv \
-g A 1,2,3,4,5,6,7,8,9,10 \
-g B 11,12,13,14,15,16,17,18,19,20 \
-g C 21,22,23,24,25,26,27,28,29,30 \
-g D 31,32,33,34,35,36,37,38,39,40 \
--outgroup D
done
```
### 绘制推断权重
我们现在可以绘制这些推断树在染色体上的权重。这可以在与之前相同的 R 脚本中完成。
再次打开 R（如果您已重新启动 R，则可能需要重新加载plot_twisst.R 脚本）。
```R
source("twisst-0.2/plot_twisst.R")
```
和之前一样，我们读入权重和窗口数据文件。这次我们将从模拟谱系中加载原始真实权重，以及我们刚刚计算的推断权重的四个文件。
```R
weights_files <- c('msms_4of10_l50k_r500_sweep.weights.tsv.gz', #true weights file from Part 1
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w20.phyml_bionj.weights.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w50.phyml_bionj.weights.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w100.phyml_bionj.weights.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w500.phyml_bionj.weights.tsv')
window_data_files <- c('twisst-0.2/examples/msms_4of10_l50k_r500_sweep.data.tsv.gz',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w20.phyml_bionj.data.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w50.phyml_bionj.data.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w100.phyml_bionj.data.tsv',
'msms_4of10_l50k_r500_sweep.seqgen.SNP.w500.phyml_bionj.data.tsv')
```
加载测量和窗口数据。请注意，当给定多个输入文件时，import.twisst 函数会将它们解释为单独的染色体。
```R
twisst_data <- import.twisst(weights_files, window_data_files,
names=c("True weights", "20 SNP windows", "50 SNP windows", "100 SNP windows", "500 SNP windows"))
```
现在我们再次绘图以比较真实权重和推断权重。 （您可能需要扩展绘图窗口才能正确显示多个绘图）。
```R
pdf("simulted_and_inferred_trees_weights_comparison.pdf", width=8, height=8)
plot.twisst(twisst_data, show_topos=FALSE, include_region_names=T)
dev.off()
```
这表明两种拓扑占主导地位。一种是拓扑 3，您会注意到它是“物种”拓扑，其中 cydno (chi) 与 timareta (txn) 组合在一起，并且两个 melpomene 种群（ros 和 ama）组合在一起。另一种是拓扑 11（粉色），其中 txn 与 ama 分组（嵌套在 melpomene 对内）。这告诉我们，基因渗入的主要方向是从 H. melpomene amaryllis 进入 H. timareta thelxinoe。这符合我们目前对这个群体的理解，在这个群体中，提马雷塔沿着安第斯山脉的东坡扩张，与墨尔波烯杂交，并获得了每个地区偏爱的警告模式等位基因。
未完待续!
页:
[1]