远缘袖蝶属杂交和基因渗入的研究方法
The Research Methods to Detect Hybridization and Introgression in the Heliconius Butterflies   

研究背景

基因渗入 (introgression) 是种间杂交产生的基因交换，是向受体物种引入供体基因组片段的一种有效途径 (Mallet, 2005; Harrison and Larson, 2014)。相对于不适应的 (maladaptive) 或中性 (neutral) 的遗传片段，对生物生存和繁殖有益的基因组片段有更大的几率被固定下来，形成适应性的基因渗入 (adaptive introgression)，是基因组中适应性遗传变异的重要来源 (Martin and Jiggins, 2017)。其研究难点在于区别基因渗入与不完全谱系分选。随着测序技术的迅速发展，已经在多个研究体系中展开基因组层面的基因渗入研究，例如在现代人类与尼安德特人之间 (Fu et al., 2015)、三刺鱼近缘种间 (Jones et al., 2012)、疟蚊近缘种间 (Fontaine et al., 2015) 等。本研究方法利用进化基因组学理论和最新技术，整合系统发生和群体遗传方法，在非模式动物南美袖蝶远缘种间建立了一套分析流程，能够鉴定得到局部的基因渗入片段并判断  其方向，且能够有效区分基因渗入和不完全谱系分选，以此深入研究远缘物种间的杂交现象。
        本文提供了详细的基于全基因组数据集，检测远缘物种间基因渗入的实验流程和数据分析方法，该方法的有效性已在鳞翅目类昆虫中得以验证，详细结果请参阅Zhang et al., 2016.

仪器设备

注：如不需要样品制备，直接对已有数据进行分析，可跳过设备1-5，只配备和使用设备6。

1. PCR仪
2. 电泳仪
3. 超声波DNA破碎仪 (Covaris S220)
4. 生物分析仪 (Agilent 2100)
5. Illumina HiSeq 2500型测序仪
6. 服务器 (catalog number: ThinkSystemSR650；系统：Gentoo；CPU：Xeon 6230R 2.1GHz 26核；内存：512 Gb)
   

软件

注：本方法所使用的软件需要在服务器预先安装或下载直接使用，具体的下载和安装流程请参考软件对应的链接。

1. bcl2fastq2 Conversion Software (v2.20) https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html
2. Trimmomatic (v0.38) http://www.usadellab.org/cms/?page=trimmomatic
   
3. Bowtie2 (v2.3.4.3) http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
4. SAMtools (v1.9) http://www.htslib.org/
5. Picard (v1.84) http://broadinstitute.github.io/picard/
6. GATK (v3.3) https://gatk.broadinstitute.org
7. VCFtools https://vcftools.github.io/index.html
8. Phylogenomic https://github.com/npchar/Phylogenomic
9. RAxML https://cme.h-its.org/exelixis/web/software/raxml/index.html
10. iTOL https://itol.embl.de/
11. PhyTime https://www.stat.auckland.ac.nz/~guindon/phytime/
12. PhyML (v3.1) http://www.atgc-montpellier.fr/phyml/
13. TreeAnnotator https://beast2.blogs.auckland.ac.nz/treeannotator/
14. Genomics_general https://github.com/simonhmartin/genomics_general
15. R (v3.6.2) https://www.r-project.org
16. Perl (v5.32) https://www.perl.org/get.html
17. Python (v3.6) https://www.python.org/downloads/
    

实验步骤

一、样品制备和测序

注：如不需要样品制备，直接对已有数据或本方法提供的测试数据进行分析，可跳过步骤一。

1. 基因组DNA提取
   以袖蝶属 (Heliconius) 蝴蝶样品为例，收集成虫个体的胸部组织，采用DNeasy Blood & Tissue Kit (catalog number: 69504，Qiagen)，参照产品说明书的步骤 (https://www.qiagen.com/us/shop/pcr/dneasy-blood-and-tissue-kit/#resources)，提取基因组DNA，采用常规琼脂糖凝胶电泳检验提取的DNA产物，保存于-20 °C。
   
2. 测序文库构建
   采用已制备的基因组DNA，使用TruSeq DNA Library Prep Kits (catalog number: FC-121-2001，Illumina)，参照产品说明书的步骤 (https://support.illumina.com/downloads/truseq_dna_sample_preparation_guide_15026486.html)，构建双末端Illumina测序文库。
   
3. 上机测序
   所建文库使用Illumina HiSeq2500测序平台进行测序。Illumina测序仪下机所得的原始数据通常为bcl格式文件 (basecall file)，但下游分析一般需要fastq格式文件。因此，在进行下游分析前，使用bcl2fastq2 Conversion Software (v2.20) 将bcl文件根据之前添加的索引 (index) 分出，并转为fastq格式文件。命令行格式为：
   /usr/local/bin/bcl2fastq –runfolder-dir <运行文件夹及路径> --output-dir <输出文件夹及路径>
   

二、数据回贴和获得基因型信息

1. 数据获得
   本方法可以使用由步骤一获得的测序数据，或者使用发表于Genome Biology (2016年) 的文章 (Genome-wide introgression among distantly related Heliconius butterfly species) ( Zhang et al., 2016) 的部分数据作为测试数据。该数据可直接在NCBI的SRA数据库下载获得 (数据登录号：PRJNA308754)。
   
2. 数据质量控制
   使用Trimmomatic (v0.38) 去除原始序列数据中的低质量碱基和接头序列。去除序列首尾碱基质量低于3的碱基，滑窗去除平均碱基质量值低于10的90 bp窗口，去除长度低于36 bp的reads，得到高质量的基因组测序数据。命令格式为：
       java -jar trimmomatic-0.38.jar PE <输入fastq压缩文件1及路径> <输入fastq压缩文件2及路径> <输出fastq压缩文件1及路径> <输出未配对fastq压缩文件1及路径> <输出fastq压缩文件2及路径> <输出未配对fastq压缩文件2及路径> ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads SLIDINGWINDOW:90:10 LEADING:3 TRAILING:3 MINLEN:36
   
3. 数据回贴至参考基因组
   使用Bowtie2 (v2.3.4.3) 及参数-very-sensitive-local将序列数据比对至H. melpomene v2参考基因组 (下载地址：http://ensembl.lepbase.org/Heliconius_melpomene_melpomene_hmel2/Info/Index)。首先需生成参考基因组的索引文件，命令格式为：
       bowtie2-build <参考基因组fasta文件名及路径> <索引文件基名及路径>
           之后对每个个体基因组重测序数据，进行回贴并生成sam文件，命令格式为：
       bowtie2 --very-sensitive-local -p 8 -x <索引文件基名及路径> -1 <fastq文件1及路径> -2 <fastq文件2及路径> -S <输出sam文件名及路径>
   
4. 数据格式转换
   使用SAMtools (v1.9) 将sam文件转为bam文件。首先需生成参考基因组的索引文件，命令格式为：
           samtools faidx <参考基因组fasta文件名及路径>
           之后对于需转为bam的sam文件进行操作，命令格式为：
         samtools view -bST <参考基因组fasta文件名及路径> -o <输出的bam文件名及路径> <输入的sam文件名及路径>
   
5. 添加数据信息和去除PCR重复
   使用Picard (v1.84) 的AddOrReplaceReadGroups.jar命令对bam文件的reads添加信息，并使用MarkDuplicates.jar去除数据中的PCR重复。每一步均会生成一个新的bam文件作为后续的输入文件，并且每一步生成新的bam文件之后需采用BuildBamIndex.jar生成bam文件的索引，即bai文件。
           使用AddOrReplaceReadGroups.jar的格式为：
           java -jar /usr/share/picard-tools-1.84/AddOrReplaceReadGroups.jar INPUT=<输入bam文件名及路径> OUTPUT=<输出bam文件名及路径> SORT_ORDER=coordinate RGID=<用户添加ID信息> RGLB=<用户添加建库信息> RGPL=<用户添加平台信息> RGSM=<用户添加样品信息> RGPU=none VALIDATION_STRINGENCY=LENIENT
           使用MarkDuplicates.jar的格式为：
        java -jar /usr/share/picard-tools-1.84/MarkDuplicates.jar INPUT=<输入bam文件名及路径> OUTPUT=<输出bam文件名及路径> METRICS_FILE=<输出metrics文件名及路径> REMOVE_DUPLICATES=true ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT
           使用BuildBamIndex.jar的格式为：
          java -jar /usr/share/picard-tools-1.84/BuildBamIndex.jar INPUT=<输入bam文件名及路径> VALIDATION_STRINGENCY=LENIENT
   
6. 使用GATK流程获得基因型数据集
   在使用GATK流程处理个体基因组重测序bam文件之前，先要使用Picard (v1.84) 的CreateSequenceDictionary.jar命令对参考基因组fasta文件生成dict文件。
           使用Picard (v1.84) 的CreateSequenceDictionary.jar命令的格式为：
           java -jar /usr/share/picard-tools-1.84/CreateSequenceDictionary.jar R=<参考基因组fasta文件名及路径> O=<参考基因组dict文件名及路径>
           之后再使用GATK (v3.3) 中的RealignerTargetCreator和IndelRealigner命令重新比对每一个bam文件序列中的插入缺失区域。
           使用RealignerTargetCreator的格式为：
           java -jar /usr/share/GenomeAnalysisTK-3.3-0/GenomeAnalysisTK.jar –T RealignerTargetCreator -nt 10 -I <输入bam文件名及路径> -R <参考基因组fasta文件名及路径> -o <输出intervals文件名及路径>
           使用IndelRealigner的格式为：
           java -jar /usr/share/GenomeAnalysisTK-3.3-0/GenomeAnalysisTK.jar -T IndelRealigner -I <输入bam文件名及路径> -R <参考基因组fasta文件名及路径> -targetIntervals <输入intervals文件名及路径> -o <输出bam文件名及路径> -maxReads 100000
          最后并使用UnifiedGenotyper命令对预处理过的所有个体基因组重测序数据的bam文件进行基因分型并生成vcf格式的数据集，本方法使用了以下参数：heterozygosity 0.01, stand_call_conf 50, stand_emit_conf 10, dcov 250。以包括四个个体bam文件为例，具体的使用格式为：
          java -jar /usr/share/GenomeAnalysisTK-3.3-0/GenomeAnalysisTK.jar -T UnifiedGenotyper -nt 10 -R <参考基因组fasta文件名及路径> -I <个体1bam文件名及路径> -I <个体2bam文件名及路径> -I <个体3bam文件名及路径> -I <个体4bam文件名及路径> --heterozygosity 0.01 -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 250 -o <输出vcf文件名及路径>
   
7. 筛选高质量的单核苷酸多态性 (SNP) 位点
   使用VCFtools筛选高质量的SNP位点 (Qual > 30) 用于后续分析。具体使用格式为：
           vcftools --vcf <输入vcf文件名及路径> --minQ 30 --out <输出vcf文件名及路径> --recode
   

三、全基因组系统发生分析

基于全基因组SNP数据的系统发生分析适用于分析分歧时间相对较短的物种，可以获得基因组水平的界定物种关系的系统发生树，为研究物种间杂交和基因渗入时选取合适的系统发生拓扑结构提供依据。

1. 采用vcftools保留并对齐在所有个体中都得到基因型信息的高质量SNP位点，具体的使用格式为：
            vcftools --vcf <输入vcf文件名及路径> --minQ 30 --max-missing-count 0 --out <输出vcf文件名及路径> --recode
   
2. 采用GATK (v3.3) 中的FastaAlternateReferenceMaker 提取SNP信息并转换为fasta文件，并进一步采用Phylogenomic软件包中的fasta2relaxedPhylip.pl脚本转换为phylip格式文件。具体使用格式为：
          java -Xmx2g -jar GenomeAnalysisTK.jar -R<参考基因组fasta文件名及路径> -T FastaAlternateReferenceMaker -o <输出fasta文件名及路径> -L <输入intervals文件名及路径> -V <输入vcf文件名及路径>
       perl ./fasta2relaxedPhylip.pl -f <输入fasta文件名> -o <输出phylip文件名>
   
3. 使用RAxML软件构建全基因组最大似然树。使用的核苷酸替换模型为GTRGAMMA，进行100次快速自展 (bootstrap) 重复。具体使用格式为：
            raxmlHPC -f a -m GTRGAMMA -p 12345 -x 12345 -# 100 -T 20 -s <输入phylip文件名> -n test
4. 使用在线工具iTOL (https://itol.embl.de/) 对系统发生树进行可视化。
5. 使用PhyTime软件对基因组水平的拓扑结构进行时间校正，以H. cydno和H. melpomene 以及和H. pachinus之间的平均分歧时间 (1.4 Mya及0.43 Mya) 作为校准点，估算其他类群之间的分歧时间，并利用TreeAnnotator对结果进行统计。具体使用格式为：
       phytime -i <输入phylip文件名> -d nt -q -m GTR --calibration <输入校正节点文件名> -t e -v e -u <输入系统发生树文件名> --r_seed 35
   

四、检验种间基因流

为了更好的在全基因组水平鉴定潜在的基因渗入位点，建议使用帕特森D统计量 (Patterson’s D-statistic) 及f统计量 (modified f-statistic，fd) 对染色体和局部基因组区域进行统计检验，并进一步整合序列分歧水平等信息，以获得可靠的渗入位点。其中：

1. D统计量通过比较满足ABBA和BABA模式的衍生等位基因分布情况，检测基因流存在与否。以含有四个分类群的拓扑结构 (((P1, P2), P3), O) 为例，ABBA和BABA模式分别代表了P3与P2和P3与P1之间共享衍生等位基因。在零假设下，两种模式的数量应该相同。如两者之间存在显著差异，即表明对应类群之间发生了基因渗入。对于含有多个个体的群体数据，使用等位基因频率代替替换数量进行计算，计算公式如下：
   
            其中，P₁、P₂、P₃和P₄为分析中使用的四个类群，P(^)_ij为SNP i在种群j中观察到的等位基因频率。对于基因组水平的基因流估计，使用长度为50 kb的滑动窗口计算D统计量，并利用R软件包bootstrap (v. 2012-04) 以刀切法 (jackknife) 计算染色体水平的均值和标准误。
            其中，D统计量可以采用Genomics_general软件包进行计算，首先采用parseVCF.py转换文件格式，进而采用ABBABABAwindows.py计算，具体使用格式为：
         python parseVCF.py -i <输入压缩vcf文件名> --skipIndels --minQual 30 --gtf flag=DP min=5 | bgzip > <输出压缩geno文件名>
       python ABBABABAwindows.py -g <输出压缩geno文件名及路径> -f phased -o <输出csv文件名> -w 50000 -m 100 -s 50000 -P1 A -P2 B -P3 C -O D -T 10 --minData 0.5 --popsFile <输入种群信息文件名> --writeFailedWindows --polarize &
           为了鉴定不同尺度的局部渗入位点，也使用5 kb、10 kb和50 kb等多种不同的窗口大小计算D统计量。当窗口长度较小，处于连锁不平衡的范围内时，相邻窗口之间可能存在相关性。在这种情况下，利用R软件包wwntests (v.1.0.1) 中的block_bootstrap.R以滑动分块自助法 (moving-block bootstrap) 计算窗口内部的均值和标准误，其中分块大小设置为n^1/3，n为数据集大小。若均值在双尾z检验中显著偏离0，则认为该位点存在基因流。由于对基因组内的大量窗口进行统计检验，使用Benjamini-Hochberg方法对概率值进行多重比较校正，以错误发现率 (false discovery rate, FDR) 0.01作为存在显著差异的阈值。
   
2. fd统计量较D统计量有更高的分辨率，并且能够更好的估计渗入的比例，因此选择使用fd统计量进行进一步的筛选。fd统计量计算公式如下：
   
           其中，P₁、P₂、P₃和O为分析中使用的四个类群，P_D代表了渗入的供体种群，为P₂和P₃中具有更高衍生等位基因频率的群体。f_d计算方法同D统计量，可以采用Genomics_general软件包进行计算。
           同样以5 kb、10 kb和50 kb的窗口大小进行统计检验。检验方法同D统计量。
   
   
   图 1. 基因组范围内帕特森D统计量和f_d统计量示意图
   
   
3. 为了排除由于不完全谱系分选造成的假阳性信号，计算类群之间的序列绝对分歧 (sequence divergence, d_xy)，作为ABBA-BABA分析的补充。对于发生基因渗入的区域，序列分歧时间小于两个类群之间的分歧时间，因此通常具有较低的dxy水平。使用50 kb的滑动窗口计算染色体水平的dxy均值，对于潜在的渗入位点，使用100 bp窗口计算dxy，并利用R软件包的曼-惠特尼秩和检验 (Mann-Whitney U-test) 与所在染色体进行比较，去除分歧水平高于染色体均值的位点。对于群体数据，dxy的计算公式如下：
   
           其中，p_x和p_y分别为x和y类群中参考等位基因的频率，各等位基因频率可使用vcftools计算获得，具体使用格式为：
           vcftools --vcf <输入vcf文件名及路径> --keep <待计算的样品名文件> --freq --out <输出文件名及路径>。
   
4. 进一步检查潜在渗入位点的测序深度，以排除序列错误回贴造成的假阳性信号。由于与参考序列差异较大的重测序数据更加难以有效回贴，因此在基因型鉴定过程中，可能会基于与参考序列相似的数据得到对等位基因频率的错误估计。这种偏差也可能导致显著的渗入信号，特别是从参考基因组H. melpomene向其他物种的基因流。这些受影响的区域也会表现出序列覆盖深度的降低，因此可以通过位点的深度进行筛选。另一方面，具有异常高覆盖的区域也可能代表了序列的比对错误，并导致对等位基因频率的错误估计。因此，过滤去除测序深度低于5和高于40的位点。可以使用vcftools的参数--site-depth以及--site-mean-depth对测序深度进行评估；使用vcftools的参数--minDP，--maxDP以及--min-meanDP，--max-meanDP按照测序深度对位点进行过滤。
5. 对候选的渗入位点，使用PhyML (v3.1) 软件构建最大似然树，并根据系统发育关系推断基因流的方向。使用的核苷酸替换模型为GTR，进行100次自展分析。具体使用格式为：
   /usr/share/PhyML-3.1/PhyML-3.1_linux64 -i <输入phylip文件名及路径> --sequential -n 1 -b 100 -m GTR -t e -v e -o tlr --print_site_lnl --run_id 1
   

五、渗入位点功能注释

1. 基于H. melpomene基因组注释，获取渗入位点内部的蛋白编码基因。可在Lepbase (http://blast.lepbase.org/) 数据库检索基因已知的结构和功能。
2. 通过blastx将基因编码序列比对到NCBI nr (https://blast.ncbi.nlm.nih.gov/Blast.cgi)、Uni-Prot (https://www.uniprot.org/) 等蛋白数据库，基于蛋白序列的同源性，推断渗入基因的潜在功能。
3. 基于GO (http://geneontology.org)、KEGG (http://www.genome.jp/kegg) 和Pfam (v33.1,http://pfam.xfam.org) 等数据库对渗入基因进行富集分析。

致谢

感谢张宇博对本文撰写提供的帮助。感谢审稿人在本文修改中提出的建议和意见。本研究得到北京市自然科学基金杰出青年项目（编号JQ19021）、国家自然科学基金面上项目（编号31871271）、启东创新基金支持。

参考文献

1. Mallet, J. (2005). Hybridization as an invasion of the genome. Trends Ecol Evol 20: 229-237.
2. Harrison, R. G. and Larson, E. L. (2014). Hybridization, introgression, and the nature of species boundaries. J Hered 105: 795-809.
3. Martin, S. H. and Jiggins, C. D. (2017). Interpreting the genomic landscape of introgression. Curr Opin Genet Dev 47: 69-74.
4. Fu, Q., Hajdinjak, M., Moldovan, O. T., Constantin, S., Mallick, S., Skoglund, P., Patterson, N., Rohland, N., Lazaridis, I., Nickel, B., Viola, B., Prüfer, K., Meyer, M., Kelso, J., Reich, D. and Pääbo, S. (2015). An early modern human from Romania with a recent Neanderthal ancestor. Nature 524: 216-219.
5. Jones, F. C., Grabherr, M. G., Chan, Y. F., Russell, P., Mauceli, E., Johnson, J., Swofford, R., Pirun, M., Zody, M. C., White, S., Birney, E., Searle, S., Schmutz, J., Grimwood, J., Dickson, M. C., Myers, R. M., Miller, C. T., Summers, B. R., Knecht, A. K., Brady, S. D., Zhang, H., Pollen, A. A., Howes, T., Amemiya, C., Broad Institute Genome Sequencing Platform & Whole Genome Assembly Team, Baldwin, J., Bloom, T., Jaffe, D. B., Nicol, R., Wilkinson, J., Lander, E. S., Palma, F. D., Lindblad-Toh, K. and Kingsley, D. M. (2012). The genomic basis of adaptive evolution in threespine sticklebacks. Nature 484: 55-61.
6. Fontaine, M. C., Pease, J. B., Steele, A., Waterhouse, R. M., Neafsey, D. E., Sharakhov, I. V., Jiang, X., Hall, A. B., Catteruccia, F., Kakani, E., Mitchell, S. N., Wu, Y. C., Smith, H. A., Love, R. R., Lawniczak, M. K., Slotman, M. A., Emrich, S. J., Hahn, M. W. and Besansky, N. J. (2015). Extensive introgression in a malaria vector species complex revealed by phylogenomics. Science 347: 1258524.
7. Zhang, W., Dasmahapatra, K. K., Mallet, J., Moreira, G. R. P. and Kronforst, M. R. (2016). Genome-wide introgression among distantly related Heliconius butterfly species. Genome Biol 17: 25.