（全文约2000字）

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

1. k-mer进行基因组调查的软件

k-mer进行基因组调查分为k-mer频数统计和基因组特征评估两步。

GCE可以分步实现两步。第一步k-mer频数统计和第二步基因组特征评估。
GCE第一步的结果sample.histo可以用在GenomeScope和其他基因组特征评估的软件上，实现第二步。

2. GCE 简介

GCE (genomic charactor estimator)是华大基因在2013年开发的一款基于贝叶斯模型的用于基因组调查的软件，在2020年发布了版本v2。

3. GCE 安装

3.1. GCE 下载

GCE主要托管在BGI的ftp站点：ftp://ftp.genomics.org.cn/pub/gce。
也可以在GCE github：https://github.com/fanagislab/GCE上找到。

3.2. GCE 安装

已编译【推荐】

git clone git@github.com:fanagislab/GCE.git

未编译

wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.2.tar.gz
tar -xzvf gce-1.0.2.tar.gz #解压缩和解包
make #编译

文件夹下有kmerfreq和gce两个命令。

老版本是kmer_freq_hash代替kmerfreq命令。

4. GCE 运行

GCE软件里包含两个主要的命令，kmerfreq用来完成第一步k-mer频数统计，gce用来完成第二步基因组特征评估。

4.1. k-mer频数统计

4.1.1. 运行

命令

/path/gce-1.0.2/kmerfreq -k 17 -t 24 -p sample input.path &> kmer_freq.log

参数

-k 17：k-mer size
-t 24：线程
-p prefix：输出文件的前缀
input.path：输入数据的路径保存在input.path文本文件

4.1.2. 输出文件

sample.kmer.freq.stat：结果文件

结果文件sample.kmer.freq.stat，记录了每个k-mer频数的统计信息，用于生成GCE的输入文件。
之前的版本，该文件只统计到第255行，第255行之后的数据合并至第255行，表示k-mer出现频数>=255的片段总数。
现在这个版本(gce-1.0.2)上限变成了65534。

kmer_freq.log：运行日志文件

老版本的日志文件还对测序数据进行了简要统计。
在该文件的最下方，统计了k-mer片段总数、k-mer种类数、k-mer平均频数、碱基总数、reads平均长度、基因组大小的粗略估计等信息。

4.2. 基因组特征评估

4.2.1. 获取参数

获取k-mer总数

less sample.kmer.freq.stat | grep "#Kmer indivdual number"

用于gce的-g参数

获取k-mer深度分布表

less sample.kmer.freq.stat | perl -ne 'next if(/^#/ || /^\s/); print; ' | awk '{print $1"\t"$2}' > sample.kmer.freq.stat.2colum

sample.kmer.freq.stat.2colum文件包含两列数据，第一列是k-mer频数，第二列是频数对应的k-mer的种类数量。
预期第二列数据是泊松分布，有明显峰值；频数最小的(比如小于5)的那几列对应的k-mer数量高是测序错误造成的，频数最大的那列对应的k-mer数量高是把所有大于该列的频数进行合计数量造成的，两端都可以忽略。

4.2.2. 运行gce

纯合模式

gce -f sample.kmer.freq.stat.2colum -g 173854609857 > gce.table 2> gce.log

杂合模式

gce -f sample.kmer.freq.stat.2colum -g 173854609857 -H 1 -c 75 > gce2.table 2> gce2.log

开发者建议

对于不能判断纯合还是杂合的数据，可以先运行纯合模式，获得初始峰值(raw_peak)，即k-mer的期望深度，用作-c参数。

然后再用-c参数和-H 1参数运行杂合模式。最后比较两种模式的结果，从而判断哪种更适用当前数据。

4.2.3. gce的参数

-f和-g是必需参数，其他都为可选参数。

-f sample.freq.stat.2colum：k-mer频数分布表。
-g 173854609857：k-mer片段总数，通过上面命令获取，或者查看kmer_freq.log获取。
-H 1：默认是0。homozygous mode 纯合模式(0)，heterozgyous mode 杂合模式(1)。
-c 75：独特的k-mer的期望深度，通过肉眼检查sample.gce.table的峰值获取。
-b 1：有bias(1)，无(0)bias，默认是0。
-m 1：评估模型，离散模型(0)，连续模型(1)，默认是0。
-M 1500：设置最大深度，默认1500，大于1500的会被忽略。
-D 1：期望值的精度，默认是1。

4.2.4. 结果

生成两个文件：gce.table和gce.log。gce.table保存了用于作图的数据，gce.log日志文件最后记录了基因组特征评估结果的统计。

gce.log文件最后记录了如下内容：

raw_peak        effective_kmer_species  effective_kmer_individuals      coverage_depth  genome_size     a[1]    b[1]
75      742400596       168346645871    75.8021 2.22087e+09     0.663012        0.271515

含义：

raw_peak：覆盖度为 75 的 kmer 的种类数最多，为主峰。
effective_kmer_species：真实的k-mer种类的总数(去除测序错误造成的低频k-mers)
effective_kmer_individuals：真实的k-mer个体的总数(去除测序错误造成的低频k-mers)
coverage_depth：估算出的真实k-mers的覆盖深度
genome_size：基因组大小。
genome_size = effective_kmer_individuals / coverage_depth
a[1]： uniqe kmers (在基因组上仅出现 1 次的 kmer ) 的种类数占总种类数的比例。
b[1]： uniqe kmers (在基因组上仅出现 1 次的 kmer ) 的个体数占总个体数的比例。该值代表着基因组上拷贝数为 1 的序列比例。

如果使用杂合模式-H 1，则会在gce.log文件最后额外得到下面信息：

a[1/2]：a[1/2]=0.223671表示在所有的 uniqe kmers 种类中，有 0.223671 比例的 kmer 属于杂合 kmer 。
b[1/2]：a[1/2]=0.326934 表示在所有的 uniqe kmers 个体中，有 0.326934 比例的 kmer 属于杂合 kmer 。
通过计算，还可以获得的信息：
k-mer种类的杂合率 kmer-species heterozygous ratio = 0.125918，0.125918 是由 a[1/2] 计算出来的。

0.125918=a[1/2]/（2−a[1/2])0.125918 = a[1/2] / （ 2- a[1/2] )0.125918=a[1/2]/（2−a[1/2])

重复序列的含量 = 1−b[1/2]−b[1]1 - b[1/2] - b[1]1−b[1/2]−b[1]

5. references

GCE github：https://github.com/fanagislab/GCE
kmerfreq：https://github.com/fanagislab/kmerfreq
GCE paper：https://arxiv.org/abs/1308.2012
GCE blog：http://www.chenlianfu.com/?p=2335

用k-mer分析进行基因组调查：（五）用GCE分步实现相关推荐

分析全基因组上的蛋白信息
背景在ncbi上查一个细菌的全基因组,搜取全基因组上的这五种信息:gene,locus_tag,note,location,product,蛋白序列的前10个氨基酸(为什么爬这个最后讲),并分析这个 ...
生信步骤|kmc+genomescope进行基因组调查
在组装未知基因组时,往往需要利用重测序数据提前进行基因组调查,以获取其基因组规模,杂合率,重复序列比例,GC含量等信息.从而更好地拟定后继测序策略.基因组调查可以采用kmers方法.kmers基因组调 ...
好程序员教程分析Vue学习笔记五
好程序员教程分析Vue学习笔记五,上次我们学习了Vue的组件,这次我们来学习一下路由的使用.在Vue中,所谓的路由其实跟其他的框架中的路由的概念差不多,即指跳转的路径. 注意:在Vue中,要使用路由, ...
iMeta视频教程 | StrainPanDA分析宏基因组共存菌株的组成和基因成分谱
点击蓝字关注我们 StrainPanDA宏基因组分析软件超详细教程 StrainPanDA主要针对宏基因组纵向数据,能同时实现共存菌株的组成和基因成分谱的获取,方便用户快速建立功能和组成之间的关联 ...
Week2 Teamework from Z.XML - 必应缤纷桌面助手 - 软件分析与用户需求调查
软件分析与用户需求调查(2013) from Z.XML 本次团队作业要求: 通过定性, 定量地分析, 总结和评定某软件是否满足了目标用户的需求,并把分析的过程和结果用博客表达出来. 选题:必应缤纷桌 ...
Bioinformatics BIB｜港城大孙燕妮组用于识别和分析宏基因组数据中噬菌体序列的网站...
PhaBOX: 用于识别和分析宏基因组数据中噬菌体序列的网站 PhaBOX: a server for identifying and characterizing phage contigs in ...
【JUC】JDK1.8源码分析之ConcurrentLinkedQueue（五）
一.前言接着前面的分析,接下来分析ConcurrentLinkedQueue,ConcurerntLinkedQueue一个基于链接节点的无界线程安全队列.此队列按照 FIFO(先进先出)原则对元素 ...
谭浩强《C程序分析》（第五版）第七章
为什么要用函数模块化程序设计可以实现编好一批常用的函数来实现不同的功能,例如sin函数实现求一个数的正弦值,用abs函数实现求一个数的绝对值,把他们保存在函数库中.需要用时,直接在程序中写上sin ...
【自用·记录】产品类补充学习·SaaS/visio/竞品分析SWOT用户体验五要素
目录一.SaaS学习 SaaS标准化意义: 降低边际成本(研发成本.推广和执行成本) 实现产品价值 1.降低边际成本 2.沉淀最佳实践 3.积累产品能力 SaaS标准化策略标准化设计难点标准化 ...
必应输入法的分析与用户需求调查报告
0.引言为响应微软团队提供的用户体验调查,我们远航1617团队分成三个小组,分别对功能.体验.自选这三部分内容对微软必应输入法进行了调查统计与分析建议. 我们团队自认为此次调查的人群比较有代表性,包 ...

用k-mer分析进行基因组调查：（五）用GCE分步实现