使用 GCE 进行基因组大小评估

1. GCE 简介

GCE(Genome Characteristics Estimation) 是华大基因用于基因组评估的软件，其参考文献为：Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects。下载地址：ftp://ftp.genomics.org.cn/pub/gce。

GCE 软件包中主要包含 kmer_freq_hash 和 gce 两支程序。前者用于进行 kmer 的频数统计，后者在前者的结果上进行基因组大小的准确估算。
2. GCE 下载和安装

$ wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
$ tar zxf gce-1.0.0.tar.gz

3. kmer_freq_hash 的使用

kmer_freq_hash 的常用例子：

$ /gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
-k 21 -l reads.list -a 10 -d 10 -t 24 -i 50000000 -o 0 -p species &> kmer_freq.log

kmer_freq_hash 的常用参数：

-k <17>
    设置 kmer 的大小。该值为 9~27，默认值为 17 。
-l string
    list文本文件，其中每行为一个fastq文件的路径。
-t int
    使用的线程数，默认为 1 。
-i int
    初始的 hash 表大小，默认为 1048576。该值最好设置为（kmer 的种类数 / 0.75）/ 线程数。如果基因组大小为 100M，测序了 40M 个 reads，reads 的长度为 100bp，测序错误率为 1%，kmer的大小为 21，则kmer的种类数为100M+40M*100*1%*21=940M，若使用24线程，则该参数设置为 i=940M/0.75/24=52222222。
-p string
    设置输出文件的前缀。
-o int
    是否输出 k-mer 序列。1: yes, 0: no，默认为 1 。推荐选 0 以节约运行时间。
-q int
    设置fastq文件的phred格式，默认为 64。该值可以为 33 或 63。
-c double
    设置k-mer最小的精度，该值位于 0~0.99，或为 -1。 -1 表示不对 kmer进行过滤。设置较高的精度，可以用于过滤低质量 kmer。精度是由 phred 格式的碱基质量计算得来的。
-r int
    设置获取 k-mer 使用到的 reads 长度。默认使用 reads 的全长。
-a int
    忽略read前面该长度的碱基。
-d int
    忽略read后面该长度的碱基。
-g int
    设置使用该数目的碱基来获取 k-mers，默认是使用所有的碱基来获取 k-mer。

kmer_freq_hash 的主要结果文件为 species.freq.stat。该文件有 2 列：第1列是kmer重复的次数，第二列是kmer的种类数。该文件有255行，第225行表示kmer重复次数>=255的kmer的总的种类数。该文件作为 gce 的输入文件。
kmer_freq_hash 的输出到屏幕上的信息结果保存到文件 kmer_freq.log 文件中。该文件中有粗略估计基因组的大小。其中的 Kmer_individual_num 数据作为 gce 的输入参数。
4. gce 的使用

Example:
./gce -f 17mer.freq >gce.table 2>gce.log
./gce -f 17mer.freq -c 20 -H 1 >gce.table 2>gce.log

gce 的常用例子：

$ ./gce-1.0.0/gce \
-f species.freq.stat -c 85 -g 4112118028 -m 1 -D 8 -b 1 > species.table 2> species.log

gce 的常用参数：

-f string
    kmer depth frequency file
-c int
    kmer depth frequency 的主峰对应的 depth。gce 会在该值附近找主峰。
-g int
    总共的 kmer 数。一定要设定该值，否则 gce 会直接使用 -f 指定的文件计算 kmer 的总数。由于默认下该文件中最大的 depth 为 255，因此，软件自己计算的值比真实的值偏小。同时注意该值包含低覆盖度的 kmer。
-M int
    支持最大的 depth 值，默认为 256 。
-m int
    估算模型的选择，离散型（0），连续型（1）。默认为 0，对真实数据推荐选择 1 。
-D int
    precision of expect value，默认为 1。如果选择了 -m 1，推荐设置该值为 8。
-H int
    使用杂合模式（1），不使用杂合模式（0）。默认值为 0。只有明显存在杂合峰的时候，才选择该值为 1 。
-b int
    数据是（1）否（0）有 bias。当 K > 19时，需要设置 -b 1 。

gce 的结果文件为 species.table 和 species.log 。species.log 文件中的主要内容：

raw_peak now_node low_kmer now_kmer cvg genome_size a[1] b[1]
84 35834245 22073804 4044916750 84.6637 4.83093e+07 0.928318 0.637648

raw_peak：覆盖度为 84 的 kmer 的种类数最多，为主峰。
now_node： kmer的种类数。
low_kmer：低覆盖度的 kmer 数。
now_kmer：去除低覆盖度的 kmer 数，此值 = （-g 参数指定的总 kmer 数） - low_kmer 。
cvg：估算出的平均覆盖度
genome_size：基因组大小，该值 = now_kmer / cvg 。
a[1]：在基因组上仅出现 1 次的 kmer 之种类数比例。
b[1]：在基因组上仅出现 1 次的 kmer 之数量比例。该值代表着基因组上拷贝数为 1 的序列比例。

如果使用 -H 1 参数，则会得额外得到如下信息：

for hybrid: a[1/2]=0.223671 a1=0.49108
kmer-species heterozygous ratio is about 0.125918

上面结果中，0.125918 是由 a[1/2] 计算出来的。 0.125918 = a[1/2] / （ 2- a[1/2] ) 。
a[1/2]=0.223671 表示在所有的 uniqe kmer 中，有 0.223671 比例的 kmer 属于杂合 kmer 。

此外，有 a[1/2] 和 b[1/2] 的值在最后的统计结果中。重复序列的含量 = 1 - b[1/2] - b[1] 。

则杂合率 = 0.125918 / kmer_size 。若计算出的杂合率低于 0.2%，个人认为测序数据应该是纯合的。这时候，应该不使用 -H 1 参数。使用 -H 1 参数会对基因组的大小和重复序列含量估算造成影响。

#Attention:
#If there is clear hybrid peak, you need to set -H and -c at the same time, -c 可以从输出的Kmer表格中看到峰值大概位置
#If you use the continous model, you need to set -D at the same time

非杂合模式示例

此处使用上述得到的k-mer频数统计结果文件“test.freq.stat”进行基因组特征评估。

输入文件“test.freq.stat”，该文件中，峰值对应k-mer频数等于58的位置，不再指定k-mer片段总数而是默认使用输入k-mer频数统计表文件计算k-mer总数，已知该物种基因组为单倍体，故不使用杂合模式，使用连续型估算模型且期望值精度设定为8，设置计算时所支持的最大k-mer频数256（若k-mer频数表由kmer_freq_hash得到，则可缺省，非kmer_freq_hash所得结果一定要设置该参数）。

程序运行完毕得到结果文件“test.gce.stat”和“test.gce.log”。

“test.gce.stat”为中间计算结果输出，一般情况下无需关注。

“test.gce.log”为程序运行的日志文件，同时记录了物种基因组特征评估结果。评估统计结果同样在该文件的最下方

raw_peak，k-mer频数统计结果的“主峰”频数；

now_node，k-mer的种类数；

low_kmer，低覆盖度的k-mer数；

now_kmer，过滤低覆盖度的k-mer数后的k-mer总数；

cvg，估算出的测序平均深度；

genome_size，估算出的基因组大小，genome_size = now_kmer / cvg ；

a[1]，仅出现1次的k-mer种类数占k-mer种类数总数的比值；

b[1]，仅出现1次的k-mer片段数量占k-mer片段总数量的比值；该值可用于表征基因组中拷贝数为1的序列比例，则1-b[1]可认定为重复序列比例。

杂合模式示例

若测序物种的基因组高杂合（如通过k-mer曲线判断，例如曲线中存在杂合峰），我们可考虑运行杂合模式，此时相较于非杂合模式的运行结果会更佳。当然，若测序物种的基因组仅为简单基因组（无杂合），则是不能使用杂合模式运行的，强行运行杂合模式结果将相当不可靠。

此处继续使用上述得到的k-mer频数统计结果文件“test.freq.stat”进行基因组特征评估。虽然该测序物种的k-mer曲线显示其为单倍体，不可以使用杂合模式；但此处作为示例主要演示软件的使用，因此请无需在意该操作的正确性。

输入文件“test.freq.stat”，同上述，该文件中峰值对应k-mer频数等于58的位置，不再指定k-mer片段总数而是默认使用输入k-mer频数统计表文件计算k-mer总数，此处使用杂合模式，并使用连续型估算模型且期望值精度设定为8，设置计算时所支持的最大k-mer频数256（若k-mer频数表由kmer_freq_hash得到，则可缺省，非kmer_freq_hash所得结果一定要设置该参数）。

与上述结果一致，程序运行完毕得到结果文件“test.gce.stat”和“test.gce.log”。然后我们关注“test.gce.log”最下方的统计结果。

结合其它K-mer分析工具一起使用

由于GCE的kmer_freq_hash程序统计k-mer频数时，支持的最大频数深度为225，出现次数大于255的k-mer数量会与出现次数等于255的k-mer数量合并，因此有时可能无法满足分析需求。因此，我们可以考虑将GCE结合其它k-mer分析工具一起使用。

如下示例，依然使用了本次的测试数据。首先使用JELLYFISH进行k-mer频数统计，之后将JELLYFISH的结果输入至GCE中，评估物种基因组大小、重复序列含量等信息。由于JELLYFISH支持最大k-mer频数为10000，因此我们可知，结合JELLYFISH（JELLYFISH+gce）一起分析的结果与只使用GCE（kmer_freq_hash+gce）分析的结果相比会更加准确。

其中，JELLYFISH使用说明可参见：http://blog.sciencenet.cn/blog-3406804-1161522.html

基因组特征最终评估结果如下所示。特别注意，由于JELLYFISH支持最大k-mer频数为10000，因此在此处的gce参数中，-M设定为10001。

作者：沙子_要变努力
链接：gce基因组特征评估 2019-07-24 - 简书
来源：简书

使用 GCE 进行基因组大小评估相关推荐

三大公有云托管 Kubernetes 服务 (EKS、GKE、AKS) 评估
EKS vs GKE vs AKS - Evaluating Kubernetes in the Cloud | StackRoxhttps://www.stackrox.com/post/2021/ ...
基因组特征评估——k-mer analysis
在进行基因组组装之前,对基因组进行简单的评估,以获取基因组大小.杂合度.重复序列等特征有助于我们概括性的了解将要获取的基因组的大致信息.K-mer analysis是快速对基因组的上述概况进行评估的常 ...
宏基因组组装质量评估新方法-MAGISTA
谷禾健康尽管地球上微生物类群的繁多,但只有一小部分得到了培养和有效命名.因为大多数菌无法在非常特定的条件下培养分离鉴定. 在过去十年中,宏基因组研究的重要性已经凸显,因为它能够评估细菌基因库并发现当 ...
2022-2028年中国快捷酒店行业市场全景评估及前瞻分析报告
[报告类型]产业研究 [出版时间]即时更新(交付时间约3个工作日) [发布机构]智研瞻产业研究院 [报告格式]PDF版本报告介绍了快捷酒店行业相关概述.中国快捷酒店行业运行环境.分析了中国快捷酒店行 ...
Apriori算法通俗详解_fpgrowth2_关联分析评估
20220317 https://blog.csdn.net/a790439710/article/details/103080674 支持度,置信度各指标再理解条件模式基:在某元素比如y出现的前提 ...
scikit-learn - 分类模型的评估 (classification_report)
20201225 分类报告输出到csv from sklearn.metrics import classification_report report = classification_report ...
TVM性能评估分析（七）
TVM性能评估分析(七) Figure 1. Performance Improvement Figure 2. Depthwise convolution Figure 3. Data Fusion ...
TVM性能评估分析（六）
TVM性能评估分析(六) Figure 1. The workflow of development PC, compile, deploy to the device, test, then mod ...
TVM性能评估分析（五）
TVM性能评估分析(五) Figure 3. A futher speed up with operator fusion Table 1. Performance issue of cuBLAS' ...

使用 GCE 进行基因组大小评估

使用 GCE 进行基因组大小评估相关推荐

最新文章

热门文章