美文网首页基因组组装
基因组survey——K-mer频谱

基因组survey——K-mer频谱

作者: 徐诗芬 | 来源:发表于2021-02-03 21:24 被阅读0次

    Kmer是从测序数据中滑窗提取出的长度为k的寡聚核苷酸序列,可以评估基因组大小、杂合度、重复序列比例等。
    在测序reads均匀分布的前提下,有以下公式:
    基因组长度 = 总碱基数 / 平均测序深度 = 总kmer数 / 平均kmer深度
    根据该公式,可以使用总kmer数和平均kmer深度估算基因组长度。K值使用满足以下公式计算的最大奇数:4 ^ K / genome > 200

    做kmer有几种软件:SOAPdenovo2的Kmerfreq模块,jellyfish和KAT(K-mer Analysis Toolkit)

    jellyfish:简书1简书2GitHub
    产生的kmer频数分布数据需要用R包GenomeScopefindGSE来进行统计估计。(http://qb.cshl.edu/genomescope/)

    这里我还是用一种后来出现但整合了jellyfish和其他分析的软件KAT( is accomplished through an integrated and modified version of Jellyfish2's counting method)
    https://doi.org/10.1093/bioinformatics/btw663
    Github: https://github.com/TGAC/KAT
    Document:https://kat.readthedocs.io/en/latest/using.html

    安装

    1. 从bioconda里安装最新版本的kat:
    bioconda install kat
    #这里我用以下这种简单快捷的方法:
    conda install -c bioconda kat
    conda activate
    kat -h
    
    1. 从GitHub安装:
      需要先安装很多依赖包,不然可能很多功能用不了:

    -GCC V4.8+
    -make
    -autoconf V2.53+
    -automake V1.11+
    -libtool V2.4.2+
    -pthreads (probably already installed)
    -zlib
    -Python V3.5+ with the tabulate, scipy, numpy and matplotlib packages and C API installed. This is optional but highly recommended, without python KAT functionality is limited: no plots, no distribution analysis, and no documentation.
    -Sphinx-doc V1.3+ (Optional: only required for building the documentation.

    然后是一系列的安装编译过程:

    git clone https://github.com/TGAC/KAT.git
    cd KAT
    ./build_boost.sh
    ./autogen.sh
    ./configure
    make
    

    二、运行

    1. 给定一个k值(kmer = 17,小基因组一般用17,大基因组才用默认的27),产生不同kmers数量的直方图
      usage:kat hist [options] (<input>)+
      options: -o 输出文件 [kat.hist];-t [1]线程数; -l [1] histogram的最低值;-h [10000] histogram的最高值;-m [27] kmer的长度
    nohup kat hist -o KAT_result/kat.hist -t 16 <(gzip -d -c NGS_10G_?.fq.gz) 2> kat_hist.log & #不支持压缩格式,需先解压
    nohup kat hist -m 17 -o KAT_result/kat.hist  -t 16 NGS_10G_1.clean.fq NGS_10G_2.clean.fq &  #kat1.hist.png
    nohup kat hist -m 17 -o KAT_result/kat.hist  -h 400 -t 16 NGS_10G_1.clean.fq NGS_10G_2.clean.fq &  #kat2.hist.png,缩短x轴看得更清楚测序深度
    
    kat1.hist.png
    kat2.hist.png
    1. Kmer的GC覆盖图
      usage: kat gcp (<input>)+
      options: -o 输出文件 [kat.hist];-t [1]线程数; -x [1] 当创建污染矩阵时,用于gc数据的bins的数量; -y [1000] 当创建污染矩阵时,用于coverage数据的bins的数量; -m [27] kmer的长度
    nohup kat gcp -m 17 -o KAT_result/kat.gcp -t 12 NGS_10G_1.clean.fq NGS_10G_2.clean.fq &
    
    kat.gcp.mx.png

    SOAPdenovo2的Kmerfreq模块

    ls NGS_10G_?.clean.fq > NGS10.clean.lst
    KmerFreq_AR -k 17 -t 12 -m 1 -p KmerFreq.10G -q 33 NGS10.clean.lst > NGS.10G.kmer.count 2>NGS.10G.kmerfreq.cerr

    freq.cz和.freq.cz.len文件是用于Corrector_AR对reads进行校正分析

    Corrector_AR -k 17 -t 12 -l 3 -Q 33 KmerFreq.10G.freq.cz KmerFreq.10G.freq.cz.len NGS10.clean.lst >corr.cout 2>corr.cerr ##没有校正的reads就不需要再跑下一步
    ls NGS_10G_1.clean.fq.cor.pair_1.fq.gz NGS_10G_2.clean.fq.cor.pair_2.fq.gz > kmer.input.fil.data
    KmerFreq_AR -k 17 -t 12 -m 1 -p KmerFreq.10G.cor -q 33 kmer.input.fil.data > NGS.10G.kmer.cor.count 2>NGS.10G.kmerfreq.cor.cerr
    

    jellyfish和genomescope

    #/usr/bin/jellyfish    jellyfish 1.1.10
    nohup /usr/bin/jellyfish count -m 17 -s 1G -t 12 -C -o jellyfish.out <(gzip -d -c NGS_10G_?.fq.gz) jellyfish.out.log & 
    nohup /usr/bin/jellyfish count -m 17 -s 1G -t 12 -C -o jellyfish.out NGS_10G_1.fq NGS_10G_2.fq 2> jellyfish.out.log & 
    /usr/bin/jellyfish histo -t 12 jellyfish.out_0 > jellyfish.histo 
    ~/Software/genomescope/genomescope.R jellyfish.histo 17 149 genomescope.result
    
    plot.png

    相关文章

      网友评论

        本文标题:基因组survey——K-mer频谱

        本文链接:https://www.haomeiwen.com/subject/agfjaktx.html