重测序（RADseq）做群体遗传分析套路

Posted 2023-03-30

参考技术A

构建的群体，或自然群体，如各地方品种。

提取DNA后，构建文库，简要步骤如下：
① 限制性内切酶TaqI酶切；
② 连接P1接头；
③ DNA随机打断片断化；
④ 目的片段回收与末端修复；
⑤ 连接P2接头；
⑥ RAD片段富集；
⑦ 上机测序。
参考： Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers

根据识别标签序列得到每个个体的测序reads，使用trimmomatic进行过滤（其他质控软件，如fastqc，multiQC等）
设置过滤参数为：SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50。 过滤标准：两端质量低于5的碱基进行切除，并以5bp为窗口进行滑动过滤，对平均质量低于20的窗口进行切除。

BWA （其他比对软件如bowtie2/soap2/MAQ等）将过滤后的个体clean reads比对到参考基因组序列上。样本比对率反映的是样本测序数据与参考基因组的相似性，覆盖深度和覆盖度能够直接反映测序数据的均一性与参考序列的同源性。

使用GATK（或samtools+bcftools)Haplotype Caller模块进行变异检测，获得群体变异集文件（VCF 格式）。对变异进行过滤：过滤参数为缺失率小于或等于0.2、杂合率小于或等于0.2、最小等位基因频率（MAF） 大于或等于0.05，最终得到高质量的基因型数据。

在获得高质量的标记数据以后，利用vcftools将vcf文件处理得到plink.ped和plink.map文件（整理为plink软件所需格式）。

使用plink 软件随机选择连锁不平衡（LD）小于0.1，且相邻间隔在300kb以上的SNP位点，最后得到一个包含3420 个SNP位点的标记集，一般是生成.bed文件。

1.祖先成分堆叠图
使用ADMIXTURE对此 SNP位点集（bed文件）进行群体结构分析（Structure），利用交叉验证过程确定确定合适的祖先数或亚群（K值）。若不知道理想的K值，可用ADMIXTURE计算，一般当cross-validation error值最低时所对应的K值为最合适的K值。

考虑到样本所归属的分类单元，即看看哪几个物种聚在一起，对合适的K值利用Structure软件（速度慢，其他软件如frappe，ADMIXTURE也可做群体结构图，并且很快）聚类图，一些R包如hapmap也是可以做群体结构图的。

2.PCA
利用GCTA对SNP数据集进行样本的PCA分析（其他软件如EIGENSOFT中的smartpca）。GCTA可以直接读取.bed , .bim , .fam文件，利用–make-grm 生成个体对之间的遗传关系矩阵，并将GRM的下三角元素保存为二进制文件.grm.id , .grm.bin ， .grm.N.bin。使用 –pca 设置要生成主成分的数目，一般来说就可以刻画出群体结构。这一步会生成 .eigenval 和 .eigenvec 两个文件。.eigenval文件为各主成分可解释遗传信息的比例，.eigenvec文件为每个样本在top4主成分上的分解值。

3.系统发育树
构树的方法有非加权分组平均法（UPGMA，已经很少用）、最小进化法（ME）、邻接法（NJ）、最大简约法（MP）、最大似然法（ML）等。

构树软件如FastTree/MEGA/cluster X/phylip，美化可以用FigTree/ggtree/treeview/GraPhIAn。

NJ法是基于最小进化原理经常被使用的一种算法，它不检验所有可能的拓扑结构，能同时给出拓扑结构和分支长度。

GWAS的群体遗传分析也是包含这三个图，RADseq毕竟是简化基因组，得到的SNP有限，做这种群体分析效果肯定没有GWAS好。