关于经验:用kmer分析进行基因组调查四用GenomeScope评估基因组特征

7次阅读

共计 7100 个字符，预计需要花费 18 分钟才能阅读完成。

（全文约 3500 字）

【举荐】用 Smudgeplot 评估物种倍性后，用组合 jellyfish+GenomeScope1.0 做二倍体物种的基因组考察，用组合 KMC+GenomeScope2.0 做多倍体物种的基因组考察。

1. k-mer 进行基因组考察的软件详情

k-mer 进行基因组考察分为 k-mer 频数统计 和基因组特色评估 两步。

  • GenomeScope 能够实现第二步基因组特色评估。
  • 须要在 jellyfish/KMC 等软件的第一步后果 k -mer 频数散布表的根底上，GenomeScope 才可实现。

举荐第一步获取 k -mer 频数散布表的命令：

  1. jellyfish

    jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o sample.jf #计算 k -mer 频率，生成 sample.jf
jellyfish histo -t 10 sample.jf > sample.histo #生成 k -mer 频数直方表 sample.histo 和 k -mer 直方图

  2. KMC

    mkdir tmp
ls *.fastq.gz > FILES
kmc -k21 -t16 -m64 -ci1 -cs10000 @FILES kmcdb tmp #计算 k -mer 频率
kmc_tools transform kmcdb histogram sample.histo -cx10000 #生成 k -mer 直方图

2. GenomeScope 详情

  • GenomeScope 能够利用第一步 jellyfish 或 KMC 等其他软件剖析失去 k -mer 频数散布表 (sample.histo 文件) 实现第二步 基因组特色评估。
  • GenomeScope1.0 在 2017 年发表，用于二倍体物种的基因组考察；2020 年又发表了 GenomeScope 2.0 版本，用于多倍体物种的基因组考察，并公布了用于判断物种倍性的 Smudgeplot。
  • GenomeScope1.0的基因组特色后果包含基因组大小(genome size)，杂合度(heterozygosity)，反复序列比例，GC 含量等；
  • GenomeScope2.0的基因组特色后果包含基因组大小 (genome size)，杂合度(heterozygosity)，反复序列比例，GC 含量，基因型比例，和基因组构造(同源 / 异源多倍体) 等。
  • GenomeScope 有网页版和 Linux 本地版，性能一样；举荐网页版，免去装置的麻烦。

3. GenomeScope1.0 网页版【举荐】—— 实用于二倍体物种

GenomeScope1.0 网页版：http://qb.cshl.edu/genomescope

  1. 上传第一步取得的 k -mer 频数散布表 sample.histo 文件；
  2. 设置参数 Kmer length 为第一步抉择的 k -mer 长度值，这里是 17；
  3. 参数 Read length 为序列读长，个别为 150；

  4. 参数 Max kmer coverage 默认是 1000。

    倡议依照物种状况批改，比方 10000，以统计更精确。

    这个参数太小，可能造成过滤过多的 Kmer，导致预计的基因组大小偏小的状况。

    这个参数太大则可能把高拷贝数量的 DNA，比方叶绿体 DNA，包含进 Kmer 的统计，造成 GenomeScope 算法的误差，所以还是不举荐应用 - 1 或太大的值。

  5. 提交后几分钟就能够失去后果，保留后果图片可用于发表。

<img src=”https://github.com/yanzhongsino/yanzhongsino.github.io/blob/hexo/source/images/omics_genome.survey_GenomeScope1.0.png?raw=true” width=80% title=”GenomeScope1.0 后果示例 ” align=center/>

<p align=”center”>Figure 1. GenomeScope1.0 后果示例 </p>

4. GenomeScope2.0 网页版【举荐】—— 实用于多倍体物种

GenomeScope2.0 网页版：http://qb.cshl.edu/genomescop…

GenomeScope2.0 版本相较于 1.0，进行了许多改良，次要是减少了多倍体物种的基因组考察，并提出 Smudgeplot 办法来预计基因组的倍性和基因组构造。

4.1. GenomeScope 2.0 应用步骤

  1. 上传第一步取得的 k -mer 频数散布表 histo 文件；
  2. 设置参数 Kmer length 为第一步抉择的 k -mer 长度值，这里是 17；
  3. 参数倍性 Ploidy 依据物种的倍性设定，默认是二倍体，设置成 2；

  4. 参数 Max k-mer coverage 默认是 -1，即不限度最大 k -mer 深度。

    倡议依照物种状况批改，比方 10000，以统计更精确。

    这个参数太小，可能造成过滤过多的 Kmer，导致预计的基因组大小偏小的状况。

    这个参数太大则可能把高拷贝数量的 DNA，比方叶绿体 DNA，包含进 Kmer 的统计，造成 GenomeScope 算法的误差，所以还是不举荐应用 - 1 或太大的值。

  5. 参数 Average k-mer coverage for polyploid genome 默认是 -1，即不进行筛选，能够依据状况调整。
  6. 提交后几分钟就能够失去后果，保留后果图片可用于发表。

4.2. GenomeScope 2.0 后果

4.2.1. 二倍体后果

二倍体的 GenomeScope 2.0 后果与 GenomeScope 1.0 后果的次要不同之处在于杂合度后果 (het) 变成了 2.0 版本的代表基因型的 aa 和 ab 的比例，其中杂合基因型 ab 的比例即为杂合度。2.0 后果中的 p 值代表设置的物种倍性。

<img src=”https://github.com/yanzhongsino/yanzhongsino.github.io/blob/hexo/source/images/omics_genome.survey_GenomeScope2.0.png?raw=true” width=80% title=”GenomeScope2.0 二倍体后果示例 ” align=center/>

<p align=”center”>Figure 2. GenomeScope2.0 二倍体后果示例 </p>

4.2.2. 辨别 异源四倍体 和同源四倍体

GenomeScope2.0 增加了参数倍性Ploidy，能够评估多倍体的基因组特色。

  1. 四倍体共有两种可能的拓扑构造，代表着同源四倍体和异源四倍体，每种拓扑蕴含三种杂合基因型和一种纯合基因型，共有五种基因型。(五倍体有五种可能的拓扑，六倍体有十六种)
  2. 依据后果中杂合基因型的分布模式能够辨别异源四倍体和同源四倍体。

<img src=”https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41467-020-14998-3/MediaObjects/41467_2020_14998_Fig2_HTML.png?as=webp” width=80% title=” 四倍体的两种拓扑构造可能性 ” align=center/>

<p align=”center”>Figure 3. 四倍体的两种拓扑构造可能性 a 异源四倍体，b 同源四倍体。图片起源：GenomeScope 2.0 paper</p>

  1. GenomeScope2.0 的 四倍体 后果

在 GenomeScope 2.0 的后果中，如果杂合基因型 aaab 的比例大于 aabb，则认为该物种是异源四倍体；如果杂合基因型 aaab 的比例小于 aabb，则认为该物种是同源四倍体。

<img src=”https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41467-020-14998-3/MediaObjects/41467_2020_14998_Fig6_HTML.png?as=webp” width=100% title=”GenomeScope2.0 多倍体后果示例 ” align=center/>

<p align=”center”>Figure 4. GenomeScope2.0 多倍体后果示例 a 异源四倍体，b 同源四倍体。图片起源：GenomeScope 2.0 paper</p>

5. GenomeScope1.0 本地版 —— 实用于二倍体物种

GenomeScope1.0 的本地版是用一个 R 脚本实现的，在 GenomeScope github 能够下载 genomescope.R 脚本，下载后把 genomescope.R 文件退出环境变量即可应用。。

Rscript genomescope.R sample.histo k-mer_length read_length output_dir [kmer_max] [verbose]

必须参数：

  • sample.histo：频数散布直方表，jellyfish 的后果。
  • k-mer_length：k-mer 长度，通常是 17，21，与 jellyfish 统一。
  • read_length：reads 长度，这里是 150bp 的 PE reads，所以是 150。
  • output_dir：输入目录，后果图和文本都输入到这个目录。

6. GenomeScope2.0 本地版 —— 实用于多倍体物种

6.1. 下载和装置

git clone https://github.com/tbenavi1/genomescope2.0.git #下载
cd genomescope2.0/
mkdir ~/R_libs #创立主目录下的 R_libs 文件夹用于装置本地 R 库
echo "R_LIBS=~/R_libs/" >> ~/.Renviron #创立 / 编辑.Renviron 文件，使得 R 在创立的 R_libs 文件夹加载库
Rscript install.R #装置

装置后把目录下的 genomescope.R 文件退出环境变量即可应用。

6.2. 应用

genomescope.R -i histogram_file -o output_dir -k k-mer_length

参数：

  • -i histogram_file：频数散布直方表，jellyfish 或 KMC 的后果。
  • -k k-mer_length：k-mer 长度，通常是 17，21，与 jellyfish/KMC 的设置统一。
  • -o output_dir：输入目录，后果图和文本都输入到这个目录。
  • -p ploidy：设置倍性。
  • -l lambda：设置测序的均匀 k -mer 覆盖率的初始猜想。
  • -n name_prefix：设置输入文件的前缀。
  • -m max_kmercov：设置从剖析中排除高频 k -mers 的截止值，依据物种状况确定，举荐 1000 或 10000。

7. GenomeScope 实践经验

  1. 理论应用中发现，GenomeScope1.0 和 2.0 经常估算差别较大。倡议二倍体还是应用 GenomeScope1.0。
  2. 在估算一个约 300Mb 的二倍体基因组时，GenomeScope1.0 估算进去 267Mb，GenomeScope2.0 估算进去 149Mb。
  3. 在估算一个约 6Gb 的四倍体基因组时发现，GenomeScope1.0 估算进去 5.5Gb，GenomeScope2.0 估算进去 2.7Gb。

8. Smudgeplot

Smudgeplot 是 2020 年与 GenomeScope2.0 一起发表的用于预计物种的倍性的软件。开发者打算接下来把 Smudgeplot 整合进 GenomeScope。

8.1. Smudgeplot 原理

Smudgeplot 从 k -mer 数据库中提取杂合 k -mer 对，而后训练杂合 k -mer 对。

通过比拟 k -mer 对覆盖度的总数 (CovA + CovB) 和绝对覆盖度(CovB / (CovA + CovB))，统计杂合 k -mers 对的数量，Smudgeplot 能够解析基因组构造。

8.2. Smudgeplot 装置

  1. 依赖

依赖是 tbenavi1/KMC 和 GenomeScope2.0。

  • 用于统计 k -mers 频数的软件。倡议 tbenavi1/KMC，外面包含一个 smudge_pairs 程序，用来找杂合 k -mer 对。也能够用 jellyfish 代替 KMC，参考 manual of smudgeplot with jellyfish：https://github.com/KamilSJaro…。
  • GenomeScope2.0
  1. 装置

conda install -c bioconda smudgeplot #conda 装置

8.3. Smudgeplot 应用

  1. 用 KMC 计算 k -mer 频率，生成 k -mer 直方图

    mkdir tmp
ls *.fastq.gz > FILES
kmc -k21 -t16 -m64 -ci1 -cs10000 @FILES kmcdb tmp #计算 k -mer 频率
kmc_tools transform kmcdb histogram sample.histo -cx10000 #生成 k -mer 频数直方表 sample.histo 和 k -mer 直方图

kmc 命令参数：

  • -k21：k-mer 长度设置为 21
  • -t16：线程 16
  • -m64：内存 64G，设置应用 RAM 的大抵数量，范畴 1 -1024。
  • -ci1 -cs10000：统计 k -mer coverages 覆盖度范畴在 [1-10000] 的。
  • @FILES：保留了输出文件列表的文件名为 FILES
  • kmcdb：KMC 数据库的输入文件名前缀
  • tmp：长期目录

kmc_tools 命令参数：

  • -cx10000：贮存在直方图文件中 counter 的最大值。
  1. 抉择笼罩阈值
  2. 能够目视查看 k -mer 直方图，抉择笼罩阈值上 (U) 下(L)限。
  3. 也能够用命令预计笼罩阈值上 (U) 下(L)限。L 的取值范畴是[20-200]，U 的取值范畴是[500-3000]。
L=$(smudgeplot.py cutoff kmcdb_k21.hist L)
U=$(smudgeplot.py cutoff kmcdb_k21.hist U)
echo $L $U # these need to be sane values
  1. 提取阈值范畴的 k -mers，计算 k -mer pairs
  2. 用 kmc_tools 提取 k -mers，而后用 KMC 的 smudge_pairs 计算 k -mer pairs。

  3. smudge_pairs比 smudgeplot.py hetkmers 应用更少内存，速度更快地寻找杂合 k -mer pairs。

    kmc_tools transform kmcdb -ci"$L" -cx"$U" reduce kmcdb_L"$L"_U"$U"
smudge_pairs kmcdb_L"$L"_U"$U" kmcdb_L"$L"_U"$U"_coverages.tsv kmcdb_L"$L"_U"$U"_pairs.tsv > kmcdb_L"$L"_U"$U"_familysizes.tsv

  • 如果没有 KMC，能够用 kmc_dump 提取 k -mers，而后用 smudgeplot.py hetkmers 计算 k -mer pairs。

    kmc_tools transform kmcdb -ci"$L" -cx"$U" dump -s kmcdb_L"$L"_U"$U".dump
smudgeplot.py hetkmers -o kmcdb_L"$L"_U"$U" < kmcdb_L"$L"_U"$U".dump
  1. 生成污点图(smudgeplot)

smudgeplot.py plot kmcdb_L"$L"_U"$U"_coverages.tsv

生成两个根底的污点图，一个 log 尺度，一个线性尺度。

8.4. Smudgeplot 后果

  • 热度图，横坐标是绝对覆盖度 (CovB / (CovA + CovB))，纵坐标是总覆盖度 (CovA + CovB)，色彩是 k -mer 对的频率。
  • 每个单倍型构造都在图上出现一个 ” 污点(smudge)”，污点的热度示意单倍型构造在基因组中呈现的频率，频率最高的单倍型构造即为预测的物种倍性后果。(比方这个图提供了三倍体的证据，AAB 的频率最高)

<img src=”https://user-images.githubusercontent.com/8181573/45959760-f1032d00-c01a-11e8-8576-ff0512c33da9.png” width=70% title=”Smudgeplot 污点图 ” align=center/>

**<p align=”center”>Figure 5. Smudgeplot 污点图。
图片起源：Smudgeplot github</p>**

8.5. Smudgeplot 的 KMC 后果用于 GenomeScope 进行基因组考察

  • 通过 KMC 取得的频数散布表后果可用于 GenomeScope 进行基因组考察

Rscript genomescope.R kmcdb_k21.hist <k-mer_length> <read_length> <output_dir> [kmer_max] [verbose]

9. references

  1. GenomeScope 1.0 github：https://github.com/schatzlab/…
  2. GenomeScope 2.0 github：https://github.com/tbenavi1/g…
  3. GenomeScope 1.0 paper：https://academic.oup.com/bioi…
  4. GenomeScope 2.0 + Smudgeplot paper：https://www.nature.com/articl…
  5. Smudgeplot github：https://github.com/KamilSJaro…

正文完
经验
发表至： 经验
2022-06-16
 0
关于经验:我在京东做研发-京东云算法科学家解析爆火的ChatGPT
关于经验:2021助理高薪拉钩产品经理高薪训练营
关于经验:用kmer分析进行基因组调查二用jellyfish进行kmer频数统计
关于经验:干货自己开发了一个starter如何轻松嵌入到使用不同版本的springboot框架的对接方