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

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

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

1. 已编译推荐

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

2. 未编译

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命令。

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

1. 命令

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

2. 参数

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

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

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

2. kmer_freq.log：操作日志文件

• 旧版本的日志文件还简要统计了测序数据。
• 本文件底部统计k-mer片段总数、k-mer种类数、k-mer平均频率，碱基总数，reads粗略估计平均长度、基因组大小等信息。

1. 获取k-mer总数

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

用于gce的-g参数

2. 获取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数量高是由于总数超过该列的所有频数造成的，两端都可以忽略。

1. 纯合模式

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

2. 杂合模式

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

3. 开发者建议

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

然后再用-c参数和-H 1参数运行混合模式。最后，比较两种模式的结果，以确定哪一种更适合当前数据。

-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。

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

1. 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_ker_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 的序列比例。

2. 如果使用杂合模式-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]

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