减法基因组学的新序列发现

摘要

该协议的目的是使用计算和工作台研究的组合来寻找无法轻易从共同净化序列中分离的新序列, 这可能只是部分已知的序列。

摘要

减法基因组学可以用于任何研究, 这些研究的目标是确定嵌入在更大基因组背景下的基因、蛋白质或一般区域的序列。减法基因组学使研究人员能够通过全面测序和减去已知的遗传元素 (参考, r) 来分离感兴趣的目标序列 (t)。该方法可用于识别线粒体、叶绿体、病毒或种系受限染色体等新序列, 在 t 不能轻易从 r 中分离时尤其有用, 该方法从综合基因组数据 (r + t) 开始, 该方法使用基本局部对齐搜索工具 (blast) 对参考序列或序列, 删除匹配的已知序列 (r), 留下目标 (t)。为了使减法工作得最好, r 应该是一个相对完整的草案, 缺少 t。由于减法后剩余的序列是通过定量聚合酶链反应 (qpcr) 进行测试的, 因此 r 不需要完整就可以工作。在这里, 我们将计算步骤与实验步骤链接到一个周期, 可以根据需要迭代, 按顺序删除多个引用序列, 并完善对 t 的搜索。减法基因组学的优点是, 即使在物理纯化困难、不可能或昂贵的情况下, 也可以识别出一个全新的目标序列。该方法的一个缺点是为减法寻找合适的参考, 并获得 t 阳性和阴性样品进行 qpcr 检测。我们描述了该方法在斑马雀的细菌系限制染色体中识别第一个基因的应用。在这种情况下, 计算过滤涉及三个引用 (r), 在三个周期内按顺序删除: 不完整的基因组组装、原始基因组数据和转录数据。

引言

这种方法的目的是从基因组上下文中识别一个新的靶向 (t) 基因组序列, 无论是 dna 还是 rna, 还是参考 (r) (图 1)。如果目标不能物理分离, 或者这样做的成本很高, 则该方法最有用。只有少数生物有完美完成的减法基因组, 因此我们方法的一个关键创新是将计算方法和板凳方法结合成一个周期, 使研究人员能够在参考不完善的情况下分离目标序列, 或者是草案来自非模型生物体的基因组。在周期结束时, qpcr 测试用于确定是否需要更多的减法。经过验证的候选 t 序列将显示 qpcr 在已知 t 阳性样本中的统计上更强的检测结果。

该方法已在发现没有宿主同源的 1^、2、3、4的新细菌药物靶点中得到应用, 并从受感染的宿主^{中识别新病毒5,6}。除了 t 的鉴定外, 该方法还能提高 r: 我们最近用这种方法从斑马雀参考基因组中鉴定了936缺失的基因, 并从仅有细菌的染色体 (t)⁷中发现了一个新的基因。当 t 可能与已知序列有极大的差异时, 或者当 t 的同一性基本不明确时, 例如在斑马雀细菌群限制的^{第7号染色体中}, 减法基因组学尤其有价值。

通过事先不需要 t 的阳性识别, 减法基因组学的一个关键优势是它是无偏的。在最近的一项研究中, readhead等人研究了阿尔茨海默氏症与四个大脑区域病毒丰度之间的关系。为了病毒识别, readhead等人建立了一个由515种病毒⁸组成的数据库, 严重限制了他们的研究能够识别的病毒制剂。减法基因组学可以用来比较健康和阿尔茨海默氏症的基因组, 以分离可能与疾病相关的新病毒, 而不管它们与已知的传染因子有何种相似之处。虽然已知有263种人类目标病毒, 但据估计, 约有167万种未被发现的病毒物种存在, 其中 631.00-827 000 种有可能感染人类⁹。

分离新病毒是减法基因组学特别有效的一个领域, 但一些研究可能不需要如此严格的方法。例如, 识别新病毒的研究使用了无偏见的高通量测序, 然后是逆转录酶和 blastx 的病毒序列⁵或丰富病毒核酸提取和逆转录病毒序列⁶. 虽然这些研究采用了重新测序和组装, 但没有使用减法, 因为目标序列是通过 blast 得到肯定的。如果这些病毒是完全新颖的, 与其他病毒没有关系 (或远近相关), 减法基因组学将是一种有用的技术。减法基因组学的好处是可以获得全新的序列。如果生物的基因组是已知的, 它可以被减去留下任何病毒序列。例如, 在我们发表的研究中, 我们通过减法基因组学从斑马雀中分离出一个新的病毒序列, 尽管这不是我们最初的意图^.

在抗生素耐药1、^2、3、4急剧上升的推动下, 减法基因组学也被证明在识别^{细菌疫苗靶点}方面很有用。为了最大限度地降低自身免疫反应的风险, 研究人员通过减去人类宿主中具有同源性的任何蛋白质来缩小潜在的疫苗目标。一项特别的研究, 研究了伪结核病的结核菌, 对几个细菌基因组的脊椎动物宿主基因组进行了减法, 以确保可能的药物靶点不会影响宿主中的蛋白质, 从而产生副作用¹. 这些研究的基本工作流程是下载细菌蛋白质组, 确定生命蛋白质, 去除多余的蛋白质, 使用 blastp 分离必需的蛋白质, blastp 对宿主蛋白质组分离任何具有宿主同源体的蛋白质^{1,2,3 个,4}. 在这种情况下, 减法基因组学确保所开发的疫苗不会对宿主^{1、2、3、4}产生任何非目标影响。

我们使用减法基因组学来识别细菌系限制染色体 (grc) 上的第一个蛋白质编码基因 (在这种情况下, t), 它存在于种系中, 但不存在于两性的体细胞组织10中.在这项研究之前, 唯一了解 grc 的基因组信息是一个重复的区域¹¹。在卵巢的 rna 序列中进行了 de novo 的组装, 从成年斑马雀的卵巢和睾丸组织 (r + t) 中进行。序列的计算消除使用已公布的体细胞 (肌肉) 基因组序列 (r 1) 12、其原始 (somatic) 读取数据 (r2) 和体细胞 (大脑) 转录组 (r3)¹³进行.每个周期的第5步 (图 2a) 的 qpcr 测试可连续使用三个引用, 这表明需要额外的筛选。通过 dna 和 rna 的 qpcr 以及克隆和测序, 证实了发现的α-snap 基因。我们在示例中表明, 此方法是灵活的: 它不依赖于匹配的核酸 (dna 与 rna), 减法可以使用由组件或原始读取组成的引用 (r) 进行。

研究方案

1.德诺沃组装启动顺序

注: 只要可以从这些数据生成程序集, 就可以使用任何下一代序列 (ngs) 数据。合适的输入数据包括 illumina、pacbio 或 oxford nanopore 读取组装成一个快速文件。对于具体性, 本节描述了一个基于照明的转录组件, 特定于我们所做的斑马雀研究⁷;但是要知道, 具体细节将因项目而异。在我们的示例项目中, 原始数据来自 miseq, 并且从每个样本中获得了大约 1, 000万个配对读取。

1. 使用三叉戟 0.32^{14 去除}illumina 适配器和低质量的碱基。在命令行上, 输入:
   java-jar 装饰-0.32.jar pe-电话 33 forward.fq.gz reverse.fq.gz baseout 质量 _ 和 _ 适配器 _ 裁剪 illumipleip:useq3-pe:2:30 leading:3 slinging:3 slingingw:4leenle:40
2. 使用 pear¹⁵ v. 0.9.6 创建高质量的合并读取从修剪输出配对读取, 使用默认参数。在命令行上, 输入:
   梨-f & lt;quality_and_adaptor_trimmed_1P.fastq>-r & lt;quality_and_adaptor_trimmed_2P.fastq>
3. 使用 rep从而在 1.1^{16 上}对通过 pear 生成的读取进行错误更正。遵循^{17 中}描述的分步协议。
4. 在默认模式下使用 trinity v.^{2.4.0 18}来组装更正后的序列。对于字符串特定的库, 请使用-ss _ lib _ type 参数。输出是一个快达文件 (您的 _ 目. fasta)。在命令行上, 输入:
   trinity-seqtype fq-ss _ libn _ likl _ livcl _ lise10g–输出 trinity _ 输出-左质量 _ and 适配器 _ 跌幅 _ 前进 _ p现度 _ reads. fq-右质量 _ and _ 适配器 _ 削减 _ 反向 _ pdq/cpr 10
   注: 输出将放置在一个新的目录中, trinity _ 输出, 程序集将被命名为 "trity-f. fasta", 如果需要, 可以将其重命名为您的 _ 程序集. fasta。有关更多详细信息, 请参阅三一重工网站: https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity。

2. 根据参考序列进行装配

注: 当引用是程序集或长读取 (如 sanger) 时, 请使用此步骤;如果它是由原始 ilumina 读取组成的, 请参阅下面的步骤 3, 以便将读取映射到查询。所有 blast 步骤都是通过版本 2.2.29 + 完成的, 尽管这些命令应适用于任何最新的 blast 版本。

1. 在命令行中创建参考序列的 blast 数据库 (核苷酸 _ 引用. fasta)。在命令行中输入以下内容:
   核苷酸中的核子。
2. blastt 将查询程序集 (在步骤1中生成) 与引用数据库匹配。若要获取输出文件, 请使用 [out _ 成果. txt]并生成表格输出 (使用 python 脚本执行后续处理步骤所需), 请使用 [-outfmt 6].这些选项可以按任意顺序组合, 因此一个示例完整命令是 [blastn-查询您的 _ 程序。如果需要电子值设置, 请使用具有适当数字的-evalue 选项, 例如 [-evalue 1e-6]。但是, 请注意, 减法周期有效地反转了讨论中描述的值设置。
3. 为了提高严格性, 请使用程序集中的蛋白质序列作为 blast 查询, 并将其转换为核苷酸 blast (tblastn), 用于6路转换 (核苷酸) 数据库。这种方法适用于大多数非模型系统, 避免了蛋白质注释不完全的问题。
   1. 使用-db _ gencode 选项, 确保为所研究的生物体选择正确的遗传代码。若要获取查询的蛋白质序列, 请运行 trandecoder. longorfs 命令 (从 trandecoder 包 v. 3.0.1), 以确定从组装的查询序列中最长的开放读取帧。命令是 [trandecoderd. longorfs-t 您的 _ 目. fasta];输出将被放置在名为 "转录. transdecder _ dir" 的目录中, 并将包含一个名为 "最长 _ orfs. pep" 的文件, 其中包含您的 _ 集合. fasta 中每个序列中最长的预测蛋白质序列。
   2. 若要使用 tblastn, 请运行命令 [tblastn 查询最长 _ orfs. pep-db 核苷酸 _ 公核苷. db out blast _ ft. txt-输出 6].如果有高质量的蛋白质参考, 使用蛋白质与 blastp 匹配, 而不是 tblastn。
   3. 建立一个关于蛋白质参考的 blast 数据库 [makeblastdb-db 类型 prot in 蛋白 _ 参考. fasta out 蛋白 _ 参考. db], 然后 [blastp-查询最长 _ 只能 _ orfs. pep-db 蛋白 _ 参考. db out 为主。请确保将结果另存为用于下游处理的文件, 并使用表格 (输出 6) 以确保 python 脚本可以正确地分析它们。

3. 地图读取到大会上

注: 如果参考数据集包含原始基因组读取, 而不是组装序列或 sanger 序列, 则可以使用此方法 (步骤 2.1)。

1. 使用 bw7-mem v 0.7.12¹⁹或 bowtie220^,将下载的原始读取 (raw _ reads. fastq) 映射到查询程序集中。输出将是. sam 格式。命令如下: 第一个索引程序集: [bwa 索引您的 _ 册. fasta], 然后映射读取 [bwa mem 您的 _ 册.fasta. 汇编. 请 > mopem. sam]。(请注意, 此处的 ">" 符号不是大于符号; 而是指示输出进入文件 mappit. sam)。

4. 使用 python 脚本删除任何匹配的序列

注: 提供的脚本使用 python 2.7。

1. 在步骤2之后, 使用减法 python 脚本, 使用命令 [./Non-matching_sequences.py您的 _ 程序. 请快速编写 blast _ 构开除. txt]。在运行脚本之前, 请确保 blast 输出文件的格式为 6 (表格)。该脚本将以名为您的 _ 集合. fasta _ 不匹配 _ 配合 _ 小时候的匹配序列 _ blast _ results.txt.fasta 和记录的匹配序列 (如您的 _ 程序集. fasta _ 匹配 _ 小时候 _ 生于 _ blast _ ast _ ast _ ast _results.txt.fasta. 非匹配文件将是最重要的, 它是潜在的 t 序列的来源, 用于测试和减法基因组学的进一步循环。
2. 在步骤3之后, 运行 python 脚本 unmappen. py 以输入步骤3.1 中的. sam 作为输入, 并标识查询序列的名称, 而不进行任何匹配的读取, 并将它们保存到新的文本文件中。使用命令 [./删除 unppol.py. mappd. sam], 输出将被映射. sam _ contigs _ 连同 _ no _ reads. txt。(程序将生成一个瘦身 sam 文件, 删除所有未映射的读取; 出于本协议的目的, 可以忽略此文件, 但可能对其他分析有用)。
3. 由于上一步的输出是一个名为 mappd. sam _ contigs _ 与 _ no _ re10s. txt 的文本文件中的序列名称列表, 用这些序列提取一个快速的文件: [./getcontig. py 您的 _ 装配。输出将是一个名为 mappd. sam _ contigs _ 与 _ no _ reads.txt.fasta 的文件。

5. 为保留的序列设计引物

注: 此时有一个包含候选 t 序列的快达文件。本节介绍 qpcr, 以实验测试它们是来自 t 还是来自先前未知的 r 区域。如果步骤4中的减法删除了所有序列, 则初始组件未能包含 t, 或者减法可能过于严格。

1. 手动使用 geneious²¹来确定最佳的引物序列。
   1. 突出显示前馈引物的候选序列 21-28 bp。避免运行4个或更多的任何基地。尝试以所有底膜的组合相当均匀的区域为目标。在 3 ' 端的单个 g 或 c 是有益的, 有助于锚定底漆。
   2. 单击屏幕右侧的"统计" 选项卡, 查看该序列的估计熔融温度 (tm), 因为候选区域将突出显示。希望获得55-60°c 之间的熔融温度, 同时避免重复和 gcc 的长时间运行。
   3. 按照步骤5.1.1。和5.1.2 选择反向底漆, 位于15-250 基对 3 ' 的正向底漆。虽然引物长度不需要匹配, 但预测的 t m 应尽可能接近前引物的 t m。请务必反转补充序列 (如果在突出显示序列时右键单击 geneious, 则为菜单选项)。
2. 使用"入门设计" 功能, 该功能位于序列窗口的顶部工具栏中。
   1. 点击"入门设计"按钮。在"目标区域" 下插入要放大的区域。
   2. 在 "特性" 选项卡下, 插入所需的大小、熔融温度 (tm) 和% gc (请参阅步骤 5.1.1)。
   3. 单击"确定"生成引物。通过自定义寡核苷酸服务订购引物。
3. 使用控制 dna (编码 t 和 r) 验证引物, 以优化 tm 和扩展时间。使用常规的 taq 和凝胶电泳来查看带的大小, 但也可以按照步骤6中的方法使用 qpcr 进行优化。
   1. 对正向和反向引物进行10倍稀释, 使引物的浓度为10μm。
   2. 使用0.5μl 的 dntp pcr 混合物, 前向底漆0.5 微米、反向底漆0.5μl、taq 聚合酶0.1μl、模板2μl、0.75μl 镁、2.5μl 缓冲液和18.15μl 水, 使每个模板的浓度为 5 ng/2。μl。
   3. 在 pcr 程序中测试不同熔融温度下的引物。通常, 最佳性能是观察到熔体温度略低于预测的 tm 的引物, 但通常不超过60°c。此外, 使用本指南测试最佳扩展时间: 每 1000 bp 1分钟 (因此, 通常 10-30, 具体取决于放大器长度)。
   4. 进行终点凝胶电泳, 以确认引物放大预期的序列。在 200 v 的 2% tae 琼脂糖凝胶上运行25μl 与 5μl 6x 甘油染料混合的 qpcr 产品, 时间为20分钟。

6.剩余序列的 qpcr验证

注: 此步骤需要在步骤5中验证引物和建立 pcr 条件。

1. 运行每个模板一式三份, 并与以下组合;250μl 的 powersybr 绿色主混料, 0.5 微米的正向底漆, 浓度为 10μm, 0.5μl 的反向底漆, 浓度为 10μm, 10.5μl 的水, 和1μl 的模板 dna (浓度为 2ng/μl), 以便每口井包含25μl 的总体积。
2. 运行 qpcr 程序, 根据步骤4中经过验证的温度和延长时间进行通知。我们设计并验证了所有引物, 使其与两级循环兼容, 95°c 为10分钟初始熔体, 然后40个周期的95°c 为 30秒, 60°c 为1分钟。然而, 三阶段 (熔融退火延伸) 程序可能更适合引物, 必要时应进行调整。我们建议至少在 qpcr 中首次使用引物来验证单个 dna 产物的扩增, 从而生成最终的变性曲线。
3. 测量与肌动蛋白 (或任何其他合适的 ' r ' 控制) 有关的 qpcr/sybr 绿色信号, 所有病例计算 2^{-(基因 ct-β-肌动蛋白 ct)}的平均值和标准偏差。
4. (可选)进行端点凝胶电泳, 以确认正确的产品尺寸检测与 qpcr。在这里, 运行25μl 的 qpcr 产品与5μl 的6xl 甘油染料在 2% tae 琼脂糖凝胶在 200 v 的20分钟。

7. 使用 "暂停数据" 的新引用重复此操作。

注: 如果步骤6验证了从 t 中识别的序列, 请在此处结束循环 (图 2a)。但是, 各种考虑因素可能会促使周期的延续, 例如, 如果文件中保留了许多 r 序列, 或者如果在步骤6中没有一个候选 t 序列得到 qpcr 的验证。

1. 获取新的引用。此步骤可实现循环的新迭代, 并可能包括原始基因组数据、原始 rna-seq 数据或其他组装的数据集。用于参考数据的宝贵资源包括国家生物技术信息中心 (https://www.ncbi.nlm.nih.gov/genome) 的基因组数据库, 该数据库储存可通过 ftp (ftp://ftp.ncbi.nlm.nih.gov/genomes/) 访问的组装基因组,和基因表达综合 (https://www.ncbi.nlm.nih.gov/geo/), 其中原始的下一代序列读取存储。基因组项目可以通过其他与项目相关的网站和数据库提供其原始序列数据。

结果

运行 blast 后, 输出文件将具有与数据库匹配的查询序列列表。在 python 减法之后, 将获得一些不匹配的序列, 并通过 qpcr 进行测试。下文将讨论这一结果和下一步工作。

负面的结果.在对参考序列进行对比后, 可以看到两个可能的负面结果。可能没有 blast 结果, 这意味着总序列没有任何类似的序列来引用。这可能是为...

讨论

虽然减法基因组学是强大的, 但它不是一种 cookie 切割器方法, 需要在几个关键步骤中进行自定义, 并仔细选择参考序列和测试样本。如果查询程序集的质量较差, 则筛选步骤可能只隔离程序集项目。因此, 使用特定项目的适当验证协议彻底验证新程序集非常重要。对于 rna-seq, 三一重工网站¹⁸提供了指导方针, 对于 dna, 可以使用 reapr23^这样的工具。使用 blast ?...

材料

Name	Company	Catalog Number	Comments
Accustart II Taq DNA Polymerase	Quanta Bio	95141
Blasic Local Alignment Search Tool (BLAST)			https://github.com/trinityrnaseq/trinityrnaseq/wiki/Transcriptome-Assembly-Quality-Assessment
Bowtie 2			https://www.python.org/download/releases/2.7/
BWA-MEM v. 0.7.12			https://github.com/BenLangmead/bowtie2
Geneious			https://blast.ncbi.nlm.nih.gov/Blast.cgi
PEAR v. 0.9.6			http://www.mybiosoftware.com/reptile-1-1-short-read-error-correction.html
Personal Computer	Biomatters		http://www.geneious.com/
PowerSYBR qPCR mix	ThermoFisher	4367659
Python v. 2.7			https://sco.h-its.org/exelixis/web/software/pear/
Reptile v.1.1			https://alurulab.cc.gatech.edu/reptile
Stratagene Mx3005P	Agilent Technologies	401456
TransDecoder v. 3.0.1			https://sourceforge.net/projects/bio-bwa/files/
Trinity v. 2.4.0			https://github.com/TransDecoder/TransDecoder/wiki