关于github:用kmer分析进行基因组调查五用GCE分步实现

（全文约 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… 上找到。

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

4. GCE 运行

GCE 软件里蕴含两个次要的命令，kmerfreq用来实现第一步 k -mer 频数统计，gce用来实现第二步基因组特色评估。

4.1. k-mer 频数统计

4.1.1. 运行

1. 命令

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

2. 参数
3. -k 17：k-mer size
4. -t 24：线程
5. -p prefix：输入文件的前缀
6. input.path：输出数据的门路保留在 input.path 文本文件

4.1.2. 输入文件

1. sample.kmer.freq.stat：后果文件
2. 后果文件 sample.kmer.freq.stat，记录了每个 k -mer 频数的统计信息，用于生成 GCE 的输出文件。
3. 之前的版本，该文件只统计到第 255 行，第 255 行之后的数据合并至第 255 行，示意 k -mer 呈现频数 >=255 的片段总数。
4. 当初这个版本 (gce-1.0.2) 下限变成了 65534。
5. kmer_freq.log：运行日志文件
6. 老版本的日志文件还对测序数据进行了简要统计。
7. 在该文件的最下方，统计了 k -mer 片段总数、k-mer 品种数、k-mer 均匀频数、碱基总数、reads 均匀长度、基因组大小的粗略预计等信息。

4.2. 基因组特色评估

4.2.1. 获取参数

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 数量高是把所有大于该列的频数进行共计数量造成的，两端都能够疏忽。

4.2.2. 运行 gce

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 参数运行杂合模式。最初比拟两种模式的后果，从而判断哪种更实用以后数据。

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 日志文件最初记录了基因组特色评估后果的统计。

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_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 的序列比例。

2. 如果应用杂合模式-H 1，则会在 gce.log 文件最初额定失去上面信息：
3. a[1/2]：a[1/2]=0.223671 示意在所有的 uniqe kmers 品种中，有 0.223671 比例的 kmer 属于杂合 kmer。
4. b[1/2]：a[1/2]=0.326934 示意在所有的 uniqe kmers 个体中，有 0.326934 比例的 kmer 属于杂合 kmer。
   通过计算，还能够取得的信息：
5. k-mer 品种的杂合率 kmer-species heterozygous ratio = 0.125918，0.125918 是由 a[1/2] 计算出来的。

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

• 反复序列的含量 = $$1 – b[1/2] – b[1]$$

5. references

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