微卫星标记的开发和数据分析流程

仪器设备

1. Windows系统电脑

软件信息

1. Java: https://www.oracle.com/java/technologies/javase-downloads.html
   
2. Perl: https://www.activestate.com/products/perl/downloads/
   
3. R: https://mirrors.tuna.tsinghua.edu.cn/CRAN/
   
4. Rstudio: https://rstudio.com/products/rstudio/download/
   
5. GMATA v2.2: https://github.com/XuewenWangUGA/GMATA
   
6. GenAlEx 6.5: http://biology-assets.anu.edu.au/GenAlEx/Download.html
   
7. Microsatellite toolkit: https://www.researchgate.net/profile/Benjamin-Barth-2/post/Where-can-I-download-the-excel-micro-satellite-toolkit/attachment/59d635dac49f478072ea3998/AS%3A273668205154304%401442258992481/download/MStools.zip
   
8. Genepop on the web: https://genepop.curtin.edu.au
   
9. MicroChecker 2.2: https://micro-checker.software.informer.com/
   
10. GelQuest: https://www.sequentix.de/gelquest/
    
11. GeneMarker: https://www.softgenetics.com/GeneMarker.php
    

实验步骤

1. 微卫星标记的开发
   对于没有任何参考序列信息的物种，常采用重复序列特异探针 + 磁珠富集的方法，该方法实验过程较为繁琐。随着测序成本的降低，我们可使用高通量测序技术，测定和组装物种基因组的随机片段，来获取微卫星位点。
           更常见的情况是，我们可以从已知序列数据集 (例如参考基因组、简化基因组、转录组) 搜索微卫星位点。这里我们以NCBI下载的参考基因组为例，介绍微卫星标记的开发流程。对于其他序列数据集，方法是类似的。
   
   1.1

   从NCBI搜索并下载fasta格式的参考基因组序列。
   1.2

   使用GMATA v2.2软件 (Wang and Wang, 2016) 进行位点开发。按安装说明书GMATA_installation.pdf的提示，安装perl, R和JAVA语言环境。配置好后双击文件夹中的GMATAv2.2.jar运行。
   1.3

   参考说明书GMATA manual_V2.01_20151218.pdf的3.2章节进行操作，依次将每一步生成的文件导入下一步。
           设置重复单元nt值的一个通用原则是：重复单元越小，重复次数越多，多态性就越高。但是重复次数越多，最终位点的获取量就越少。对于简化基因组或基因组草图数据，因为序列gap较多，为了保证产出足够的微卫星位点，建议多次尝试不同的值。
           下面给出了SSR identification和Marker designing的参考值。如果DNA质量不佳，建议将-Smax值降低为220，因为低质量的微卫星分型结果经常会出现在220 bp以上的等位基因中。
           Min-length(nt): 2, Max-length(nt): 4, Min. repeat-times: 8
           -smin: 90, -Smax: 300, -tm: 60, -fl: 150, -l: 500
   1.4

   通过最后的eMapping可以帮助我们评估引物的特异性。打开生成的以.eMap.frg结尾的文件，里面包含了模拟PCR得到的等位基因片段的大小，如：
           >MK1 3 175/176,176/176,177/176
           >MK10 3 264/265,265/265,266/265
   大多数位点可能会得到不止1对等位基因，这是由参考基因组序列组装不完善或者转座元件造成的正常现象。但是，如果某个位点的片段大小存在明显差异，表明发生了非特异性扩增，那么这个位点就应该被舍弃。比如：
           >MK1033 4 220/136,220/136,135/136,136/136
   1.5

   在后缀名为.ssr.mk的文件中可以找到全部引物的信息。然后我们需要在不同参考序列或染色体序列中广泛地选择一些位点，保证微卫星位点尽可能均匀地分布在基因组中，以避免位点连锁。
   1.6

   (可选步骤) 非编码区的微卫星位点多是中性的，多态性相对较高。编码区的微卫星位点与性状关联更强，多态性较弱。对于有参考基因组的物种，我们可以根据研究目的对选择位点所在的基因区域进行检查，以初步筛选。
   1) 可以通过基因组的Feature_table确定SSR标记的位置。Feature_table可以在NCBI特定物种的基因组主页面进行下载。下载后解压并使用Excel打开，根据.ssr.mk文件第一列提供的参考序列名称在Excel中进行查找。然后根据Excel的位点位置进行检查，以确定候选的位点位于编码区还是非编码区。
   2) 也可以通过微卫星数据库 (Avvar et al., 2020) http://data.ccmb.res.in/msdb/ 搜索并下载物种的微卫星的注释信息：Download.tsv。用Excel打开，然后根据.ssr.mk文件中微卫星位点的位置进行检查。
   
2. 微卫星序列的扩增
   
   2.1

   引物合成。微卫星标记的引物通常需要在一条引物的5' 端进行荧光标记。常用的荧光标记物是FAM, HEX, TAMRA, ROX。
   2.2

   进行PCR扩增。使用引物进行普通的PCR扩增即可。荧光引物和PCR产物都应避光保存。
   注：对SSR引物直接进行荧光标记的成本较高，有一些节约实验成本的策略：1) 使用通用荧光引物的三引物反应体系。请读者参见Blacket et al. (2012)。2) 进行多重PCR，即多对引物同时进行扩增。需要进行多次预实验探索引物可行的组合。
3. 微卫星基因分型
   
   3.1

   对于只有基因型定性需求的研究 (例如亲子关系鉴定)，可以使用聚丙烯酰胺凝胶电泳。该方法不需要对引物进行荧光标记。
   3.2

   常用的基因分型方法是基于荧光标记的毛细管电泳，只需要把扩增样品送到提供STR检测服务的生物公司即可。
           为了降低成本，一般我们会要求公司进行混检，即将同一样本不同SSR位点的PCR产物混合进行检测。我们需要提供一个尽可能不会混淆产物识别的混合方案，主要考虑荧光颜色和片段长度两个因素。荧光颜色相同，长度接近的产物不应该混在一起。
           通常20个位点可以被分为3至4组进行检测。如何尽可能地减少分组的数量需要读者根据自己的位点进行综合考虑。如果读者对扩增可能出现的片段长度不了解，建议初次检测时避免高度混合的方案。
   3.3

   生物公司分型后会返回给我们pdf格式的峰图、片段长度读数的表格、以及原始数据文件。将pdf峰图所示等位基因记录在Excel中，数据保存格式见4.1。对于二倍体生物，单峰为纯合子，记录两遍；双峰为杂合子，分别记录，两个峰值之差应为该位点重复单元大小的整数倍。
   3.4

   我们建议使用分型软件，如GeneMarker (Holland and Parson, 2011)，从原始数据对读数重新判读和校正，以尽可能排除无效等位基因并拯救部分缺失值。该软件有一个月的试用期。另外一款可替代的软件是GelQuest。图1中我们给出了微卫星峰图常见的瑕疵类型，以及正确的读数位置供读者参考。对于特定的SSR位点，往往会在不同的等位基因间有稳定的峰形特征。图1中A1、B、D、E等情形，一般能在SSR引物开发验证时发现，可在后续试验中不使用这些分型困难的位点。
           以GeneMarker为例：Open Data-Add-选择同一组的所有文件-OK-Run-左边选择公司使用的内参 (常用ABI或ABI-5-Color)-Next-Next-Ok-Ok。待处理结束后，左边的列表的文件图标会成为绿色。工具栏可以关闭多余窗口。峰图中，左键向右框选为放大、向左框选为返回初始视图，右键拖动。
   
   
   图1. 微卫星峰图常见瑕疵类型 灰色指示了正确的读数位置。(A) 序列端部氢键脱落。应读右侧。这里给出了等位基因邻近 (A1) 和不邻近 (A2) 两种情况。(B)影子带。由PCR扩增时滑链造成，常在左侧出现连续降低的峰，应读右侧最高峰。与邻近杂合双峰 (A1) 的区别特征是杂合峰通常两峰高度接近，且每个峰左侧底部会带一个小的影子峰。(C) 负值峰。软件对来自其它荧光干扰的校正过度。指示值无实际意义，但暗示引物荧光标记脱落或污染。(D) 峰顶部无法对齐到整数。电泳误差，全局统一偏大或偏小读数即可。(E) 单峰竞争优势。邻近的两个峰高差别较大，但弱势峰也表现出微卫星峰的特征 (出现影子带)。建议两个峰都读。
   
   
4. 微卫星数据格式
   
   4.1

   微卫星数据常见的有genepop格式、fstat格式和STRUCTURE格式。其中最主流的是genepop格式。在最初数据手动录入Excel时，推荐两个Excel插件：GenAlEx 6.5 (Peakall and Smouse, 2006) 和Microsatellite toolkit。二选一使用。GenAlEx的输入格式如图2所示。Microsatellite输入格式见插件目录中的MS-data.xls。
   注：GenAlEx导出的genepop格式会默认将等位基因的读数重新编码为两位数，因此不建议使用它的这个功能。
   4.2

   从Excel转换为genepop格式。以Microsatellite为例：打开目录中的MS-tools.xla将插件加载到Excel。Excel加载项-Microsatellites-Microsatellite Toolkit-输入格式选择Diploid_One_column format，不勾选Check-OK-如不筛除种群和位点，就直接点OK-选择Genepop中的3-digitformat，勾选黑色选框还可以计算一些遗传多样性参数-GO！-提示输入一个genepop格式的名字，可以不输入。得到的GenePop3D工作表即为标准的Genepop 3位数字编码的微卫星格式。将它复制并保存到一个新建的文本中即可用于下游分析。
   4.3

   可以使用PGDSpider (Lischer and Excoffier, 2012) 可以将genepop格式转换为其它软件需要的格式。
   
   
   图2. GenAlEx的输入格式 第一行依次为：位点数、个体数、种群数、每个种群的个体数，缺失值表示为0。
   
   
5. 数据质控
   
   5.1

   使用MicroChecker (van Oosterhout et al., 2004) 进行数据质控。首先打开genepop格式文件，如果提示格式错误，请删除文件里POP行后面多余的空格和分隔符。
   5.2

   从左到右依次操作：指定每个位点的重复单元类型-输入预期最大等位基因-检查异常值-分析。注意，MISS值，即000，也会被识别为错误，点击Analyse后忽略即可。
   5.3

   结果会显示等位基因的观测频率和期望频率，下方会提示是否存在stuttering、large allele dropout和 null allleles。
   5.4

   菜单栏的Tools中有两个工具：第一个工具Estimated Frequencies可以给出观测到的等位基因频率和预期的4种算法估计的"真实"频率，并且在Adjusted Genotypes中给出校正的基因型数据。但是这个功能会打乱个体的顺序，因此实际意义有限。第二个工具Null Allele Estimates用于判断位点是否存在无效等位基因。
   5.5

   使用Genepop (Rousset, 2008) 进行哈迪温伯格平衡 (HWE) 检验和连锁不平衡 (LD) 检验。推荐使用Genepop的网页版：Genepop on the web (https://genepop.curtin.edu.au/)。
   1) 点击第一个进入HWE测试。进去后我们选择第一个在种群水平测试杂合子缺失。输出选项中，选择HTML输出。输入文件-浏览-选择genepop文件。然后submit。然后以同样的操作测试杂合子过剩。
   2) 统计测试的原假设为数据服从HWE，即不存在杂合子缺失或过剩，因此结果p < 0.05表明存在杂合子缺失或过剩。由于是多重比较，比较前需要对p值进行Bonferroni多重校正。将得到的p值乘以用于比较的项目数，然后再与0.05显著性水平进行比较。
   3) LD检验的操作也是类似的，但不需要进行Bonferroni多重校正。得到的位点-位点的结果中，如果p值小于0.05，则表明位点存在显著的连锁。
   5.6

   质控可以加强我们对数据质量的了解。但是我们应该谨慎对待和调整我们的数据值。通常数据总是会存在一些瑕疵，如果某些位点有异常，而我们又不愿意剔除这些位点时，我们首先可以回去检查峰图的质量并保守地进行调整，其次是设法给予生物学的解释 (比如杂合子缺失)。HWE偏离在理论上会影响依赖它的软件，如STRUCTURE。另外，增加样本量会改善HWE和LD测试的结果。
6. 微卫星数据载入R语言环境
   R语言提供了丰富的工具包，能够执行几乎所有的微卫星数据分析。常见的工具包比如：hierfstat ( Goudet, 2004)——计算各种统计值 (比如等位基因丰富度)、adegenet (Jombart, 2008) 和pegas (Paradis, 2010)——综合性工具包。读者可以参考Kamvar et al. (2016) 的在线教程https://popgen.nescent.org/，去探索R包中丰富的功能。作为引入，这里我们介绍了将微卫星数据读入R环境的方法。
   
   install.packages("adegenet") ; library("adegenet") #安装并载入包
           #下面给出三种格式的读取例子，载入后的数据类型均为genind。
   
           #载入genepop格式：
           genepopData <- read.genepop("yourdata.gen", # genepop格式文件后缀名应为.gen
           ncode= 3L)#数据中等位基因应该是3位数字编码，否则请设置ncode=2L
   
           #载入STRUCTURE格式：
           structureData <- read.structure("yourdata.str", n.ind=个体数, n.loc=位点数,
           onerowperind=FALSE，#每两行一个个体，否则设置为TRUE
           col.lab=1, #个体名字设置在第一列
           col.pop=2, #种群名字设置为第二列
           ask=FALSE)
   
           #载入fstat格式:
           fstatData <- read.fstat("yourdata.dat")
   
   

致谢

感谢张润志研究员对手稿的审查和指导，感谢中国科学院战略性先导科技专项 (XDA19050204) 提供支持。

参考文献

1. Avvaru, A. K., Sharma, D., Verma, A., Mishra, R. K. and Sowpati, D. T. (2020). MSDB: a comprehensive, annotated database of microsatellites.Nucleic Acids Res 48(D1): D155-D159.
2. Blacket, M. J., Robin, C., Good, R. T., Lee, S. F. and Miller, A. D. (2012). Universal primers for fluorescent labelling of PCR fragments—an efficient and cost-effective approach to genotyping by fluorescence. Mol Ecol Resour 12(3): 456-463.
3. Goudet, J. (2004). HIERFSTAT, a package for R to compute and test hier-archical F-statistics. Mol Ecol Notes 5(1), 184–186.
4. Holland, M. M. and Parson, W. (2011). GeneMarker® HID: a reliable software tool for the analysis of forensic STR data. J Forensic Sci 56(1): 29-35.
5. Jombart, T. (2008). Adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics 24(11): 1403-1405.
6. Kamvar, Z. N., López-Uribe, M. M., Coughlan, S., Grünwald, N. J., Lapp, H., & Manel, S. (2016). Developing educational resources for population genetics in R: An open and collaborative approach.Mol Ecol Resour 17(1): 120-128.
7. Lischer, H. E. and Excoffier, L. (2012). PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28(2): 298-299.
8. Peakall, R. O. D. and Smouse, P. E. (2006). GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes 6(1): 288-295.
9. Paradis, E. (2010). Pegas: an R package for population genetics with an integrated–modular approach. Bioinformatics 26(3): 419-420.
10. Rousset, F. (2008). Genepop'007: a complete re-implementation of the genepop software for Windows and Linux.Mol Ecol Resour 8(1): 103-106.
11. van Oosterhout, C., Hutchinson, W. F., Wills, D. P. and Shipley, P. (2004). MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes 4(3): 535-538.
12. Wang, X. and Wang, L. (2016). GMATA: an integrated software package for genome-scale SSR mining, marker development and viewing. Front Plant Sci 7: 1350.