0

用k-mer分析进行基因组调查(genome survey) —— 用GCE分步实现

 3 months ago
source link: https://yanzhongsino.github.io/2022/06/07/omics_genome.survey_GCE/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

【推荐】用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 下载

3.2. GCE 安装

  1. git clone [email protected]: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 #编译

文件夹下有kmerfreqgce两个命令。
老版本是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

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

4.1.2. 输出文件

  1. sample.kmer.freq.stat:结果文件
  • 结果文件sample.kmer.freq.stat,记录了每个k-mer频数的统计信息,用于生成GCE的输入文件。
  • 之前的版本,该文件只统计到第255行,第255行之后的数据合并至第255行,表示k-mer出现频数>=255的片段总数。现在这个版本(gce-1.0.2)上限变成了65534。
  1. kmer_freq.log:运行log文件
  • kmer_freq.log程序运行的日志文件。
  • 老版本的日志文件还对测序数据进行了简要统计。在该文件的最下方,统计了k-mer片段总数、k-mer种类数、k-mer平均频数、碱基总数、reads平均长度、基因组大小的粗略估计等信息。

4.2. 基因组特征评估

4.2.1. 获取参数

  1. 获取k-mer总数
    less sample.kmer.freq.stat | grep "#Kmer indivdual number" #获取k-mer总数,用于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 #获取k-mer深度分布表
  • 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 的序列比例。
  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] ) 。
  • 重复序列的含量 = 1 - b[1/2] - b[1] 。

5. references


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK