方法文章

RNA下一代测序及生物信息学分析流程在位点特异性水平鉴定表达的LINE-1元件

11.4K 次观看

⸱

DOI：

10.3791/59771

⸱

2019年5月19日

 ,  ,  ,  ,  , 
1Tulane Cancer Center, Tulane University, 2Department of Pathology, Tulane University, 3Department of Structural and Cellular Biology, Tulane University, 4Department of Epidemiology, Tulane University

本文内容

摘要

本文介绍了一种生物信息学方法及相关分析，用于在特定位点水平上鉴定 LINE-1 的表达。

摘要

长散在重复序列-1（Long INterspersed Elements-1，LINEs/L1s）是一类可自我复制并随机插入基因组中的重复元件，可导致基因组不稳定性和诱变效应。在单个位点水平上解析L1元件的表达模式，有助于深入理解这一诱变元件的生物学特性。该自主性转座元件在人类基因组中占比显著，拷贝数超过50万个，但其中99%的拷贝存在截短或缺陷。然而，由于L1元件数量庞大且绝大多数为缺陷型拷贝，因此从作为其他基因组成部分而表达的L1相关序列中鉴定出真实表达的完整L1元件极具挑战性。此外，由于这些元件具有高度重复性，确定具体哪一个L1位点发生表达也十分困难。为克服上述难题，本文提出一种基于RNA-Seq的生物信息学分析方法，用于在位点特异性水平上鉴定L1的表达。简而言之，我们提取细胞质RNA，富集带有多聚腺苷酸尾的转录本，并采用链特异性RNA-Seq分析，将测序读段唯一比对至人类参考基因组中的L1位点。对于每个具有唯一比对读段的L1位点，我们进行可视化人工审阅，以确认其转录起始来源于自身的启动子，并根据各个L1位点的可比对性对映射的转录本读段数量进行校正。该方法已应用于前列腺肿瘤细胞系DU145，以验证本方案检测少数全长L1元件表达的能力。

引言

逆转录转座子是一类重复的DNA元件，可通过RNA中间体以“复制-粘贴”机制在基因组中“跳跃”。其中一类逆转录转座子称为长散在核元件-1（LINEs/L1s），在人类基因组中占比约六分之一，拷贝数超过500,0001。尽管其数量庞大，但大多数拷贝存在缺陷或被截短，估计仅有80至120个L1元件具有活性2。一个完整的L1序列长约6 kb，包含5’和3’非翻译区、一个内部启动子及相关的反义启动子、两个不重叠的开放阅读框（ORFs），以及一个信号序列和polyA尾3,4,5。在人类中，L1由不同进化年龄的亚家族组成，较古老的家族随时间积累的特异性序列突变较多，而最年轻的亚家族为L1HS6,7。L1是唯一具有自主活性的人类逆转录转座子，其开放阅读框编码逆转录酶、内切核酸酶以及具有RNA结合和分子伴侣活性的核糖核蛋白复合物（RNPs），这些成分共同介导逆转录转座过程，并通过一种称为靶位点引物逆转录的机制插入基因组8,9,10,11,12。

已有报道称，L1的逆转座可通过多种机制导致人类种系疾病，包括插入突变、靶位点缺失以及基因组重排13,14,15,16。最近有假说提出，L1可能在肿瘤发生和/或肿瘤进展中发挥作用，因为在多种上皮性癌症中已观察到该致突变元件的表达水平升高及其插入事件增多17,18。据估计，每200次出生中就有一例新的L1插入事件19。因此，深入理解具有活性表达的L1的生物学特性至关重要。然而，由于L1具有重复性特征，且在其他基因的转录本中存在大量缺陷拷贝，使得这一层面的分析极具挑战性。

幸运的是，随着高通量测序技术的出现，研究人员已取得进展，能够在位点特异性水平上解析并鉴定真实表达的L1元件。利用RNA新一代测序技术鉴定表达型L1元件的方法存在不同的理念。目前仅提出两种较为合理的方法用于在位点特异性水平上比对L1转录本。其中一种方法仅关注那些能够通读L1多聚腺苷酸化信号并延伸至侧翼序列的潜在转录事件20。我们的方法则利用L1元件之间微小的序列差异，仅比对那些唯一映射到单一基因组位点的RNA-Seq读段21。这两种方法在转录本水平的定量方面均存在一定局限性。通过为每个L1位点引入“唯一可比对性”校正因子21，或采用更复杂的算法重新分配那些无法唯一比对至特定基因组位点的多重比对读段22，可能有助于提高定量准确性。本文将逐步详细介绍RNA提取、新一代测序及生物信息学分析流程，以在位点特异性水平上鉴定表达的L1元件。我们的方法最大限度地利用了对功能性L1元件生物学特性的已有认知，包括：功能性L1元件必须由L1启动子驱动产生，转录起始于L1元件的起始位置，其转录本应在细胞质中被翻译，且其转录本应与基因组序列共线性排列。简而言之，我们收集新鲜的细胞质RNA，富集带有poly(A)尾的转录本，并采用链特异性的RNA-Seq分析方法，将测序读段唯一地比对至人类参考基因组中的L1位点。随后，这些比对结果仍需进行大量人工审阅，以确认转录本读段是否真正起源于L1启动子，方可将某一位点判定为真实表达的L1元件。我们以DU145前列腺肿瘤细胞系样本为例，展示该方法如何从大量失活的L1拷贝中识别出少数正在活跃转录的L1成员。

方案

1. 胞质RNA提取

  1. 通过以下方法获取细胞。
    1. 从融合度为2.75%–100%的T-75培养瓶中收集活细胞。
      1. 用5 mL冷PBS洗涤培养瓶2次，最后一次洗涤时刮下细胞并转移至15 mL锥形管中。在1,000 × g、4 °C条件下离心2分钟，小心移除并弃去上清液（材料表）。
    2. 从组织样本中收集细胞。
      1. 在组织解剖后1小时内进行细胞质RNA提取的组织制备，并始终将组织置于冰上。如需长期保存，按照制造商说明书使用RNA抑制剂溶液保存组织，解剖后最多可保存72小时（材料表）。
      2. 将10 µm3的组织样本切碎，在无菌Dounce匀浆器中用5 mL冷PBS匀浆新鲜样本，转移至15 mL锥形管中，在1,000 × g、4 °C条件下离心2分钟，小心移除并弃去上清液（材料表）。
  2. 向细胞沉淀中加入2 mL裂解缓冲液，混匀并在冰上孵育5分钟。
    1. 使用150 mM NaCl、50 mM HEPES（pH 7.4）和25 μg/mL皂苷（digitonin）新鲜配制裂解缓冲液（材料表）。
    2. 由于不同细胞类型所需的裂解缓冲液中皂苷最低浓度可能不同，需通过显微镜确认经裂解缓冲液处理的细胞已失去质膜但核膜保持完整。
    3. 使用前立即加入1,000 U/mL RNase抑制剂（材料表）。
  3. 在1,000 × g、4 °C条件下离心1分钟，收集上清液。
  4. 将上清液加入预冷的7.5 mL Trizol和1.5 mL氯仿中。所有涉及氯仿的操作均须在洁净的化学通风橱内进行（材料表）。
  5. 在3,220 × g、4 °C条件下离心35分钟。
  6. 将水相（上层）转移至新的预冷15 mL离心管中。
  7. 加入4.5 mL氯仿，涡旋振荡。
  8. 在3,220 × g、4 °C条件下离心10分钟。
  9. 将水相转移至新的预冷离心管中。
  10. 加入4.5 mL异丙醇，充分混匀，并在-80 °C下孵育过夜（材料表）。
  11. 在3,220 × g、4 °C条件下离心45分钟。
  12. 弃去异丙醇，加入15 mL 100%乙醇（材料表）。
  13. 在3,220 × g条件下离心10分钟。
  14. 弃去乙醇，倒置离心管沥干并干燥约1小时。
    1. 使用无菌棉签吸除残留的乙醇（材料表）。
  15. 根据沉淀大小，用100至200 μL无RNase的水重悬样本（材料表）。
  16. 根据制造商说明书，采用电泳技术对样本进行分馏，以测定样本的质量和浓度23（材料表）。
    1. 若RIN > 8，则样本符合RNA-Seq分析要求24。

2. 下一代测序

  1. 将细胞质RNA样本提交至下一代测序平台进行测序，目标是生成至少5000万条成对末端、100 bp的读长。
  2. 选择带有聚腺苷酸化修饰的RNA并进行链特异性测序。

3. 创建注释（如果已有现有注释，则此步骤可选）

  1. 创建全长L1注释或下载全长L1注释（补充文件1a-b）。
    1. 使用UCSC基因组浏览器的表格浏览器工具（https://genome.ucsc.edu/cgi-bin/hgTables）下载LINE-1元件的Repeat Masker注释。指定哺乳动物分支、人类基因组、hg19组装版本（或使用更新的hg38版本），并在“Class Name”下筛选“LINE1”。将其下载为.gtf文件，并命名为FL-L1-BLAST.gtf。
    2. 在人类基因组中对包含启动子区域的L1.3全长L1元件的前300 bp进行本地BLAST搜索，并向下游延伸6,000 bp，以确定L1坐标的末端，然后将该信息添加至注释文件中。保存为gtf文件，并命名为FL-L1-RM.gtf。
    3. 使用bedtools将RepeatMasker注释与基于启动子的L1注释取交集，并命名为FL-L1-BLAST_RM.txt （软件包）。
      1. 在Linux终端中使用以下命令：bedtools intersect -a FL-L1-BLAST.gtf -b FL-L1-RM.gtf > FL-L1-BLAST_RM.txt。
    4. 根据正链和负链对交集后的FL-L1注释进行分离。
      1. 将FL-L1-BLAST_RM.txt复制到电子表格软件中，按“minus”和“plus”链排序，然后按染色体位置排序。
      2. 创建两个新的电子表格文档，一个包含负链上的全长L1交集坐标，另一个包含正链上的全长L1交集坐标，并分别保存为FL-L1-BLAST_RM_minus.xls和FL-L1-BLAST_RM_plus.xls。
      3. 将这两个新文档另存为.txt文件。
    5. 使用mac2unix程序将.txt文件转换为正确的注释文件格式（软件包）。
      1. 在终端中使用以下命令：mac2unix.sh FL-L1-BLAST_RM_minus.gff。
      2. 在终端中使用以下命令：mac2unix.sh FL-L1-BLAST_RM_plus.gff。
      3. 将新文件保存为.gff扩展名。
    6. 或者，使用AWK命令筛选与+链和–链相关的行。
      1. 使用以下命令获取+链：awk ‘/+/’ FL-L1_BLAST_RM.gtf > FL-L1_BLAST_RM_plus.gtf。
      2. 使用以下命令行获取-链：awk ‘/-/’ FL-L1_BLAST_RM.gtf > FL-L1_BLAST_RM_minus.gtf。

4. 读段比对流程以鉴定表达的 L1

选项描述
–p此参数指定计算机在执行比对时应使用的线程数。更大的计算机内存可支持更多线程，具体数值应通过实验确定。
–m 1此参数指示程序仅接受在基因组中具有唯一最佳匹配（优于其他任何基因组匹配）的读段。
–y此为“tryhard”开关，使比对程序搜索所有可能的匹配位置，而不因达到固定数量的匹配而提前终止。
–v 3此参数限制程序仅对与基因组比对错配数不超过3个的读段使用内存。
–X 600此参数仅允许相互距离在600个碱基以内的成对读段被比对。这确保读段对在基因组中具有共线性，并排除涉及加工后RNA分子的结构。
–chunkmbs 8184此命令为每个与L1相关的读段可能产生的大量比对结果分配额外内存。

表1：Bowtie的命令行选项。

  1. 使用 Bowtie 将配对末端测序的 fastq 文件与感兴趣的 RNA-Seq 样本进行比对。
    注意：必须使用 Bowtie1 而非 Bowtie2，因为唯一比对所需的参数仅在 Bowtie 的此版本中提供（软件包）。相较于 STAR 等可识别剪接的比对工具，Bowtie 更适用于评估与 L1 生物学及表达更相关的共线性、连续性读段。
    1. 在 Linux 终端中使用以下命令行：bowtie -p 10 -m 1 -S -y -v 3 -X 600 --chunkmbs 8184 hg_X_Y_M_index -1 hg_sample_1.fq -2 hg_sample_2.fq | samtools view -hbuS - | samtools sort – hg_sample_sorted.bam。Bowtie 命令行选项的说明见 表 1。
  2. 使用 samtools（软件包）及以下 Linux 命令对输出的 bam 文件进行链特异性分离。请注意，若未采用标准的下一代测序流程，实际的标志值可能有所不同。
    1. 使用以下命令行选择正义链（top strand）：samtools view -h hg_sample_sorted.bam | awk 'substr($0,1,1) == "@" || $2 == 83 || $2 == 163 {print}' | samtools view -bS - > hg_sample_sorted_topstrand.bam。
    2. 使用以下命令行选择反义链（bottom strand）：samtools view -h hg_sample_sorted.bam | awk 'substr($0,1,1) == "@" || $2 == 99 || $2 == 147 {print}' | samtools view -bS - > hg_sample_sorted_bottomstrand.bam。
  3. 使用 bedtools（软件包）基于 L1 位点的注释生成读段计数。
    1. 使用以下命令行为正义链上的 L1 生成正义方向的读段计数：bedtools coverage -abam FL-L1-BLAST_RM_plus.gff -b hg_sample_sorted_topstrand.bam > hg_sample_sorted_bowtie_tryhard_plus_top.txt。
    2. 使用以下命令行为反义链上的 L1 生成正义方向的读段计数：bedtools coverage -abam FL-L1-BLAST_RM_minus.gff -b hg_sample_sorted_bottomstrand.bam > hg_sample_sorted_bowtie_tryhard_minus_bottom.txt。
  4. 对步骤 5.1.1 中的 bam 文件进行索引，以便在 Integrative Genomics Viewer (IGV)25（软件包）中查看。
    1. 使用以下命令行：samtools index hg_sample_sorted.bam
  5. 若需采用批处理模式以同时处理更多 RNA-Seq 样本，请使用超级计算机脚本完成步骤 4.1（脚本名为 human_bowtie.sh），用于完成步骤 4.2–4.3 的脚本为 human_L1_pipeline.sh，用于完成步骤 4.4 的脚本为 bam_index.sh。这些脚本可在 补充文件 2 中找到，并附有用于运行脚本的相应超级计算机命令。

5. 人工审校

  1. 创建一个电子表格，用于记录比对到每个注释 L1 位点的读段数。
    1. 复制步骤 4.3.2 中生成的 hg_sample_sorted_bowtie_tryhard_minus_bottom.txt 文件，并将该工作表标签命名为“minus-bottom”。
      1. 根据 J 列中读段数从高到低对所有列进行排序。
    2. 复制步骤 4.3.1 中生成的 hg_sample_sorted_bowtie_tryhard_plus_top.txt 文件，并在另一个工作表中将其标签命名为“top-plus”。
      1. 根据 J 列中读段数从高到低对所有列进行排序。
    3. 创建第三个标签为“combined”的工作表，并添加来自“minus-bottom”和“top-plus”工作表中读段数不少于十个的所有位点。
      1. 根据 J 列中读段数从高到低对所有列进行排序。
    4. 将以下文件加载到 IGV25（软件包）中：1）感兴趣的参考基因组，用于可视化注释基因；2）FL-L1-BLAST_RM.gff，用于可视化 L1 注释；3）hg_sample_sorted.bam，用于可视化来自目标样本的比对转录本；4）hg_genomicDNA_sorted.bam，用于评估基因组区域的可比对性。
    5. 移除与每个 bam 文件关联的覆盖度和剪接连接行。
    6. 压缩 hg_sample_sorted.bam 和 hg_genomicDNA_sorted.bam，以便所有 IGV 轨道能够显示在同一屏幕内。
  2. 手动审校。
    1. 利用“combined”工作表中列出的位点坐标，在 IGV25（软件包）中查看所识别的位点。
    2. 若在 L1 方向上游 5 kb 范围内无任何读段，则可判定该位点为由其自身启动子真实表达。
      1. 将该行标记为绿色，并注明其为真实表达 L1 的原因。
        注：若 L1 上游区域不可比对，则为本规则的例外情况。若属此类情况，将该行标记为红色，并注明无法评估 L1 启动子上游区域的表达，因此无法可靠判断该 L1 的表达状态。
    3. 若在上游 5 kb 范围内存在读段，则判定该位点并非由其自身启动子真实表达。
      1. 将该行标记为红色，并注明其非真实表达 L1 的原因。
      2. 若某位点位于一个同向表达基因的内含子内且 L1 上游存在读段，或位于一个同向表达基因的下游且 L1 上游存在读段，或表现为未注释的表达模式且 L1 上游存在读段，则将其审校为假阳性。
        注：若存在一种例外情况，即在 L1 启动子起始位点正上方仅有少量读段重叠，但略位于 L1 上游，且该 L1 上游无其他读段，则可认为该 L1 为真实表达。将该行标记为绿色，并注明其为真实表达 L1 的原因。
    4. 若比对到该位点的读段分布模式与该特定 L1 区域的可比对性模式不一致，则将其审校为很可能为假阳性。
      注：例如，若某 L1 区域高度可比对，但读段仅在 L1 内某一紧凑区域堆积，则其更可能并非源自该 L1 自身启动子的表达，而更可能来自未注释的来源（如外显子或 LTR）。此类情况下，将该位点标记为橙色，并注明其可疑原因。通过在 UCSC 基因组浏览器中检查 L1 位置，验证可疑堆积信号的来源。
    5. 若某位点位于基因组中零星表达的未注释区域环境中，则判定其非真实表达。
      注：例如，L1 上游 10 kb 处可能存在表达读段，且每隔约 10 kb 均有比对读段，其中部分读段与 L1 对齐。此类 L1 更可能并非由其自身启动子驱动表达，而更可能是由于基因组表达的未注释模式导致读段比对至该位点。此类情况下，将该位点标记为橙色，并注明其可疑原因。

6. 评估参考基因组中可比对性的读段比对策略（若已有现成的基因组DNA比对数据集，则此步骤可选）

  1. 下载全基因组DNA序列文件并转换为.fq文件
    1. 访问NCBI网站：https://www.ncbi.nlm.nih.gov/sra
    2. 在搜索框中输入 WGS HeLa paired end。
    3. 在 Results by taxon 下选择 Homo sapiens。
    4. 选择一个双端测序（paired end）样本，且读长为100 bp或以上，例如以下样本：https://www.ncbi.nlm.nih.gov/sra/ERX457838[accn]
    5. 通过点击 Run，然后选择 Metadata 来确认读长，示例如下：https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=ERR492384
    6. 在Linux终端中输入以下命令以下载全基因组DNA序列数据：sratoolkit.2.9.2-mac64/bin/prefetch -X 100G ERR492384
      注：SRA工具包的prefetch功能会从NCBI网站下载编号为“ERR492384”的数据（见Software Packages）。参数“100G”将下载数据量限制为100吉字节。
    7. 在Linux终端中输入以下命令：fastq-dump --split-files ERR492384
      注：该命令将下载的基因组DNA数据集拆分为两个fastq文件。
  2. 使用Bowtie进行比对
    1. 在Linux中使用以下命令进行比对：bowtie -p 10 -m 1 -S -y -v 3 -X 600 --chunkmbs 8184 hg_X_Y_M_index -1 hg_genomicDNA_1.fq -2 hg_genomicDNA_2.fq | samtools view -hbuS - | samtools sort – hg_genomicDNA_sorted.bam。
      1. 参见步骤4.1以了解Bowtie比对中所用参数的含义（见Software Packages）。
      2. 基因组比对后的bam文件可应作者请求获取，用于评估可比对性（mappability）。
  3. 使用samtools对步骤4.2.1生成的bam文件进行索引，以便在IGV25（见Software Packages）中查看，从而进一步支持人工审校。
    1. 在Linux中使用以下命令：samtools index hg_genomicDNA_sorted.bam
  4. 评估每个L1位点的可比对性（mappability）
    1. 使用bedtools程序、FL-L1注释文件以及比对后的基因组序列数据，确定唯一比对到L1位点的读段数量（见Software Packages）。
      1. 在Linux中使用以下命令：bedtools coverage -abam FL-L1-BLAST_RM.gtf –b hg_genomicDNA_sorted.bam > L1_Mappability_hg_genomicDNA.txt。
    2. 当一个L1位点上有400条唯一比对读段时，将其定义为具有完全覆盖的可比对性。
    3. 计算将每个L1位点的基因组DNA比对读段数标准化至400所需的缩放因子。
    4. 为根据各个L1位点的可比对性对表达量进行标准化，将步骤6.4.3中计算出的缩放因子乘以第4–5节中确定的真实表达L1位点的RNA转录本比对读段数。

结果

上述步骤及图1中的图示说明已应用于人前列腺肿瘤细胞系DU145。RNA样本经细胞质提取后，采用poly-A选择性、链特异性、双端测序方法进行高通量测序。使用Bowtie软件比对双端测序数据，仅保留唯一比对结果，即双端序列读段比对到某一基因组位置的匹配度优于其他任何基因组位置。将DU145测序数据比对至人类参考基因组，生成bam文件，该文件可应作者要求获取。利用bedtools从DU145链分离的bam文件中提取比对至全长L1序列的读段数量信息，并将这些读段按数量从大到小排序，随后通过IGV软件人工审阅每个L1位点周围的基因组环境，以确认其真实性（补充表1）。若某样本被判定为真实表达，则在最右侧列以绿色标注，并说明其被接受的理由。符合方法部分所述标准、被判定为真实表达的L1位点示例如图2a-b所示。若某样本被判定为非真实表达，则在最右侧列以红色标注，并注明拒绝理由。根据方法部分所述标准，因表达来源于非自身启动子而被拒绝的L1位点示例详见图2c-e。

本研究仅针对具有完整启动子区域的全长L1序列。如果不作此区分，将会引入大量来源于截短型L1序列的转录噪声。在DU145细胞中，截短型L1的示例见图3a-b，这些序列通过唯一比对的RNA-Seq读段被识别出来。然而，在IGV中可以明显看出，这些转录本并非起始于截短型L1本身，而是由于L1序列被包含在某个基因内，或位于一个已表达基因的下游所致。

总体而言，在DU145细胞中，经过人工审校后被排除为真实表达的L1位点及其测序读数的比例约为50%（补充表2），这表明若不进行人工审校，大量比对到L1的转录本读数将被误判为假阳性。具体而言，在DU145细胞中共有114个全长L1位点具有唯一比对到正义链的测序读数，总计3,152条读数；但经过人工审校后，仅有60个位点被确认为由自身启动子驱动表达，共包含1,879条读数（补充表1）。即使已通过选择细胞质mRNA来减少与L1生物学无关的表达，这一情况仍然存在。需要注意的是，在DU145中比对到转录本最多的L1位点因不属于真实表达的L1而被排除（图4）。总体来看，经人工审校后被接受或排除为真实表达的L1位点，其比对到的转录本数量在两者之间分布相似（图4）。

经过人工审校后，能够唯一比对到DU145细胞中真实表达的特定L1位点的读段数范围为175条至人为设定的最低阈值10条（图5）。这种通过唯一比对转录本读段来鉴定L1位点的方法限制了对表达水平进行准确定量的能力。为校正这一偏差，我们为每个位点基于其可比对性（mappability）建立了校正因子。为构建该校正因子，首先使用bedtools从HeLa基因组bam文件中提取比对至所有全长L1位点的唯一比对读段数，并将这些位点按读段数从高到低作图（补充图1）。我们人为设定，HeLa基因组测序样本中具有400条读段的L1位点具有完全的可比对性。将各L1位点在HeLa基因组测序样本中所能比对到的读段数相对于400进行归一化处理，所得归一化数值再乘以在DU145细胞中比对到各真实表达L1位点的读段数（补充表2）。如预期所示，可比对性校正分数较高的L1元件主要来自较年轻的亚家族，如L1PA2（补充表2）。在根据各L1位点的可比对性分数对读段数进行校正后，大多数位点的表达定量值有所提高（图6）。经可比对性校正后，唯一比对到DU145细胞中真实表达的特定L1位点的读段数范围为612至4条，且高表达至低表达位点的排序发生了变化（图6）。

图解RNA测序流程、基因组比对及在IGV中逐步进行数据校正的流程图
图1：实验流程示意图。
图示描述了在人类样本中鉴定表达型L1元件的各个步骤。请注意，如果已有适当的文件，步骤1和步骤2无需重复进行。这些适当文件可从补充文件1a-b和补充文件2中下载获得。红色框标示的步骤表示使用bedtools coverage程序统计比对到L1元件且方向一致的测序读段数量。这些具有同向比对读段的基因组位点即为需要人工校正的L1位点。请点击此处查看该图的放大版本。

带有FL-L1注释的RNA-Seq样本；基因组比对示意图；hg19参考基因组；可比对性评估。
图2：DU145细胞中 curated L1位点的示例。
在IGV中加载了参考基因组、与参考基因组版本匹配的全长L1 gff注释文件（补充文件1）、DU145的bam文件，以及最后用于评估可比对性的HeLa基因组bam文件，所有数据均可应作者要求获取。图中添加了箭头以帮助可视化注释L1的方向。红色箭头和测序读段表示序列方向从右向左；蓝色箭头和读段表示序列方向从左向右。a) 在IGV中，该L1位点似乎由其自身的启动子驱动表达，因为在L1上游超过5 kb范围内无正义链方向的测序读段。该L1可比对性较低，不位于任何基因内，并且存在预期的反义链启动子活性证据26。b) 在IGV中，该L1位点似乎由其自身的启动子驱动表达，因为在L1上游超过5 kb范围内无正义链方向的测序读段。该L1可比对性较低，且位于一个方向相反的基因内部。c) 在IGV中，该L1位点被排除为表达的L1，因为在5 kb范围内存在同方向的上游读段。该L1位于一个同方向基因内部，因此转录本读段很可能来源于该表达基因的启动子。d) 在IGV中，该L1位点被排除为表达的L1，因为在5 kb范围内存在同方向的上游读段。该L1位于一个同方向的高表达基因下游，因此转录本读段很可能来源于该表达基因的启动子，并延伸超过了正常的基因终止子。e) 在IGV中，该L1位点被排除为表达的L1，因为在5 kb范围内存在同方向的上游读段。该L1不在参考基因注释中任何基因的内部或附近，因此这些在L1元件内部及上游的转录本提示可能存在未被注释的启动子。请点击此处查看该图的放大版本。

基因组定位图，hg19参考基因组，截短L1注释，RNA-Seq分析，可比对性。
图3：背景噪音同样来源于截短型L1元件。
我们的L1注释未包含截短型L1，因为它们是背景噪音的主要来源。图中添加了箭头以帮助可视化注释L1的方向。蓝色的箭头和测序读段表示在序列中从左至右的方向。a) 展示了一个属于L1MB5亚家族的截短L1实例，长度为2706 bp。在IGV中可以明显看出，这些读段来源于一个表达基因的下游延伸区域。b) 展示了另一个截短L1的例子。该L1属于L1PA11，长度为4767 bp。在IGV中可以明显看出，唯一比对到该L1的读段来源于其所在基因的已表达外显子。请点击此处查看该图的放大版本。

L1表达条形图；转录本比对读段数与L1位点；基因表达分析；DU145细胞。
图4：在人基因组中比对到所有全长完整L1位点的唯一转录本读段，来源于DU145前列腺肿瘤细胞系。
黑色表示经人工审校后确认为真实表达的特定L1位点，红色表示经人工审校后排除为真实表达的特定L1位点。灰色表示每个位点比对读段数少于10条的L1位点。由于这些位点所占转录本读段比例极小，未进行人工审校。横轴刻度标记每100个全长完整L1位点。另有约4,500个位点因无任何比对读段而未在图中显示。请点击此处查看该图的放大版本。

转录本读段的柱状图；DU145细胞系；基因表达分析；L1元件研究
图5：在DU145前列腺肿瘤细胞系中唯一比对到真实表达的全长完整L1元件的转录本读段。
图中显示的是经过人工审校后，在DU145细胞中比对到特定基因座的转录本读段数量。请点击此处查看该图的放大版本。

显示 DU145 数据集分析中原始与可比对性校正后转录本读段的柱状图。
图 6：经可比对性校正后映射到真实表达 L1 的读段数。
图中显示的是经位点特异性可比对性评分校正后，映射到人工审编的 DU145 细胞中 L1 位点的转录本读段数量。请点击此处查看该图的放大版本。

补充文件 1：根据方向性标注的全长完整人类 L1 序列。a) FL-L1-BLAST_RM_minus.gff。 b) FL-L1-BLAST_RM_plus.gff。 请点击此处下载该文件。

补充文件2：用于自动化第4节所述生物信息学流程的超级计算机脚本。 请点击此处下载该文件。

补充图1：用于确定L1可比对性的基因组DNA样本。
图中显示了来自HeLa细胞系样本的基因组转录本reads中，唯一比对到基因组中全部5,000个全长L1位点的reads数量。当有400条reads比对到某个L1时，即判定该L1具有完全覆盖的可比对性。请点击此处下载该图。

补充表 1：DU145 中 L1 的人工审编。 请点击此处下载该表格。

补充表 2：经可比对性校正的 DU145 细胞中筛选的 L1 序列。 请点击此处下载该表格。

讨论

L1活性已被证明可导致遗传损伤和基因组不稳定，从而促进疾病的发生27,28,29。在约5000个全长L1拷贝中，仅有几十个进化上较年轻的L1元件贡献了大部分的逆转座活性2。然而，有证据表明，即使是一些较古老且不具备逆转座能力的L1元件，仍可能产生导致DNA损伤的蛋白质30。要全面理解L1在基因组不稳定性和疾病中的作用，必须明确L1在单个基因座水平上的表达情况。然而，大量与L1相关的序列作为背景被整合到其他与L1逆转座无关的RNA中，这为识别真实的L1表达带来了重大挑战。此外，由于L1序列具有重复性，许多短读长测序序列无法唯一比对到单个特定基因座，这也阻碍了对单个L1基因座表达模式的识别与研究。为克服这些挑战，我们开发了上述利用RNA-Seq数据鉴定单个L1基因座表达的方法。

我们的方法通过多个步骤，过滤掉与L1逆转座无关的L1序列所产生的高水平（超过99%）转录噪声。第一步是制备细胞质RNA。通过选择细胞质RNA，可显著减少在细胞核内表达的含内含子mRNA中出现的L1相关读段。在测序文库构建过程中，另一个用于降低与L1无关的转录噪声的步骤是选择带有polyA尾的转录本，从而去除存在于非mRNA物种中的L1相关转录噪声。此外，采用链特异性测序以识别并剔除反义链L1相关转录本。在鉴定映射到L1的RNA-Seq转录本数量时，使用包含功能性启动子区域的全长L1注释信息，也可消除原本来源于截短L1的背景噪声。最后，消除与L1逆转座无关的L1序列转录噪声的关键步骤是对鉴定出具有映射RNA-Seq转录本的全长L1进行人工审校。该人工审校过程包括在基因组周围环境背景下可视化每一个经生物信息学鉴定为表达的L1位点，以确认其表达确实起源于L1自身的启动子。该方法已应用于前列腺肿瘤细胞系DU145。即使采取了所有降低背景噪声的制备步骤，在DU145中经生物信息学鉴定出的L1位点仍有约50%被判定为来源于其他转录来源的L1背景噪声（图4），这凸显了获得可靠结果所需的高度严谨性。尽管这种结合人工审校的方法工作量较大，但在本分析流程的建立过程中，对于评估和理解全长L1周围的基因组环境而言是必要的。下一步工作包括通过自动化部分审校规则来减少所需的人工审校工作量；然而，由于基因组表达特性尚未完全明确、参考基因组中存在未注释的表达来源、低可比对区域，以及参考基因组构建过程中存在的复杂因素，目前尚无法完全实现L1审校的自动化。

通过测序鉴定单个L1位点表达的第二个挑战与重复性L1转录本的比对有关。在此比对策略中，要求转录本必须能够唯一且共线性地比对到参考基因组，才能被成功定位。通过筛选那些一致性比对的成对末端序列，可比对到参考基因组中L1位点的转录本数量得以增加。这种唯一比对策略能够确保比对到单个L1位点的读段具有较高的可信度，但可能会低估每个被确认为真实表达的重复性L1的表达量。为了大致校正这种低估，我们开发并应用了一种基于各L1位点“可比对性”（mappability）的评分，将其应用于唯一比对的转录本读段数量（图6）。需要注意的是，理想情况下，应根据匹配的全基因组测序（WGS）样本，在全长L1上对完整覆盖度的读段进行可比对性评分。在此，我们使用HeLa细胞的WGS数据来确定每个L1位点的可比对性评分，从而对DU145前列腺肿瘤细胞系中比对到L1位点的读段数进行上调或下调校正。该可比对性计算是一种粗略的校正方法，但所选定的“完全覆盖可比对性”阈值为400个读段，这一数值是在考虑肿瘤细胞系动态特性的基础上确定的。从补充图1中可以观察到，部分L1位点在HeLa WGS中比对到的读段数极高，这些可能来源于HeLa细胞中存在但未包含在参考基因组中的染色体重复序列，因此这些位点未被选为完全可比对性覆盖的代表。相反，根据补充图1，100%读段覆盖的平均水平出现在约400个读段处，因此我们假设该平均值同样适用于DU145前列腺肿瘤细胞系。

利用RNA-Seq技术产生的100-200 bp读长进行比对的这一策略，倾向于优先选择参考基因组中进化上更古老的L1元件，因为较古老的L1随时间累积了独特的突变，使其更易于比对。因此，该方法在识别最年轻的L1以及非参考序列的多态性L1时灵敏度有限。为了鉴定最年轻的L1，我们建议采用5’ RACE方法富集L1转录本，并结合使用PacBio等长读长测序技术21。这种方法可实现更特异的比对，从而可靠地鉴定出表达的年轻L1。联合使用RNA-Seq和PacBio技术，可以获得更全面的真实表达L1列表。要鉴定真实表达的多态性L1，下一步的关键步骤包括构建多态性序列并将其插入参考基因组中。

尽管研究重复序列存在诸多生物学和技术上的挑战，但通过上述严格的实验流程，利用RNA测序技术去除与逆转座无关的L1序列的转录噪声，我们能够逐步筛选出高水平的转录背景噪声，从而在单个基因座水平上可靠且严格地鉴定L1的表达模式及其表达量。

披露

作者无任何利益冲突需要披露。

致谢

我们感谢董岩博士提供DU145前列腺肿瘤细胞。我们感谢Nathan Ungerleider博士在编写超级计算机脚本方面提供的指导和建议。本研究的部分工作由美国国立卫生研究院（NIH）资助，项目编号分别为PD的R01 GM121812、VPB的R01 AG057597以及TK的5TL1TR001418。我们还感谢Cancer Crusaders和杜兰大学癌症中心生物信息学核心团队提供的支持。

材料

本文使用的材料清单
姓名公司目录编号评论
1 M HEPESAffymetrixAAJ16924AE
5 M NaClInvitrogenAM9760G
Agilent bioanalyzer 2100Agilent technologies
Agilent RNA 6000 Nano KitAgilent technologies5067-1511
bedtools.26.0https://bedtools.readthedocs.io/en/latest/content/installation.html
bowtie-0.12.8https://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12.8/
细胞刮刀Olympus plastics25-270
氯仿FisherC298-500
皂苷Research Products International Corp50-488-644
乙醇FisherA4094
Gibco（磷酸盐缓冲液）Invitrogen10-010-049
匀浆器Thomas ScientificBBI-8541906
IGV 2.4https://software.broadinstitute.org/software/igv/download
异丙醇FisherA416-500
mac2unixhttps://sourceforge.net/projects/cs-cmdtools/files/mac2unix/
棉签Fisher23-400-122
RNA later 溶液InvitrogenAM7022
RNaseZap RNase 去污染溶液InvitrogenAM9780
samtools-1.3https://sourceforge.net/projects/samtools/files/
sratoolkit.2.9.2https://github.com/ncbi/sra-tools/wiki/Downloads
SUPERase·In RNase 抑制剂InvitrogenAM2694
TrizolInvitrogen15-596-018
水（无 DNASE、RNASE）FisherBP2484100

参考文献

  1. International Human Genome Sequencing. Initial sequencing and analysis of the human genome. Nature. 409, 860(2001).
  2. Brouha, B., et al. Hot L1s account for the bulk of retrotransposition in the human population. Proceedings of the National Academy of Sciences of the United States of America. 100 (9), 5280-5285 (2003).
  3. Dombroski, B. A., Mathias, S. L., Nanthakumar, E., Scott, A. F., Kazazian, H. H. Isolation of an active human transposable element. Science. 254 (5039), 1805(1991).
  4. Swergold, G. D. Identification, characterization, and cell specificity of a human LINE-1 promoter. Molecular and Cellular Biology. 10 (12), 6718-6729 (1990).
  5. Speek, M. Antisense promoter of human L1 retrotransposon drives transcription of adjacent cellular genes. Molecular and Cellular Biology. 21 (6), 1973-1985 (2001).
  6. Deininger, L., Batzer, M. A., Hutchison, C. A., Edgell, M. H. Master genes in mammalian repetitive DNA amplification. Trends in Genetics. 8 (9), 307-311 (1992).
  7. Boissinot, S., Chevret, P., Furano, A. L1 (LINE-1) Retrotransposon Evolution and Amplification in Recent Human History. Molecular Biology and Evolution. 17 (6), 915-918 (2000).
  8. Khazina, E., Weichenrieder, O. Non-LTR retrotransposons encode noncanonical RRM domains in their first open reading frame. Proceedings of the National Academy of Sciences of the United States of America. 106 (3), 731-736 (2009).
  9. Martin, S. L., Bushman, F. D. Nucleic acid chaperone activity of the ORF1 protein from the mouse LINE-1 retrotransposon. Molecular and Cellular Biology. 21 (2), 467-475 (2001).
  10. Feng, Q., Moran, M. H., Kazazian, H. H., Boeke, J. D. Human L1 Retrotransposon Encodes a Conserved Endonuclease Required for Retrotransposition. Cell. 87 (5), 905-916 (1996).
  11. Mathias, S. L., Scott, A. F., Kazazian, H. H., Boeke, J. D., Gabriel, A. Reverse transcriptase encoded by a human transposable element. Science. 254 (5039), 1808(1991).
  12. Luan, D. D., Korman, M. H., Jakubczak, J. L., Eickbush, T. H. Reverse transcription of R2Bm RNA is primed by a nick at the chromosomal target site: A mechanism for non-LTR retrotransposition. Cell. 72 (4), 595-605 (1993).
  13. van den Hurk, J. A. J. M., et al. Novel types of mutation in the choroideremia (CHM) gene: a full-length L1 insertion and an intronic mutation activating a cryptic exon. Human Genetics. 113 (3), 268-275 (2003).
  14. Miné, M., et al. A large genomic deletion in the PDHX gene caused by the retrotranspositional insertion of a full-length LINE-1 element. Human Mutation. 28 (2), 137-142 (2007).
  15. Solyom, S., et al. Pathogenic orphan transduction created by a nonreference LINE-1 retrotransposon. Human Mutation. 33 (2), 369-371 (2012).
  16. Hancks, D. C., Kazazian, H. H. Roles for retrotransposon insertions in human disease. Mobile DNA. Mobile DNA. 7, 9-9 (2016).
  17. Tubio, J. M. C., et al. Mobile DNA in cancer. Extensive transduction of nonrepetitive DNA mediated by L1 retrotransposition in cancer genomes. Science. 345 (6196), 1251343-1251343 (2014).
  18. Ewing, A. D., et al. Widespread somatic L1 retrotransposition occurs early during gastrointestinal cancer evolution. Genome Research. 25 (10), 1536-1545 (2015).
  19. Beck, C. R., Garcia-Perez, J. L., Badge, R. M., Moran, J. V. LINE-1 elements in structural variation and disease. Annual Review of Genomics and Human Genetics. 12, 187-215 (2011).
  20. Philippe, C., et al. Activation of individual L1 retrotransposon instances is restricted to cell-type dependent permissive loci. eLife. 5, e13926(2016).
  21. Deininger, P., et al. A comprehensive approach to expression of L1 loci. Nucleic Acids Research. 45 (5), e31-e31 (2017).
  22. Jin, Y., Tam, O. H., Paniagua, E., Hammell, M. TEtranscripts: a package for including transposable elements in differential expression analysis of RNA-seq datasets. Bioinformatics. 31 (22), 3593-3599 (2015).
  23. Agilent RNA 6000 Nano Kit Guide. , Agilent. (2017).
  24. Mueller, O. L., Schroeder, A. RNA Integrity Number (RIN) –Standardization of RNA Quality Control. , Agilent Technologies. (2016).
  25. Robinson, J. T., et al. Integrative genomics viewer. Nature Biotechnology. 29, 24(2011).
  26. Speek, M. Antisense promoter of human L1 retrotransposon drives transcription of adjacent cellular genes. Molecular Cellular Biology. 21 (6), 1973-1985 (2001).
  27. Belancio, V. P., Deininger, L., Roy-Engel, A. M. LINE dancing in the human genome: transposable elements and disease. Genome Medicine. 1 (10), 97-97 (2009).
  28. Iskow, R. C., et al. Natural Mutagenesis of Human Genomes by Endogenous Retrotransposons. Cell. 141 (7), 1253-1261 (2010).
  29. Scott, E. C., et al. A hot L1 retrotransposon evades somatic repression and initiates human colorectal cancer. Genome Research. 26 (6), 745-755 (2016).
  30. Kines, K. J., Sokolowski, M., deHaro, D. L., Christian, C. M., Belancio, V. P. Potential for genomic instability associated with retrotranspositionally-incompetent L1 loci. Nucleic Acids Research. 42 (16), 10488-10502 (2014).

重印与许可

申请许可以重复使用本 JoVE 文章的文本或图表

申请许可

标签

LINE 1 RNA RNA Seq

相关文章