生信学习笔记：使用SNP data做基因渗入分析 (2)

下载测试数据包

根据教程的提示，下载测试数据，还有对应的脚本来计算统计量D和Fd。

wget https://raw.githubusercontent.com/mmatschiner/tutorials/master/analysis_of_introgression_with_snp_data/data/NC_031969.f5.sub1.vcf.gz

git clone https://github.com/simonhmartin/genomics_general.git

为了确保Python能够识别此存储库中包含的库，使用以下命令将存储库的路径添加到常规Python路径：

export PYTHONPATH=/scratch/pawsey0149/hhu/test/introgress/new/genomics_general

可以通过以下命令测试python能否识别存储库中对应的python包：

python -c 'import genomics'

如果没错，就不会有任何报错信息。

测试数据预处理

由于该教程所提供的脚本都是使用更加轻量型的“基因型（.geno）”文件格式，首先我们使用脚本parseVCF.py，将VCF文件NC_031969.f5.sub1.vcf.gz转换为基因型文件格式。

python genomics_general/VCF_processing/parseVCF.py -i NC_031969.f5.sub1.vcf.gz | gzip > NC_031969.f5.sub1.geno.gz

简单查看一下VCF格式和基因型格式的异同，你可以发现，基因型格式的文件同样包含VCF格式中所有样本相同基因型信息，以及染色体ID和其位置，不同的是少了一些该位点的比对质量等信息。

###zcat NC_031969.f5.sub1.vcf.gz |tail -n3

NC_031969       38037706        .       A       G       62186.6 PASS    AC=2;AF=0.239;AN=4;DP=4172;ExcessHet=-0;FS=0;InbreedingCoeff=0.8559;MLEAC=214;MLEAF=0.246;MQ=44.24;QD=32.26;SOR=0.96      GT      ./.     1/1     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.       ./.     0/0     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.
NC_031969       38037812        .       C       T       4662.43 PASS    AC=4;AF=0.021;AN=30;BaseQRankSum=-0.366;ClippingRankSum=0;DP=4333;ExcessHet=0;FS=0;InbreedingCoeff=0.4767;MLEAC=17;MLEAF=0.019;MQ=34.14;MQRankSum=0.967;QD=30.21;ReadPosRankSum=-0.366;SOR=1.179  GT      0/0     0/0     ./.     ./.     ./.       ./.     0/0     0/0     ./.     ./.     0/0     ./.     0/0     0/0     0/0     ./.     1/1     1/1     ./.     ./.     0/0     0/0     ./.     0/0       0/0     0/0     ./.     ./.
NC_031969       38038142        .       A       G       71277.4 PASS    AC=2;AF=0.523;AN=2;BaseQRankSum=1.65;ClippingRankSum=0;DP=4119;ExcessHet=-0;FS=0;InbreedingCoeff=0.8646;MLEAC=422;MLEAF=0.533;MQ=54.55;MQRankSum=1.65;QD=27.96;ReadPosRankSum=1.65;SOR=0.418      GT      1/1     ./.     ./.     ./.     ./.       ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.     ./.       ./.     ./.     ./.     ./.

###zcat NC_031969.f5.sub1.geno.gz |tail -n 3

NC_031969       38037706        N/N     G/G     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     A/A       N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N
NC_031969       38037812        C/C     C/C     N/N     N/N     N/N     N/N     C/C     C/C     N/N     N/N     C/C     N/N     C/C     C/C     C/C     N/N       T/T     T/T     N/N     N/N     C/C     C/C     N/N     C/C     C/C     C/C     N/N     N/N
NC_031969       38038142        G/G     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N       N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N     N/N 

要将样本分配给物种以计算D统计量，我们需要编写一个包含两列的表，其中第一列包含样本ID，第二列包含相应的物种ID。尽管我们在ABBA-BABA测试的第一次应用中仅使用了四种物种，但你已经可以将该数据集的所有样本和种类包括在此表中;这将有助于后来与其它四种物种的分析。因此，将以下文本写入名为的新文件samples.txt：

###vi samples.txt

  IZA1  astbur
  IZC5  astbur
  AUE7  altfas
  AXD5  altfas
  JBD5  telvit
  JBD6  telvit
  JUH9  neobri
  JUI1  neobri
  LJC9  neocan
  LJD1  neocan
  KHA7  neochi
  KHA9  neochi
  IVE8  neocra
  IVF1  neocra
  JWH1  neogra
  JWH2  neogra
  JWG8  neohel
  JWG9  neohel
  JWH3  neomar
  JWH4  neomar
  JWH5  neooli
  JWH6  neooli
  ISA6  neopul
  ISB3  neopul
  ISA8  neosav
  IYA4  neosav
  KFD2  neowal
  KFD4  neowal

为了便于以后用其他四种物种进行分析，我们将使用四种物种ID的变量。我们将首先关注Neolamprologus cancellatus（“neocan”）的潜在杂交性质，这意味着我们将测试该物种的渐渗，因此我们将其用作“spc2”。原因是它已经被过往的研究推测为一个混杂的物种，接着我们把饰圈沼丽鱼（“telvit”）作为渗入的潜在的捐赠者，“SPC3”的位置。最后将（“astbur”）该数据集中包含最远的物种看作一个外群。你可以决定使用哪种物种作为第二种可能有助于杂交（发生基因渗透）的物种Neolamprologus或Altolamprologus fasciatus（“altfas”）属中的任何一种(上面列表所提到的），作为"spc1"。然后，将下面命令中的“XXXXXX”替换为您为“spc1”选择的物种的ID，然后执行该命令。

spc1=XXXXXX; spc2=neocan; spc3=telvit; spc4=astbur

然后，使用freq.py存储库中的脚本，计算四种物种的等位基因频率：

python genomics_general/freq.py -g NC_031969.f5.sub1.geno.gz -p ${spc1} -p ${spc2} -p ${spc3} -p ${spc4} --popsFile samples.txt --target derived | grep -v nan | gzip > NC_031969.f5.sub1.pops1.tsv.gz

在上面的命令中，指定了基因型格式的输入文件，-g NC_031969.f5.sub1.geno.gz用-p选项分别指定了要分析的四个物种，并提供了含有样本名称的文件--popsFile samples.txt。该选项--target derived规定，对于四种物种ID中的最后一种（此处为“spc4”）观察到的等位基因应被视为祖先等位基因，因此备选等位基因是衍生的等位基因。忽略第四物种中等位基因未固定的位点。使用| grep -v nan输出传递给程序grep，该命令将删除对应于无法计算等位基因频率的位点所在的行。之后，将输出传递给程序gzip进行压缩| gzip，然后导入到一个名为的新文件中NC_031969.f5.sub1.pops1.tsv.gz（这里，“pops1”用于区分以后的分析与其他四种物种的集合，然后你可以称之为“pops2”，“pops3”等）。要查看脚本的其他选项freq.py，您可以输入python genomics_general/freq.py -h

如果你要计算D统计量，Fd统计量，还有其P值，这都基于一个假设，D统计量是0，现在会用到一个脚本calculate_abba_baba.r。其中在计算P value时，这个脚本需要你输入每一个染色体的长度，这里要继续新建一个文件chr_lengths.txt，并将对应染色体的长度输入进去。

  NC_031965 38372991
  NC_031966 35256741
  NC_031967 14041792
  NC_031968 54508961
  NC_031969 38038224
  NC_031970 34628617
  NC_031971 44571662
  NC_031972 62059223
  NC_031973 30802437
  NC_031974 27519051
  NC_031975 32426571
  NC_031976 36466354
  NC_031977 41232431
  NC_031978 32337344
  NC_031979 39264731
  NC_031980 36154882
  NC_031981 40919683
  NC_031982 37007722
  NC_031983 31245232
  NC_031984 36767035
  NC_031985 37011614
  NC_031986 44097196
  NC_031987 43860769
  UNPLACED  141274046

接着使用R脚本calculate_abba_baba.r来计算对应的D统计量，Fd统计量还有P-value。

Rscript calculate_abba_baba.r NC_031969.f5.sub1.pops1.tsv.gz pops1.abba_baba.txt ${spc1} ${spc2} ${spc3} ${spc4} chr_lengths.txt

注意：该脚本calculate_abba_baba.r假定传递给它的第一个命令行参数是具有等位基因频率的文件的名称，第二个是输出文件，第三个到第六个是“spc1”，“spc2”，“spc3”的ID ，和“spc4”按此顺序排列，第七个是有关染色体长度信息的文件。

简单看看生成的结果文件：

###less pops1.abba_baba.txt

Species 1: altfas
Species 2: neocan
Species 3: telvit
Species O: astbur

Number of sites: 80964
Number of BBAA sites: 13331.140625
Number of ABBA sites: 7419.078125
Number of BABA sites: 913.078125
D statistic: 0.780830292278784
fd statistic: 0.473484314376555
p-value: 0

你将看到分析中使用的四种物种ID的列表，然后是分析中使用的站点总数以及“BBAA”，“ABBA”和“BABA”位点的数量。将Neolamprologus cancellatus（“neocan”）称为“spc2”，将Telmatochromis vittatus（“telvit”）称为“spc3”，“ABBA”位点的数量约为7,500，但其他数字将取决于您使用的物种"SPC1"的种类。如果“BBAA”位点的数量不大于“ABBA”站点的数量，则可能需要calculate_abba_baba.r在切换“spc1”和“spc3”后使用脚本重复分析以获得正确的D -statistic。

接着在“spc1”的位置用不同的物种重复最后的步骤，找到具有最强基因渗入信号的物种，用D-统计量法测量。因此，再次使用以下命令，这次使用不同的物种ID而不是“XXXXXX”：然后，重复后续步骤，但在每个命令中将生成文件的命名“pops1”替换为“pops2”，“pops3”等，以避免覆盖先前的结果。

问题1：能否确定具有最大值的D统计量的物种（请注意，D统计量只有在“BBAA”位点的数量大于其他两个数字时才有效）？

答案：
在“spc1”的位置，Altolamprologus fasciatus（“altfas”）应观察到最高的D-统计量。实际上，这种物种与Neolamprologus cancellatus（“neocan”）相比具有比任何其他物种更多的遗传变异，Neolamprologus cancellatus（“neocan”）和Altolamprologus fasciatus（“altfas”）共享衍生物的“BBAA”位点的数量等位基因超过13,300，相比之下，Neolamprologus cancellatus（“neocan”）和Telmatochromis vittatus（“telvit”）分享了衍生等位基因的大约7,500个位点。该0.78的D统计量非常高，相应地，没有发生基因渐渗的假设的p值恰好为0，这意味着可以安全地拒绝该假设(D统计量为0）。值得注意的是，估计混合的基因组比例的Fd统计值约为0.47，接近50％。

在继续下一步之前，请确保现在指定了四个物种，以便将导致最高D统计量的物种放置在“spc1”的位置。因此，在相应更换“altfas”后再次执行以下操作：

spc1=altfas; spc2=neocan; spc3=telvit; spc4=astbur

然后，计算穿过染色体的滑动窗口中的Fd统计量，看看是否可以识别已经接收特别高或低渐渗程度的区域。这些区域可分别表示渗入等位基因的适应性或针对它们的选择。要在滑动窗口中计算f d，使用ABBABABAwindows.py Python脚本：

python genomics_general/ABBABABAwindows.py -g NC_031969.f5.sub1.geno.gz -o pops1.abba_baba.windows.txt -P1 ${spc1} -P2 ${spc2} -P3 ${spc3} -O ${spc4} --popsFile samples.txt -f phased -w 500000 --windType coordinate --minSites 500

在上面的命令中，再次指定基因型格式的输入文件-g NC_031969.f5.sub1.geno.gz，输出文件用-o和命名pops1.abba_baba.windows.txt，并且指定“spc1”，“spc2”，“spc3”和“spc4”的种类ID -P1，-P2，-P3，和-P4，分别。此外，为物种分配样品的文件提供了选项--popsFile samples.txt，选项-w 500000指定半个Mb作为窗口的大小，--windType coordinate指定窗口大小应以坐标而不是SNP的数量来测量，并--minSites 500指定SNP 的最小数量在窗口内。

然后将Fd统计量的统计结果可视化，可以使用该R脚本plot_abbababa_windows.r：

Rscript plot_abbababa_windows.r pops1.abba_baba.windows.txt chr_lengths.txt pops1.abba_baba.windows.pdf

问题：根据该图我们可以看到什么Fd统计量在染色体中的规律？

它在几乎所有的窗口保持0.4和0.6之间。没有进一步的证据，如果将选择作为任何高峰的驱动因素可能是对该模式过度的解释。你可以使用不同的窗口大小（例如100,000bp或1Mb）重复最后两个步骤，看看这是否会改变染色体上Fd统计量的变异模式。