linux怎么估算基因组大小

ALLPATHS-LG的使用

一、ALLPATH简介

ALLPATHS-LG是一个基因组组装软件，适合于组装short reads数据，由Computational Research and Development group at the Broad Institute开发。ALLPATHS-LG是现在行业内公认进行基因组De novo组装效果最好的软件。

二. 基础注意事项

1. 不能只使用一个library数据进行组装； 2. 必须有一个"overlapping"的片段文库的paired-reads数据。比如，reads长度~ 100bp，插入片段库长度~180bp3. 必须有jumping library数据； 4. 基因组组装需要100x或以上基因组覆盖度的碱基，这个覆盖度是指raw reads数据(在 error correction和filtering之前)的覆盖度； 5. 可以使用PacBio数据； 6. 不能使用454数据和Torrent数据。主要是这两者测序太贵，如果什么时候价格降低，有 需求的话，会写出相应的代码来满足要求； 7. 官方提供了测试用数据； 8. 不支持在整个计算机集群上进行运算； 9. 需要消耗的内存峰值大约是1.7bytes每个碱基，即输入10G的碱基数据量，大约需要17 G内存； 10. 对于试探性的参数，比如K，原则上可以调整。但是我们不会自行调整，并也不推荐。AL LPATHS-LG不像其它De novo一样，Kmer大小的参数K和read大小之间没有直接的联系， ALLPATHS-LG会在运行过程中运用一系列的K值。

三. ALLPATHS-LG使用方法

1. 基础的使用方法和命令

使用RunAllPathsLG这个命令来运行。虽然有很多参数，但是在没有指导的情况下不要随意使用，使用默认设置即可。其使用方法为：

$ RunAllPathsLG arg1=value1 arg2=value2 ...

参数主要是设置程序辨别的一些目录，在程序的运行过程，会输入相应目录中的数据，将结果输入到指定的目录。一个简单的命令使用例子：

#!/bin/sh # ALLPATHS-LG needs 100 MB of stack space. In 'csh' run 'limit stacksize 100000'. ulimit -s 100000 # ALLPATHS-LG命令的写法与一般的linux参数写法不是很一样。采用 ‘参数=值’ 的方法，并使之成每行一个参数，使用'\'来连接各个参数，这样看起来直观易懂。初始接触的人可能会不适应。 RunAllPathsLG \ PRE=$PWD\ REFERENCE_NAME=species.genome\ DATA_SUBDIR=data\ RUN=run\ SUBDIR=test\ EVALUATION=STANDARD\ TARGETS=standard\ OVERWRITE=True\ MAXPAR=8 | tee -a assemble.out

2. 详细的参数说明

必须的参数 PRE (String) 程序运行的根目录，所有的其它目录全在该目录下REFERENCE_NAME (String) 参考基因组目录名称，位于PRE目录下。如果有一个参考基因组，可将参考基因组放到该 目录中；若没有，则创建该文件夹用于基因组组装DATA_SUBDIR (String) DATA子目录名称，位于REFERENCE_NAME目录下。程序从该目录中读取数据。 RUN (String) 运行目录名称，位于DATA_SUBDIR下。程序将生成的中间文件和结果文件存储于该目录 。比如组装结果是一个名为ASSEMBLES的目录，位于该目录下。 部分可选参数： SUBDIR (String) default: test 子目录名，在REF/DATA/RUN/ASSEMBLIES目录下创建的存放基因组组装结果的目录 名。 K (int) default: 96 核心Kmer大小，只有K=96能很好地运行。 EVALUATION (String: {NONE,BASIC,STANDARD,FULL,CHEAT})default:BASIC 给定一个参考基因组，pipeline能在基因组组装的不同阶段对组装过程和结果进行评估。 BASIC:基础评估，不需要参考基因组； STANDARD:使用参考基因组来运行评估模块； FULL:在某些组装模块下打开in-place评估，不会影响组装结果； CHEAT:稍微使用参考基因组指导组装，产生更详细的分析，能对组装结果产生小的(好方 向的)改变。REFERENCE_FASTA (String) default: REF/genome.fasta 评估中使用的参考基因组。 MAXPAR (int) default: 1 有些模块的运行是独立的，不相互依赖，能同时运行。该参数设定能同时运行的模块的最 大数目。由于pipeline中的绝大部分模块都能多线程运行，因此将该值设定大于1，效果不明 显。 THREADS (String) default: max 有些模块能多线程程运行，默认使用最大线程数运行。 OVERWRITE (Bool) default: False 是否覆盖存在的文件。可以设置该选项为True，在每次运行程序的时候设定RUN参数为 一个新的目录名，则比较好。 TARGETS (vec) default: standard pipeline会生成一系列的文件，不同的文件的生成需要call不同的模块。如果某文件 已经存在了并且是最新的，则跳过相应的模块的运行。本参数指定生成哪些拟定的目标文件(p seudo targets)。若目标文件没有相应的模块能生成，则会得到报错。 none:没有拟定的目标文件，仅仅生成指定的目标文件； standard:生成组装文件和选定的评估文件； full_eval:生成组装文件和额外的评估文件。TARGETS_REF (String) 在ref_dir目录中生成的目标文件。 多个目标文件的书写方法为： TARGETS_REF="{target1,target2,target3}" 。 TARGETS_DATA (String) 在data目录中生成的目标文件。 TARGETS_RUN (String) 在run目录中生成的目标文件。 TARGETS_SUBDIR (String) 在subdir中生成的目标文件。FORCE_TARGETS (Bool) default: False 生成目标文件，即使文件已经存在并且看起来是很新的。

3. 输入文件与目录的准备

两个文库：插入片段长度为180bp和3000bp，illumina测序文件结果为fastq格式。以此为例来准备ALLPATHS-LG运行所需的文件和目录。

(1) 准备 in_groups.csv 和 in_libs.csv 文件。

这两个文件内容由逗号隔开，in_groups.csv文件内容如下：

group_name, library_name, file_name firest, Illumina_180bp, seq/species_500bp_read?.fastq second, Illumina_3000bp, seq/species_3000bp_read?.fastq

in_groups.csv文件的解释：

group_name:数据独特的代号,每一份数据有一个代号； library_name:数据所属文库的名字，体现出该； filename:数据文件所存放位置。可以为相对位置，文件名可以包含'*'和'?'(但是扩展名 中不能有该符号，因为要根据扩展名识别文件类型)，从而代表paired数据。支持的文件类型有 '.bam','fasta','fa','fastq','fq','fastq.gz'和'fq.gz'。

in_libs.csv文件内容如下：

library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end Illumina_180bp, species, species.genome, fragment, 1, 180, 10, , , inward, 0, 0 Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0

in_libs.csv文件的解释：

library_name:和in_groups.csv中的相匹配； project_name:project的名字； organism_name:测序物种的名字； type:仅仅只是一个信息； paired:0:Unpaired reads1:paired readsfrag_size:小片段文库插入片段长度的均值； frag_stddev:小片段文库的插入片段长度估算的标准偏差； insert_size:大片段文库插入片段长度的均值； insert_stddev:大片段文库插入片段长度估算的标准偏差； read_orientation:reads的方向，小片段文库为inward，大片段文库为outward； genomic_start:reads从该位置开始，读入数据，如果不为0，之前的碱基都被剪掉； genomic_end:reads从该位置开始，停止读入数据，如果不为0，之后的碱基都被剪掉。

(2) 使用PrepareAllPathsInputs.pl来对数据进行转换

ALLPATHS-LG接受的输入数据要求如下：

1. ALLPATHS-LG的输入数据支持小片段文库(fragment library)、大片段文库(jum ping library)和超大片段文库(long jumping library)。并且前两种文库至少各有 一个才能进行基因组组装。超大片段文库是只插入片段>20kb的文库，其测序方向和小片段文 库一致，为inward。 2. ALLPATHS-LG的输入数据放置在//文件夹下，包含3种文件：碱基文件，质量文件和配 对信息文件 frag_reads_orig.fastb frag_reads_orig.qualb frag_reads_orig.pairs jump_reads_orig.fastb jump_reads_orig.qualb jump_reads_orig.pairs 以下是可选的超大插入片段文库对应的数据文件（非必须）： long_jump_reads_orig.fastb long_jump_reads_orig.qualb long_jump_reads_orig.pairs

使用PrepareAllPathsInputs.pl来将fastq等格式的测序结果转换成ALLPATHS-LG可接受的文件。以下是该程序的参数：

DATA_DIR 将转换后的数据文件放到此文件夹下。 PICARD_TOOLS_DIR 若输入数据为bam格式，则需要用到Picard软件，该参数Picard的路径 IN_GROUPS_CSV 输入的in_groups.csv文件名 IN_LIBS_CSV 输入的in_libs.csv文件名INCLUDE_NON_PF_READS default: 1 1:包含non-PF reads；0:仅仅只包含PF reads. PHRED_64 default: 0 0:碱基质量是ASCII的33到126，一般情况下Illumina数据的最低碱基质量是'B'1:碱基质量的ASCII码是从64到126，一般情况下Illumina数据的最低碱基质量是'#'。 PLOIDY 生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型，2表 示双倍体型。 HOSTS 列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3. host2,4.host3"表示使用本地机器的2个CPU线程，host2机器的3个CPU线程和host3机 器的4个CPU线程。 以下是不常用的参数，主要用来选择转换的数据量的大小。当测序数据量太多，而只想使用其 中一部分数据的时候，可以用到 FRAG_FRAC 使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值，则不能同时设定 FRAG_COVERAGE。 JUMP_FRAC 使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值，则不能同时设定 JUMP_COVERAGE。 LONG_JUMP_FRAC 使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值，则不能同时 设定LONG_JUMP_COVERAGE。 GENOME_SIZE 估计的基因组大小，用来计算对应覆盖度所对应的reads数 FRAG_COVERAGE 所期望的小片度库的覆盖度，比如 45. 要求GENOME_SIZE有设定 JUMP_COVERAGE 所期望的大片度库的覆盖度，比如 45. 要求GENOME_SIZE有设定 LONG_JUMP_COVERAGE 所期望的超大片度库的覆盖度，比如 1. 要求GENOME_SIZE有设定

分子层面对生物的研究，在个体水平上主要是看单个基因的变化以及全转录本的变化（RNA-seq）；在对个体的研究的基础上，开始了群体水平的研究。如果说常规的遗传学主要的研究对象是个体或者个体家系的话，那么群体遗传学则是主要研究由不同个体组成的群体的遗传规律。

在测序技术大力发展之前，对群体主要是依靠表型进行研究，如加拉巴哥群岛的13中鸟雀有着不同的喙，达尔文认为这是自然选择造成的后果 。达尔文的进化论对应的观点可以简单概括为“物竞天择，适者生存”，这也是最为大众所接受的一种进化学说。直到1968年，日本遗传学家提出了中性进化理论[2]，也叫中性演化理论。中性理论的提出很大程度上是基于分子生物化学的发展。可以这样理解中性理论：一群人抽奖，在没有内幕的情况下，每个人抽到一等奖的概率是相等的，这个可能性和参与抽奖的人的身高、年龄、爱好等因素都没有关系。中性理论常作为群体遗传研究中的假设理论（CK）来计算其他各种统计指标。

群体遗传学，研究的单位是群体，比如粳稻、籼稻、野生稻，就能够构成不同的群体；我们国内的各省份的水稻也可以作为一个个群体。 群体遗传学大概可以分为群体内的研究和群体间的研究。比如研究云南元阳的水稻的遗传多样性；如果研究是的云南元阳的水稻和东北的水稻，那就可以算成是群体间的研究。群体间和群体内的研究是相互的。

测序价格的急剧下降[3]使得大规模的群体测序得以实现。

常见的变异类型有SNP、IdDel、SV、CNV等。重测序中最关注的是SNP，其次是InDel。其他的几种结构变异的研究不是太多。

有参考基因组的物种的全基因组测序叫做重测序，没有参考基因组的物种的全基因组测序则需要从头组装。随着测序价格的降低，越来越多物种的参考基因组都已经测序组装完成。 plant genomes [4]网站实时显示全基因组测序已经完成的植物，其中2012年以后爆发式增长。在群体遗传学研究中更多的是有参考基因组的物种，尤其是模式物种，植物中常见的是拟南芥、水稻和玉米。

主要的分析流程见下图。现在的测序公司基本上都会帮客户完成整个的分析流程，因为主要耗费的资源是计算资源。我认为在整个分析的流程中最重要的是Linux目录的构建，混乱的目录会导致后续的分析频频出问题，重测序分析会生成很多的中间文件，良好的目录管理会使得项目分析流程井然有序。

该部分涉及到的软件的安装和基础的Linux基础知识就不详细说明了。

正选择似乎可以更好地用自然选择来解释。就是一个基因or位点能够使个体有着更强的生存力或者是育性，这样就会使得这个个体的后代更多，如此一来，这个基因or位点在群体中就越来越多。

正选择能够使有利的突变基因or位点在群体中得到传播，但是与此同时却降低了群体的多态性水平。也就是说原先该位点周围的核苷酸组成是多样性的，在经过正选择之后，这个位点周围核苷酸的多样性就渐渐的趋于同质化了。这就好比一块田，里面本来有水稻和稗草及其他杂草，由于稗草的适应性增强，稗草在逐渐增多，水稻慢慢变少，最后甚至是只剩下了稗草。

我们将这种选择之后多态性降低的情况叫做选择扫荡（Selective Sweep)。检测选择扫荡的软件有SweeD[7]。选择扫荡有可能是人工选择的结果，如2014年 Nature Genetics关于非洲栽培稻的文章就使用了SweeD来检测非洲栽培稻基因组上受人工选择的区域[8]。

负选择和正选择刚好是相反的。简单理解成群体中的某个个体出现了一个致命的突变，从而自己或者是后代从群体中被淘汰。这也导致群体中该位点的多态性的降低。就好比我有10株水稻，其中一株在成长过程中突然不见了，那么对我的这个小的水稻群体来说，这个消失的水稻的独有的位点在群体中就不见了，整体的多态性就降低了。

平衡选择指多个等位基因在一个群体的基因库中以高于遗传漂变预期的频率被保留，如杂合子优势。

平衡选择检测的算法有BetaScan2[10]，这是个Python脚本，输入文件只需要过滤好的SNP数据即可。

计算公式为：

其中 是有效群体大小， 是每个位点的突变速率。 但是群体大小往往是无法精确知道的，需要对其进行估计。

分离位点数 是 的估计值，表示相关基因在多序列比对中表现出多态性的位置。计算公式为：

其中 为分离位点数量，比如SNP数量。

为个体数量的倒数和：

指的是核苷酸多样性，值越大说明核苷酸多样性越高。通常用于衡量群体内的核苷酸多样性，也可以用来推演进化关系[11]。计算公式为：

可以理解成现在群体内两两求 ，再计算群体的均值。计算的软件最常见的是 vcftools ，也有对应的R包 PopGenome 。通常是选定有一定的基因组区域，设定好窗口大小，然后滑动窗口进行计算。

3KRGP文章就计算了水稻不同亚群间4号染色体部分区域上的 值[12]，能够看出控制水稻籽粒落粒性的基因 Sh4 位置多态性在所有的亚群中都降低了。说明这个基因在所有的亚群中都是受到选择的，这可能是人工选择的结果。

Tajima's D是日本学者Tajima Fumio 1989年提出的一种统计检验方法，用于检验DNA序列在演化过程中是否遵循中性演化模型[14]。计算公式为：

D值大小有如下三种生物学意义：

叫固定分化指数，用于估计亚群间平均多态性大小与整个种群平均多态性大小的差异，反映的是群体结构的变化。其简单估计的计算公式为：

的取值范围是[0,1]。当 时，表明亚群间有着明显的种群分化。

在中性进化条件下， 的大小主要取决于遗传漂变和迁移等因素的影响。假设种群中的某个等位基因因为对特定的生境的适应度较高而经历适应性选择，那该基因的频率在种群中会升高，种群的分化水平增大，使得种群有着较高的 值。

值可以和GWAS的结果一起进行分析， 超过一定阈值的区域往往和GWAS筛选到的位点是一致的，如2018年棉花重测序的文章[15]：

ROD可以基于野生群体和驯化群体间核苷酸多态性参数 的差异识别选择型号，也可以测量驯化群体和野生型群体相比损失的多态性。计算公式为：

和 一样，ROD也可以和GWAS结合起来：

群体结构分析可以简单理解成采样测序的这些个体可以分成几个小组，以及给每个个体之间的远近关系是怎么样的。群体结构分析三剑客， 分别是 进化树 、 PCA 和 群体结构图 。

进化树就是将个体按照远近关系分别连接起来的图。

常用的绘图软件是 Phylip 和 Snpphylo 。进化树修饰的软件有 MEGA ， ggtree 等，推荐网页版工具 iTOL ，无比强大。

外群定根法：当群体的个体的差异很小时，可以引入其他物种作为根。如在对三叶草建树时可以引入水稻的序列作为根进行建树。

PCA是很常见的降维方法，如微生物研究中常用来检验样品分群情况。PCA计算的软件很多，plink可以直接用vcf文件计算PCA，R语言也可以进行PCA计算。

PCA图在群体重测序中有如下几种作用：

进化树和PCA能够看出来群体是不是分层的，但是无法知道群体分成几个群合适，也无法看出群体间的基因交流，更无法看出个体的混血程度。这时候就需要群体分层图了。

可以将进化树和群体分层图结合进行展示，如下图：

先了解下概念，此处借鉴基迪奥生物网站的解释[22]。

要理解 LD 衰减图，我们就必须先理解连锁不平衡（Linkage disequilibrium，LD）的概念。连锁不平衡是由两个名词构成，连锁 + 不平衡。前者，很容易让我们产生概念混淆；后者，让这个概念变得愈加晦涩。因此从一个类似的概念入手，大家可能更容易理解 LD 的概念，那就是基因的共表达。

基因的共表达，通常指的是两个基因的表达量呈现相关性。比较常见的例子就是：转录组因子和靶基因间的关系。因为转录因子对它的靶基因有正调控作用，所以转录因子的表达量提高会导致靶基因的表达量也上调，两者往往存在正相关关系。这个正相关关系，可以使用相关系数来度量，这个数值在 - 1~1 之间。总而言之，相关性可以理解为两个元素共同变化，步调一致。

类似的，连锁不平衡（LD）就是度量两个分子标记的基因型变化是否步调一致，存在相关性的指标。如果两个 SNP 标记位置相邻，那么在群体中也会呈现基因型步调一致的情况。比如有两个基因座，分别对应 A/a 和 B/b 两种等位基因。如果两个基因座是相关的，我们将会看到某些基因型往往共同遗传，即某些单倍型的频率会高于期望值。

参照王荣焕等[23]的方法进行LD参数计算：

随着标记间的距离增加，平均的LD程度将降低，呈现出衰减状态，这种情况叫LD衰减。LD衰减分析的作用：

GWAS(genome-wide association study)，全基因组关联分析，常用在医学和农学领域。简单理解成将SNP等遗传标记和表型数据进行关联分析，检测和表型相关的位点，然后再倒回去找到对应的基因，研究其对表型的影响。这些被研究的表型在医学上常常是疾病的表型；在农学上常常是受关注的农艺性状，比如水稻的株高、产量、穗粒数等。GWAS思想首次提出是在心肌梗塞的治疗上[24]，首次应用是在2005年的文章上[25]。

目前使用最广泛的模型是混合线性模型[26]：

所有的参数软件（如Emmax）会自动完成计算。

GWAS结果文件通常只有两个图，一个是曼哈顿图，另外一个是Q-Q图。一般是先看Q-Q图，如果Q-Q正常，曼哈顿图的结果才有意义。

MSMC（multiple sequentially Markovian coalescent）[27]，底层算法很复杂，类似于PSMC。MSMC的主要功能是推断有效群体大小和群体分离历史。

这样看起来更直观：

LAMP(Local Ancestry in Admixed Populations，混杂群体的局部族源推断)，用于推断采用聚类的方法假设同时检测的位点间不存在重组情况，对每组相邻的 SNP 进行检测分析[28]，在运算速度和推断准确度上都有了质的飞跃。

用于推断群体分离和混合[29]。图是这样的：

测序方案关系到后续的分析，不同的样本量对应不同的测序方法和分析方法。

[1]. 自然选择(维基百科)

[2]. Kimura, Motoo. "Evolutionary rate at the molecular level." Nature . 217.5129 (1968): 624-626 .

[3]. 测序价格变化趋势

[4]. plant genomes

[5]. DePristo, Mark A., et al. "A framework for variation discovery and genotyping using next-generation DNA sequencing data." Nature Genetics . 43.5 (2011): 491.

[6]. Biswas, Shameek, and Joshua M. Akey. "Genomic insights into positive selection." ** TRENDS in Genetics. 22.8 (2006): 437-446.

[7]. Pavlidis, Pavlos, et al. "Sweed: likelihood-based detection of selective sweeps in thousands of genomes." Molecular biology and evolution 30.9 (2013): 2224-2234.

[8]. Wang, Muhua, et al. "The genome sequence of African rice (Oryza glaberrima) and evidence for independent domestication." Nature Genetics 46.9 (2014): 982.

[9]. Bamshad, Michael, and Stephen P. Wooding. "Signatures of natural selection in the human genome." Nature Reviews Genetics 4.2 (2003): 99.

[10]. Siewert, Katherine M., and Benjamin F. Voight. "BetaScan2: Standardized statistics to detect balancing selection utilizing substitution data." BioRxiv (2018): 497255.

[11]. Yu, N.Jensen-Seaman MIChemnick LRyder OLi WH (March 2004). Genetics . 166 (3): 1375–83.

[12]. Wang, Wensheng, et al. "Genomic variation in 3,010 diverse accessions of Asian cultivated rice." Nature 557.7703 (2018): 43.

[13]. Li, C., Zhou, A. &Sang, T. Rice domestication by reducing shattering. Science 311, 1936–1939 (2006).

[14]. Tajima, Fumio. "Statistical method for testing the neutral mutation hypothesis by DNA polymorphism." Genetics 123.3 (1989): 585-595.

[15]. Du, Xiongming, et al. "Resequencing of 243 diploid cotton accessions based on an updated A genome identifies the genetic basis of key agronomic traits." Nature Genetics 50.6 (2018): 796.

[16]. Lu, Kun, et al. "Whole-genome resequencing reveals Brassica napus origin and genetic loci involved in its improvement." Nature communications . 10.1 (2019): 1154.

[17]. Zhou, Z., Jiang, Y., Wang, Z. et al. Resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean. Nat Biotechnol 33, 408–414 (2015).

[18]. Liang, Z., Duan, S., Sheng, J. et al. Whole-genome resequencing of 472 Vitis accessions for grapevine diversity and demographic history analyses. Nat Commun 10, 1190 (2019).

[19]. Alexander, D.H., Lange, K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12, 246 (2011).

[20]. Francis, Roy M. "pophelper: an R package and web app to analyse and visualize population structure." Molecular ecology resources 17.1 (2017): 27-32.

[21]. http://www.royfrancis.com/pophelper/articles/index.html .

[22]. https://www.omicshare.com/forum/thread-878-1-1.html .

[23]. WANG Rong-Huan, WANG Tian-Yu, LI Yu. Linkage disequilibrium in plant genomes[J]. HEREDITAS , 2007, 29(11): 1317-1323.

[24]. Ozaki, K., Ohnishi, Y., Iida, A. et al. Functional SNPs in the lymphotoxin-α gene that are associated with susceptibility to myocardial infarction. Nat Genet 32, 650–654 (2002).

[25]. Klein, Robert J., et al. "Complement factor H polymorphism in age-related macular degeneration." Science 308.5720 (2005): 385-389.

[26]. Yu, Jianming, et al. "A unified mixed-model method for association mapping that accounts for multiple levels of relatedness." Nature genetics 38.2 (2006): 203.

[27]. Schiffels, Stephan, and Richard Durbin. "Inferring human population size and separation history from multiple genome sequences." Nature genetics 46.8 (2014): 919.

[28]. Sankararaman, Sriram, et al. "Estimating local ancestry in admixed populations." The American Journal of Human Genetics 82.2 (2008): 290-303.

[29]. Pickrell, Joseph K., and Jonathan K. Pritchard. "Inference of population splits and mixtures from genome-wide allele frequency data." PLoS genetics 8.11 (2012): e1002967.

[30]. Chen, Jia-hui, et al. "Genome-wide analysis of Cushion willow provides insights into alpine plant divergence in a biodiversity hotspot." Nature communications 10.1 (2019): 1-12.

[31]. 孙宽，侯一平。法医族源推断的分子生物学进展 [J]. 法医学杂志 ，2018,34 (03):286-293.

[32]. genek.tv

---