二倍体线虫基因组组装实战（2）——全基因组Survey

29 Oct 2024 6645字 23分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

全基因组 Survey

基因组测序数据组装之前，我们还要做一个全基因组survey。主要是为了减少盲目性，先做低深度的基因组分析，也是初步了解物种基因组特征的有效方法，比如评估基因组大小和杂合情况，为后续全基因组de novo组装策略指定提供指导。

全基因组Survey 分析：

  • 基因组Survey基于小片段文库的低深度测序数据( 50X-100X ) ;
  • 通过K-mer分析 ,有效的评估基因组大小、GC含量、杂合度以及重复序列的含量等信息;
  • 是全面了解某一物种基因组特征的有效方法;
  • 为后续的全基因组 de novo 测序的组装策略的制定提供理论依据。

1

基因组复杂程度的经验性标准：

  • 简单基因组: 单倍体；或纯合二倍体；或杂合度低于0.5%, 且重复序列低于50%, 且GC含量在35%-65%的二倍体。
  • 复杂基因组: 杂合度在0.5%~1.2%之间，或重复序列高于50%，或GC含量异常(<35%或>65%)的二倍体，或者多倍体。复杂基因组可以采用“2+3”即二代和三代测序技术相结合，加之Hi-C辅助组装的组装策略。
  • 高复杂基因组: 杂合度>1.2%；或重复序列占比大于65%。

有条件的话，也可以用流式细胞仪对基因组大小做个预估。我这里只有二代基因组测序数据，因此用基因组二代测序数据做全基因组survey。当然，这里要注意一点，做全基因组survey的样本和后续de novo组装的样本要来自同一个个体，避免个体间基因组特征的差异。

补充知识（一）： 基因组复杂程度预估

2

补充知识（二）：根据K-mer图确认物种倍型

二倍体：杂合峰：主峰：重复峰 = 1:2:4（比值为横坐标峰的比值）

3

三倍体：正常情况下杂合峰：主峰：重复峰 = 1:2:3（左图）。主峰和重复峰深度低则可能重叠在一起：杂合峰：主峰：重复峰 = 1:2（右图）]

4

异源四倍体：2个峰，呈现1:2的关系

5

同源四倍体：同源四倍体的峰就是1 : 2 : 3 : 4 ,其中3和4经常重叠在一起

7

1.原始数据质控

如果是二代数据，可以进行一下数据的质控，因为我这边是三代HIFI数据，其本身可以直接用于组装，就不质控了，下面给出了二代数据质控的大致流程。

1.1 fastqc生成质控报告

使用fastqc来生成质检报告

fastqc *.fq.gz -o ./

1.2 trim-galore数据过滤

过滤一遍，把末端接头adapter序列过滤掉。

trim_galore -q 25 -phred33 -length 100 -stringency 1 -paired -o clean_data 1_raw_1.fq.gz 1_raw_2.fq.gz

软件的参数可以参考这位博主的博客转录组数据分析笔记（1）——如何用fastqc和trim-galore做测序数据质控

2 k-mer分析

K-mer 概念

  • 定义：K-mer 是指将测序读取（reads）拆分成长度为 k 的序列片段。例如，对于核酸，若 k=17，则会提取包含 17 个碱基的序列。K-mer 分析常用于评估基因组特征及复杂性。

  • 基因组大小估计：

    • 基因组大小可以通过公式估算： \(\text{基因组大小} = \frac{\text{总 K-mer 数量}}{\text{K-mer 期望测序深度}}\)

  • K-mer 分布曲线特征：

    • 主峰位置：K-mer 分布曲线的主峰横坐标通常对应期望的测序深度。
    • 理想条件：在测序覆盖均匀、无误差及无重复序列的情况下，K-mer 分布曲线符合 泊松分布。

  • 基因组类型的 K-mer 分布特征：

    • 单倍体/纯合基因组：K-mer 分布曲线通常只有一个主峰。

    • 杂合二倍体基因组

      ：K-mer 分布曲线有两个峰：

      • 杂合峰：主峰一半处，深度约为主峰深度的一半。
      • 纯合峰：主峰位置。

    • 重复序列：当基因组中重复序列含量较高时，主峰后会形成重复峰（在主峰的两倍处）或拖尾。

  • 选择 K 值：

    • 一般选择 17-mer 进行基因组大小评估，因为长度为 17 的核酸序列，理论上有 4174^{17}417 种组合，足以覆盖正常基因组的复杂性。
    • 为避免回文序列的影响，K-mer 分析中选择的 K 值通常为奇数。

K-mer 分析是评估基因组特征的重要工具，通过对 K-mer 的分布特征进行分析，可以推测基因组大小、测序深度及基因组的杂合性、重复性等信息。这种方法的有效性和准确性使其在基因组学研究中得到了广泛应用。

2.1 安装 Jellyfish

K-mer 分析是一个计算资源密集型的过程，尤其是处理大规模的测序数据（如十几个 GB 的 reads）。如果我们自己编写脚本实现 K-mer 分析，需要将数据分割成不同长度的片段并统计每个 K-mer 的出现次数，这不仅耗时且容易出错。因此，使用专门的软件工具来进行 K-mer 分析是更为高效的选择。

Jellyfish 是一款广泛使用的 K-mer 统计软件，它具有以下优势：

  • 运行速度快：能够迅速处理大规模的 DNA 序列。
  • 内存消耗低：设计优化，适合在资源有限的环境中使用。
  • 支持并行计算：可以利用多核处理器提高处理效率。

安装注意事项

尽管可以通过 conda 直接安装 Jellyfish，但根据我的经验，这种方式可能会导致长时间无结果的情况。我曾尝试使用 conda 安装 Jellyfish，结果运行了整整一天仍然没有输出，也没有任何报错。为了排除参数设置的问题，我决定从 GitHub 上重新下载源代码，进行编译和安装。最终，使用这种方法我在不到 10 分钟内就得到了结果。这让我对两种安装方式的差异感到疑惑，因此在这里记录下我的经历，以帮助其他用户。

数据处理

需要注意的是，Jellyfish 不支持直接处理 .gz 压缩文件，因此在进行 K-mer 分析之前，必须使用 gunzip 命令将过滤后的 clean reads 解压缩。

conda install -c bioconda jellyfish  # 可以用conda安装，我运行的时候出了问题，暂未解决，不推荐

jellyfish的github下载地址

我们用本地安装的方式，先下载tar.gz的源码包，tar -zxvf解压后进入jellyfish-2.3.0文件夹。

我是集群登录的，下面讲的步骤都是在集群上操作（非root账户）

# 第一步检测。本质上是一个shell脚本，根据系统环境产生合适的makefile文件或者C的头文件（.h结尾的文件），非root账户下--prefix后面接上自己账户的绝对路径。
./configure --prefix=/public/home/wlxie
# 第二步编译。对源代码包进行编译，如果有错误自己看是否有依赖库的缺失，主要是这个问题。
make
# 第三步安装。如果前面没有指定自己账户的路径，这一步是会报错没有权限的（用户不能向系统目录写入文件）。
make install
# 第四步自检。
make check

动态链接库问题及解决方法

在使用 make 和 make check 安装软件时，可能会遇到由于动态链接库命名不同而导致的错误，尤其是在使用集群运行程序时，某些动态库模块无法调用。以下是我总结的解决方法，以帮助您顺利解决这些问题。

环境变量配置

在运行 configure 脚本后，软件会在您的家目录下生成 bin、lib 和 share 目录。这里需要特别关注的是 bin 和 lib 目录：

  • bin 目录包含可执行文件，您需要将其路径添加到 PATH 环境变量中，以便能够运行相关命令。
  • lib 目录包含动态链接库，您需要将其路径添加到 LD_LIBRARY_PATH 环境变量中，以确保程序能够找到所需的库文件。

为了配置这些环境变量，您可以在家目录下的 .bashrc 文件中加入以下内容：

shell复制代码export PATH="/public/home/wlxie/bin:$PATH"
export LD_LIBRARY_PATH="/public/home/wlxie/lib:$LD_LIBRARY_PATH"

保存并退出后，使用以下命令刷新系统环境变量：

source ~/.bashrc

解决动态库找不到的问题

在我遇到的报错中，libcrypto.so.1.0.0 和 libstdc++.so.6 这两个动态库无法找到。通过 locate 命令，我发现这两个库文件在系统目录 /lib64/ 下是存在的。为了解决问题，我将这两个动态库文件直接复制到家目录的 lib 文件夹中，结果所有问题都解决了。

检查动态库版本

如果遇到 libstdc++.so.6 报错某版本的文件不存在，可以在动态库目录下使用 strings 命令查看动态库中是否有对应的文件版本：

strings /lib64/libstdc++.so.6 | grep CXXABI

例如，如果找不到 GLIBCXX_3.4.26，可以检查动态库中是否存在这个版本。如果不存在，您可能需要更新动态库；如果存在但找不到，建议直接将其拷贝到自己的 lib 目录下。

检查安装结果

执行 make check 后，系统会生成一个日志文件 test-suite.log。如果没有任何项目出现失败（fail），则说明软件安装成功，安装过程中没有问题。

8

2.2 k-mer频数分布

# k-mer计数
jellyfish count -m 17 -s 100G -t 100 -C -o 17-mer.hifi.jf ../../Cae_elegans_data/hifi_Cae_elegans.fq

jellyfish count: 这是 Jellyfish 的计数命令，用于统计 k-mer 的出现次数。

-m 17: 指定 k-mer 的长度为 17 个碱基。这意味着在读取序列中提取长度为 17 的所有可能子串作为 k-mer。

-s 100G: 分配 100 GiB 的内存给 Jellyfish 用于计数。这是 Jellyfish 在整个计算过程中可以使用的总内存量。

-t 100: 指定使用 100 个线程 来并行处理数据。使用多个线程可以加快计算速度。

-C: 表示生成的 k-mer 计数是 可重用的（canonical），即将重复的 k-mer 进行合并，以减少内存使用。

-o 17-mer.hifi.jf: 指定输出文件的名称为 17-mer.hifi.jf，这个文件将保存 k-mer 计数的结果。

../../Cae_elegans_data/hifi_Cae_elegans.fq: 这是输入文件的路径，指向要进行 k-mer 计数的 HiFi 数据文件。在这个例子中，是一个 FASTQ 格式的文件，包含了 HiFi 测序的数据。

# k-mer频数统计
jellyfish histo -t 4 17-mer.hifi.jf > 17-mer.hifi.histo
# 统计总k-mer数和特征k-mer数等
jellyfish stats 17-mer.hifi.jf -o 17-mer.hifi.counts_stats.txt

jellyfish histo: 这是 Jellyfish 的命令，用于生成 k-mer 频数的直方图。

-t 4: 指定使用 4 个线程 来并行处理统计任务。多线程可以加快处理速度。

17-mer.hifi.jf: 这是之前生成的 k-mer 计数文件，包含了长度为 17 的 k-mer 及其计数。

>: 将命令的输出重定向到指定文件。

17-mer.hifi.histo: 指定输出文件的名称为 17-mer.hifi.histo，该文件将包含 k-mer 的频数直方图数据。

jellyfish stats: 这是 Jellyfish 的命令，用于统计 k-mer 计数文件中的各种统计信息。

17-mer.hifi.jf: 这是之前生成的 k-mer 计数文件。

-o 17-mer.hifi.counts_stats.txt: 指定输出文件的名称为 17-mer.hifi.counts_stats.txt，该文件将包含 k-mer 计数的统计信息。

k-mer频数统计是在计数结果文件上进一步统计各个k-mer出现的次数，频数统计结果文件17-mer.histo将k-mer从1统计到10000，最后一行是10001以后对应的总频次。17-mer.hifi.counts_stats.txt是总的统计结果，包括k-mer总数（Total），特异的k-mer数目（Distinct）只出现过一次的k-mer数量（Unique），频数最高的k-mer数量（Max_count）四项。

9

有了频数统计结果文件17-mer.hifi.histo就可以用R作图了：

# 读取 k-mer 数据
kmer <- read.table('17-mer.hifi.histo')

# 筛选只取 5-500bp 长度的 k-mer 统计频次
kmer <- subset(kmer, V1 >= 5 & V1 <= 1000)

# 提取 k-mer 频率和对应的数量
Frequency <- kmer$V1
Number <- kmer$V2

# 计算总的 k-mer 数量和期望测序深度
total_kmer_count <- sum(Frequency * Number)    # 总的 k-mer 数量
expected_depth <- sum(Frequency * Number) / sum(Number)  # 期望测序深度

# 输出统计信息到控制台
cat("Total k-mer count:", total_kmer_count, "\n")
cat("Expected sequencing depth:", expected_depth, "\n")

# 绘制并保存 k-mer 分布图
png('kmer_plot.png')
plot(Frequency, Number, type = 'l', col = 'blue', xlab = "Frequency", ylab = "Number",
     main = "k-mer Frequency Distribution")
dev.off()  # 保存图像

10

2.3 基因组大小、重复率、杂合率估算

​ 横坐标表示 k-mer 深度，纵坐标表示 k-mer 数量。通过 k-mer 分布曲线，可以判断样本的遗传特征。我测的这个基因组数据这边应该是杂合峰和纯合峰重叠了，可能是我这边粗略直接画的原因，有些测序错误的序列没有去除，所以区分不是很明显，不过应该影响不大。运行后Total k-mer count: 28852166052 ，Expected sequencing depth: 233.7991，Total repeated k-mer count: 123405820。

  • k-mer分布曲线未中出现明显异常峰，说明二代测序提取的DNA纯度还不错
  • 根据（K-mer 数量）/（K-mer 期望测序深度）估算基因组大小约为123.45M。
  • 根据（重复序列的k-mer总数）/（K-mer 期望测序深度）估计重复序列大小为528 kb，即重复率0.428%
  • 要计算杂合率，需要统计非重复k-mer的总数，也就是计算杂合峰面积，建议还是用软件或者在线工具比如genomescope2.0

jellyfish + GenomeScope是一套应用非常广泛的基因组survey方法，GenomeScope2.0适合用于分析二倍体生物。

11

以上是 GenomeScope2.0 的网页版界面，只需提供由 jellyfish 生成的 .histo 结果文件，并设置相关参数即可进行分析：

  • k-mer length: 输入 k-mer 的长度。
  • Ploidy: 指定染色体的倍性。
  • Max k-mer coverage: 默认值为 -1，表示不限制最大 k-mer 深度。
  • Average k-mer coverage for polyploid genome: 默认值为 -1，表示不进行筛选。我这里限制了10000

提交这些参数后，几分钟内即可生成可用于发表的图表和报告。

12

可以看到估计的基因组大小约为 59.07 M，杂合率为 0.496%，而杂合峰的覆盖度（深度）为 227。以下是对这个图的解读：

  • 蓝色区域: 表示实际观测到的 k-mer 分布值。
  • 黑色拟合线: 表示在去除错误（errors）后剩余的 k-mer 分布，认为这些是正确的数据，基于此评估基因组大小。
  • 黄色拟合线: 代表非重复区域的 k-mer 分布，反映理想情况下的预期分布。
  • 橙色拟合线区域: 表示低深度的错误 k-mer，这些 k-mer 被认为是由于测序错误而引入的。
  • 黑色虚线: 显示 k-mer 分布中的几个峰值。

估计的基因组大小比之前的估计小，原因在于去除错误的标准不同。之前的标准仅是简单去除了 k-mer 深度为 1-4 的错误序列，而这里采用了构建模型的方法来选择错误序列，因此更为准确。

在网页版的结果部分，还提供了总的统计结果，这样可以方便地计算重复率，整体信息一目了然，具体细节这里就不再赘述。

13

本文参考：

0基础学习基因组三代测序组装（3）——全基因组Survey)

Survey数据质控、NT比对