Methods for Obtaining Full-length DNA Barcodes Using High-throughput Sequencing

研究背景

在过去十年里，随着高通量测序技术的发展，我们见证了生物多样性研究中方法学和应用的巨大转变，例如采用标准脱氧核醣核酸 (Deoxyribonucleic acid，DNA) 序列用于快速准确地物种鉴定；利用高通量测序技术来分析复杂的环境样品 (如混合样品、环境DNA (eDNA)、无脊椎动物来源的DNA (iDNA) 等)。保存有规范记录信息的DNA条形码参考数据库也在全球科研人员的共同努力下逐渐建立起来。生命条形码数据库BOLD (The Barcode of Life Data systems, http://v4.boldsystems.org/) 已经拥有了大约八百万条条形码序列，覆盖了包含动物、植物、真菌等生物在内的约30万种物种 (2020年2月截止)。丰富的参考条形码数据为物种鉴定、系统进化关系的构建、种间交互作用和群落结构的研究，以及加深对生物多样性的理解提供了坚实的基础。
        全球条形码参考数据库已经为多个生态系统的研究提供了重要帮助。早期的生态及生物多样性研究经常利用Sanger测序方法来进行条形码测序，随着高通量测序技术 (High-Throughput Sequencing Technologies，HTS) 的推广，DNA宏条形码和线粒体宏基因组方法也被越来越多的研究者接受并使用。但是当研究某些新的环境样品时，条形码数据库中参考数据的缺乏仍会使得基于HTS的宏条形码组学研究遭遇瓶颈：得到的序列不能被有效分配到具体的物种，从而无法反应真实的生物多样性组成，使研究者难以深入了解环境的生态关系。
        以昆虫CO1条形码为例，标准条形码获取成本较十多年前显著降低。在不包括样品收集及处理的成本情况下，传统的Sanger测序的平均生产成本约10美元，如果要构建1亿个样品的条形码数据集，则需要10亿美元。人类基因组在最初构建时，大约花费了30亿美元，但随着高通量技术的发展，目前一个基因组的测序成本已经降低为不到600美元。
        先前基于高通量测序技术获取DNA条形码的方法也存在着各自的优缺点 (图1)，比如：基于罗氏454测序平台构建单样品条形码的研究可以通过拼接获得CO1全长条形码，但成本较高；测序读长较短的部分Illumina测序平台也逐渐用于条形码数据的获取，有研究者使用长度为313 碱基对 (base pair，bp) 的CO1条形码进行分析，该方法虽然可以一次性测序获得全长序列，但由于条形码长度较短，因此特异性会显著降低；也有研究者利用两次聚合酶链式反应 (polymerase chain reaction，PCR) 增加标签序列和测序引物，继而通过测序来直接获取全长条形码，但相关操作比较繁琐。在基因组领域，短数据可以利用位置关系被拼接为准确的长片段；根据类似的原理，所在项目组已经专门开发过针对条形码的组装算法，可以被用于组装混合样品中的全长条形码。根据当时的研究背景，我们开发了基于Illumina Hiseq平台和PacBio 平台的HIFI-Barcode方法，可以通过编码引物的方法，得到带有特定标签序列的PCR产物，从而实现一次性获得96个样品的DNA条形码的目标。该方法发表于Liu et al., 2017。HIFI-Barcode 方法有很高的准确性和效率，但是涉及到较为复杂的分析过程，所以我们认为基于高通量的条形码技术仍有提升的空间。随着国产测序仪MGISEQ 2000 推出单端400 bp测序的模块，我们进一步在此平台上开发出了HIFI-SE的方法，可以不需要进行打断，在其后的数据分析阶段，只需要进行扩增子两端序列的简单拼接，就可得到完整的COI条形码序列。该方法发表于Yang et al., 2020。
        以上三种方法相互独立，测序及分析策略不同，但是在前期样品准备，DNA提取和PCR操作等方面基本一致，所以在此操作手册中，我们将三种方法的实施步骤及相关比较纳入其中，方便研究者选择合适的建库测序平台和后续的分析方法。


图1. 高通量获取标准DNA条形码的策略比较，不同实验方案和测序策略 *代表本研究研制的方法策略。图片修改自Yang et al., 2020。A-E 分别代表的研究为 Meier et al., 2016; Shokralla et al., 2015; Cruaud et al., 2017; Liu et al., 2017; Hebert et al., 2018; Yang et al., 2020。

材料与试剂

1. 96孔PCR板 (Axygen，PCR-96M2-HS-C)，1.5 ml离心管 (Axygen)，2.5 μl，100 μl，200 μl，1 ml枪头 (Axygen Filter Tips)
2. 96方孔深孔板 (2 ml，PP MASTERBLOCK®, 96 Well)
3. 封口膜 (AxysealTM sealing film)
4. 琼脂糖 (NET WEIGHT, REGULAR AGAROSE G-10)
5. exTaq DNA聚合酶 (TaKaRa Ex Taq™)
6. dNTP mix (TaKaRa Ex Taq™)
7. 10x primer buffer (TaKaRa Ex Taq™)
8. ddH2O (广州誉维生物科技仪器有限公司Unique超纯水机制)
9. 合成引物 (上海生工)
10. 昆虫裂解液 Insect Lysis Buffer (见溶液配方)
11. 吸附缓冲液Binding Buffer (见溶液配方)
12. 洗脱液Wash Buffer (见溶液配方)
13. 吸附混合液Binding Mix (见溶液配方)
    注：对于DNA提取可根据自身实验室选择合适的方法和试剂主要试剂，本方法主要介绍玻璃纤维板DNA提取法 (Glass Fiber Plate method) 用昆虫腿提取DNA的方法。所用主要试剂见10-13。
    

仪器设备

1. 0.2-2 μl, 2-20 μl, 20-200 μl, 100-1000 μl移液枪 (Eppendrof)
2. 台式高速冷冻离心机 (Beckman，AllegraTM 25R Centrifuge)
3. 水浴锅 (DK-8D型，上海精宏实验设备有限公司) 或者恒温震荡仪 (MS-100 THERMO-SHAKER)，用于组织裂解
4. 96孔PCR仪 (Thermo Fisher)
5. 电泳仪、电泳槽 (DYY-6C型，北京市六一仪器厂)
6. 凝胶成像系统 (BIO-RAD)
   

软件

1. HIFI-barcode-hiseq，适用于Illumina ，MGISEQ等二代测序平台https://github.com/comery/HIFI-barcode-hiseq。
2. HIFI-barcode-pacbio，适用于Pacbio测序平台https://github.com/comery/HIFI-barcode-pacbio。
3. HIFI-barcode-SE400，适用于MGISEQ 2000 SE400 测序平台https://github.com/comery/HIFI-barcode-SE400。

实验步骤

本实验流程包含了三种高通量条形码获取的流程，三种方法的测序前实验环节方法一致 (对应步骤1至3)，测序文库构建、测序及分析环节有差异。在完成PCR产物的混合后，可选择一种平台进行后续的建库测序。流程主要如图1所示，方法1 HIFI-Barcode-Hiseq/MGISEQ可选择平台及测序技术为Hiseq或MGISEQ PE150；方法2 HIFI-Barcode-Pacibio可选择平台 Pacbio；方法3 HIFI-Barcode-SE400可选择MGISEQ SE400。后续根据标签数目的增加，单条形码的测序成本还可以进一步降低。

表1. 三种方法的多维度比较


图2. 实验操作及分析流程图

1. 动物全基因组DNA提取
   
   1.1

   简介
   本流程采用利用玻璃纤维板DNA提取法 (Glass Fiber Plate method) 进行昆虫全基因组DNA的提取。该方法是加拿大DNA条形码中心提供的批量提取DNA的方法，研究者可以根据自身项目需求采取不同的基因组DNA提取方法。
   
   1.2

   主要流程
   

   1)

   获取昆虫组织样品：
   a)利用酒精清洁实验台，并准备好相关材料和试剂。
   b)在96孔板的每个孔中加入30 μl的无水乙醇防止昆虫组织因为静电蹦跳。
   c)利用干净的镊子夹取足量昆虫组织按照顺序放入96孔板中。
   d)取样后加上封盖密封，进行短暂离心；第五：离心后打开盖子，并37 °C加热蒸发酒精，获得用于提取DNA的昆虫组织样品。

   2)

   昆虫组织裂解：在带裙边的96孔板的每个孔中加入50 μl昆虫裂解液，并将1) 中获取的组织样品依次添加到添加了裂解液的孔中，在56 °C温浴6 h使样品充分裂解。

   3)

   裂解后，首先将板子进行离心 (转速为1500 × g, 15 s) 以移除盖子上的冷凝水。然后，加100 μl的binding mix溶液到每个孔中，1000 × g震荡离心20 s。

   4)

   换板：移除封盖，将板中溶液加入到玻璃纤维板 (简称GF板) 中，GF板放在方形盒上，用密封膜封盖GF板。

   5)

   吸附：7000 × g离心5 min，基因组DNA将会被吸附在GF板的玻璃纤维膜 (简称GF膜) 上。

   6)

   两次洗脱：第一次：每孔加入180 μl的PWB洗脱液，用新密封膜密封，5000 × g离心2 min。第二次：每孔加入750 μl的WB洗脱液，重新用密封膜密封，5000 × g离心5 min。

   7)

   去除残留酒精：揭开封盖后，将GF板放置在56 °C中温浴30 min。

   8)

   孵育：在GF板下放置收集板；56 °C预热ddH2O后添加45 μl到每个孔的GF膜上，室温孵育1 min，并密封。

   9)

   收集DNA：封装好的板子放到干净的方孔板并5000 × g离心5 min；移走并丢弃GF板；DNA已经提取并收集好。
2. 条形码序列扩增与检测
   
   2.1

   引物标签设计
   

   1)

   标签 (index) 设计考虑的因素：由于需要使用高通量测序技术对昆虫条形码COI序列进行测序和分析，因此在设计引物标签时需要综合考虑以下因素：a) 标签可识别性：标签差异碱基数目、标签数量等；b) 成本相关因素：测序读长、序列长度、标签长度；c) 扩增效率：是否会影响扩增效率，是否容易造成非特异性扩增等。

   2)

   标签设计：首先，选定标签长度：利用Barcrawl程序 (v.30May2012) 生成5-10 nt长度的标签若干，在保证碱基差异大于或等于2情况下且标签数目大于或等于96条的情况下，对标签数目和长度进行比较，后选定为5 nt长的标签序列集合。对应命令行为"barcrawl -l 5 -m 2 -g 70 -a 5 -b 5 -p 5 -c 30 –nd''，命令行对应的主要条件为：两标签序列之间差异大于或等于2，GC含量在30%-70%之间。共生成了160条标签，从生成的标签集合中随机选取96条作为后续研究使用的引物标签 (附录A表6)。
   2.2

   标签引物合成
   按照标签表 (附录A表6) 设计出含标签序列的条形码引物：在标准引物的上 (下) 游的5'端添加正向 (反向) 标签序列，形成新的带有标签序列的96对COI基因标准引物。加标签后的引物由上海生工有限公司合成。
           标准引物序列为：
           上游引物LCO1490：5'-TAAACTTCAGGGTGACCAAAAAATCA-3'。
           下游引物HCO2198：5'-GGTCAACAAATCATAAAGATATTGG-3'。
           例如：编号001的引物对序列为：
           上游引物LCO1490：5'-AAAGCTAAACTTCAGGGTGACCAAAAAATCA-3'。
           下游引物HCO2198：5'-AAAGCGGTCAACAAATCATAAAGATATTGG-3'。
   2.3

   PCR
   

   1)

   实验组：在96孔PCR板中依次加入带有标签的引物、昆虫DNA和扩增需要的试剂；每个孔包含1个PCR扩增体系 (表2)。
   
   表2. PCR体系
   

   组分	每个PCR体系加入量
   ddH2O	16.2 μl
   10 μM 正向标签引物	1 μl
   10 μM 反向标签引物	1 μl
   10x primer buffer	3 μl
   dNTP mix	2.5 μl
   exTaq	0.3 μl
   DNA样本或无菌水	1 μl (1 μl DNA样本中的DNA含量约200 ng)

   
   

   2)

   空白对照：95个昆虫DNA分别对应96孔板的95个孔，剩余的1个作为空白对照，空白对照孔中需要将反应体系中的1 μl DNA替换为1 μl ddH2O，其余成分不变。

   3)

   反应条件：反应条件与常规COI条形码PCR反应一致，具体为：94 °C预变性1 min，然后进入5个循环：94 °C变性30 s、45 °C退火40 s、72 °C延伸1 min，5个循环反应结束后进入35个循环：94 °C变性30 s、51 °C退火40 s、72 °C延伸1 min，循环结束后，72 °C延伸10 min，12 °C 保持。
   2.4

   电泳，分装与质量检测
   

   1)

   电泳：PCR扩增完成后，电泳检测PCR扩增产物。

   2)

   分装：从PCR板的每个孔 (对应一个PCR反应产物25 μl) 中吸取5 μl进行混合，获得一份480 μl的混合结果后，再将混合物均分为5份储存，每份包含96 μl的混合产物。

   3)

   分别选取一份扩增混合产物送去进行DNA质量检测，包括Qubit检测和电泳检测，并确认是否能够进行建库测序。
3. 高通量测序技术选择与测序
   

   得到PCR产物之后，可以根据选择的方法对测序平台进行选择，然后进行相应的标准文库的构建和测序。

   3.1

   HIFI-Barcode-Hiseq/MGISEQ
   PE150测序。PCR产物质量检测合格后，选取一份混合产物送至测序公司进行建库与测序。建库插入片段为250 bp小片段。建库后，使用Illumina Hiseq 4000 (或BGISEQ系列平台) 进行双端测序，测序读长为150 bp，数据量约2-5 Gb。插入片段长度和测序读长为固定参数，不能随意修改。
   
   3.2

   HIFI-Barcode-PacBio
   PacBio测序。选取一份扩增产物混合液进行PacBio RSII平台质量检测、建库与测序，数据量约2-5 Gb。
   注：由于三代测序所需的DNA质量和总量都比较高，所以三代测序可能需要3-4管 (96 μl/管)。
   3.3

   HIFI-Barcode-SE400
   SE400测序。选取一份扩增产物混合液进行MGI 2000 平台的SE400平台建库与测序，数据量约2-5 Gb。
   

测序数据分析

1. 简介
   针对相同的PCR产物，在不同的平台上进行建库测序，最终根据数据特点的不同，我们设计了对应的三种分析程序来还原DNA条形码序列。
2. 基于二代测序数据的条形码获取- HIFI-Barcode-Hiseq/MGISEQ
   获取每个PCR板的二代测序数据后，分别对两份数据进行分析。首先进行原始数据的过滤，再进行条形码数据的组装与分析；条形码组装的详细流程见图3。为了方便使用，已将流程打包为软件包HIFI-barcode (https://github.com/comery/HIFI-barcode-hiseq)。
   
   
   图3. HIFI-Barcode方法条形码组装流程 HIFI-Barcode程序包含数据过滤和拆分、排序和聚类、双端拼接、间隙填充等步骤。
   
   
   2.1

   数据过滤
   按照以下条件对原始测序数据进行质量控制和数据过滤：
   

   1)

   去除有接头污染的序列：至少能比对到接头15碱基且最多有3个错配。

   2)

   去除N数目大于10 的序列。

   3)

   去除低质量reads：30%的碱基质量值低于33的序列。
   2.2

   序列拆分和过滤
   

   1)

   拆分：由于96个样品分别对应一对已知的标签序列，因此可以结合测序序列双末端所包含的标签序列和引物序列的信息，将所有reads序列进行拆分。可以分为三类：带有96种前端 (条形码5'末端) 标签引物的reads，带有96种尾端 (条形码3'末端) 标签引物的reads和没有引物序列的reads。

   2)

   去冗余：将拆分得到的含有标签引物的单端序列按照98%的相似度聚类，减少冗余数据。

   3)

   丰度排序与过滤：根据2) 的聚类结果将序列按照丰度高低进行排序，并记录相关丰度信息；根据丰度信息保留三类序列进行下游分析：a) 保留聚类结果中最高丰度的序列作为昆虫目标条形码序列进行后续的连接和组装；b) 与最高丰度序列差异大于2%的序列也被保留，作为候选结果；c) 选择聚类结果中丰度不小于相同标签的最高丰度的1/10的序列，作为疑似共生物 (如潜在的寄生虫、沃尔巴克氏体感染或肠道微生物等) 条形码序列进行后续分析。
   2.3

   双端拼接
   

   1)

   拼接：首先，根据reads的双端信息把属于同一条片段的两端序列挑选出来；然后，将每一对PE reads按照以下条件拼接起来：a) 序列之间具有高于95%的相似性重叠区；b) 重叠区域达到一定的重叠长度：前两类有引物信息的reads拼接长度需要在130-270 bp间，第三类无引物信息的reads拼接后的序列长度限制在170-270 bp之间。

   2)

   深度过滤：统计每种标签对应的拼接序列的各个碱基的深度，将覆盖深度不足5的碱基删除；从而获得可以代表标签对应的昆虫条形码前后两端的一对序列。
   2.4

   间隙 (gap) 填充
   

   1)

   输入数据：带标签的双端拼接结果作为两端信息输入，中间无标签的拼接结果统一作为中间信息输入。

   2)

   组装算法：采用SOAPBarcode算法 (图4) 对条形码序列进行填充，从而获得高精度的完整的昆虫COI基因条形码。a) 每一对两端数据的前端为起始点，末端为终点，并用kmer构建de bruijin图，从起点到终点查找潜在的连接路径；b) 采用以下策略保证连接正确性：第一，删除在分叉处kmer丰度小于平均kmer丰度10%的路径 (蓝色)；第二，如果经过第一步过滤后还有多条出度 (out degree，OD) 存在，则针对不同ODs和位于最后分叉前的kmer间进行reads数目统计，后删除丰度小于平均reads丰度10%的ODs (黑色)；第三，删除超出了预先设定的长度的路径 (红色)。
   
   
   图4. SOAPBarcode组装算法(引自Liu et al., 2013)
   

   3)

   组装结果输出：每个标签选取分值最好且长度与标准条形码相似的序列作为对应昆虫的条形码输出；根据前期丰度过滤后作为备选序列的双端序列所填充出来的条形码则作为非目标条形码输出。
   

   命令行：
   $ python3 HIFIBarcode.py all -outpre hifi -index 5 -q1 test1.fq.gz -q2 test2.fq.gz -primer indexed_primer.txt
   注：“indexed_primer.txt”文件详见附录B
   综上，可获得结合了二代测序技术和生物信息分析技术得到的目标和非目标条形码。

3. 基于三代测序数据的条形码获取- HIFI-Barcode-PacBio
   三代测序技术具有较长的读长，可以直接获得全长条形码序列，但由于三代测序技术的错误率较二代测序更高，因此需要对数据进行比较过滤后再获取全长条形码数据。分析方法如下：
   
   1)

   提取一致性循环序列 (consensus circular sequence，CCS)：从测序公司获得过滤后的h5格式文件后，直接使用Pacbio官方提供的流程PacbioSmartAnalysis pipeline将下机文件转换成FASTA格式。
   2)

   丰度过滤：三代测序结果中条形码全长可以被完全测通，且会产生多个循环，因此，本环节将循环数小于15次的序列直接过滤。
   3)

   引物序列定位：利用动态规划算法与通用引物的序列，定位到引物序列对应的位置，以此定位到条形码的双端；由于三代测序错误率较二代测序更高，我们在引物定位环节允许2个碱基的错配和1个碱基的插入或缺失。
   4)

   标签序列定位与样品匹配：将通用引物两端各自向外延伸5 bp，获得标签序列。由于标签序列之间至少有2个差异，因此若标签序列中有连续4 bp或以上的序列能够匹配到某一样本对应的标签序列，则判定该序列可能为对应样品条形码。
   5)

   确定样品条形码：a) 由于存在PCR扩增错误和PacBio测序错误，因此需要选取循环数交高的序列；b) 循环数最高且能够成功翻译为蛋白的序列作为目标条形码序列输出；c) 其它输出的条形码作为非目标条形码输出。
   命令行：
   # step1, 从H5文件中提取CCS序列
   $ source /path/PicBio/smrtanalysis/current/etc/setup.sh
   $ fofnToSmrtpipeInput.py my_inputs.fofn > my_inputs.xml
   $ smrtpipe.py --params=settings.xml xml:input.xml
   # step2, 从H5文件中提取序列循环数信息
   $ python bin/ccs_passes.py data/*.ccs.h5 >ccs_passes.lst
   # step3, 通过循环数信息过滤低质量序列
   $ awk '$2>=15{print $1}' ccs_passes.lst >ccs_passes_15.lst
   $ perl ./bin/fish_ccs.pl ccs_passes_15.lst
   data/reads_of_insert.fasta >ccs_passes_15.fa
   # step4, 根据引物序列拆分序列到对应样品
   $ perl ./bin/1.primer_like_extract.pl -p experiment_data/primer.fa -index experiment_data/index.xls -fa ccs_passes_15.fa -cm 2 -cg 1
   # 通过聚类得到丰度最高的目标序列
   $ cd 02.assignment/
   $ perl ../bin/2.cluster_count_passes_length.pl -ccs ccs.successfully_assigned.fa -pattern check.ccs_passes_15.fa.log -passes ../ccs_passes.lst
   $ perl ../bin/change_name-location.pl cluster.top1.fas >hifi-barcode-pacbio.cluster.top1.fa
   *具体输入文件格式见github页面：https://github.com/comery/HIFI-barcode-pacbio
4. 基于SE400测序数据的条形码获取-HIFI-Barcode-SE400
   SE400 平台是数据可以应用Python包HIFI-SE 来完成分析。
   
   
   图5. HIFI-SE 数据分析流程
   
   
   1)

   数据过滤，去掉包含M个N的序列，删除expected error > 10 (expected error 定义见Yang et al., 2020)的序列，或者可以根据quality的参数设置阈值。
   2)

   序列拆分，由于96个样品分别对应一对已知的标签序列，因此可以结合测序序列双末端所包含的标签序列和引物序列的信息，将所有reads序列进行拆分。主要可以分为两类：带有96种前端 (条形码5'末端) 标签引物的reads，带有96种尾端 (条形码3'末端) 标签引物的reads。其他未分配的reads主要是因为测序错误引起的引物或标签序列的错配情况，或者是PCR产生的嵌合体。
   3)

   组装，可分为两种模式，cluster模式和consensus模式。在cluster模式下，对应每个单独的样品，会对5'和3'的序列分别进行聚类，然后选取丰度最高的两个cluster序列，根据前后端的overlap信息将其拼接起来。而在第二种模式下，不对序列进行过滤，而是用所有序列直接生成consensus序列，用前后端分别得到的consensus序列直接拼接。一般来讲，如果一个样品中target物种的丰度很高且没有污染，可以考虑用第二种模式进行组装，速度会比较快。但是如果孔增产物比较复杂时，建议使用第一种模式进行组装。
   4)

   物种鉴定，HIFI-SE程序包含了一个taxonomy的程序可以允许在本地进行COI序列的物种鉴定，具体实现是将序列提交到BOLD服务器进行物种鉴定，然后在返回结果显示。可以极大的提高物种注释的效率。
   命令行：
   $ python3 HIFI-SE.py all -outpre hifi -trim -e 5 -raw test.raw.fastq -index 5 -primer index_primer.list -mode 1 -cid 0.98 -oid 0.95 -seqs_lim 50000 -threads 4 -tp 2

溶液配方

表3. DNA提取试剂、材料

表4. 预备液配置表

表5. DNA提取试剂

致谢

本研究获得科技部科技基础资源调查专项《中国东部传粉昆虫资源调查与评估》 (2018FY100403) 以及深圳市科创委基金 (NO. JCYJ20170817150755701) 的资助

竞争性利益声明

无经济或非经济性竞争性利益。

参考文献

1. Ivanova, N. V., Dewaard, J. R. and Hebert, P. D. N. (2006). An inexpensive, automation‐friendly protocol for recovering high‐quality DNA. Mol Ecol Notes 6(4): 998-1002.
2. Liu, S., Li, Y., Lu, J., Su, X., Tang, M., Zhang, R., Zhou, L., Zhou, C., Yang, Q., Ji, Y., Yu, D. W. and Zhou, X. (2013). SOAPBarcode: revealing arthropod biodiversity through assembly of Illumina shotgun sequences of PCR amplicons. Methods Ecol Evol 4(12): 1142-1150.
3. Liu, S., Yang, C., Zhou, C. and Zhou, X. (2017). Filling reference gaps via assembling DNA barcodes using high-throughput sequencing—moving toward barcoding the world. GigaScience 6(12): gix104.
4. Yang, C., Zheng, Y., Tan, S., Meng, G., Rao, W., Yang, C., Bourne, D. G., O’Brien, P. A., Xu, J., Liao,, S., Chen, A., Chen, X., Jia, X., Zhang, A. and Liu, S. (2020). Efficient COI barcoding using high throughput single-end 400 bp sequencing. BMC Genomics 21(1): 1-10.