（全文约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.gztar -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

  1. 参数
  2. -k 17：k-mer size
  3. -t 24：线程
  4. -p prefix：输入文件的前缀
  5. 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参数

  1. 获取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

  1. 杂合模式

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

  1. 开发者倡议

对于不能判断纯合还是杂合的数据，能够先运行纯合模式，取得初始峰值(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 的序列比例。
  1. 如果应用杂合模式-H 1，则会在gce.log文件最初额定失去上面信息：
  2. a[1/2]：a[1/2]=0.223671示意在所有的 uniqe kmers 品种中，有 0.223671 比例的 kmer 属于杂合 kmer 。
  3. b[1/2]：a[1/2]=0.326934 示意在所有的 uniqe kmers 个体中，有 0.326934 比例的 kmer 属于杂合 kmer 。
    通过计算，还能够取得的信息：
  4. 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=...