kmer估计基因组大小

生信喵

一、利用 kmer 估计基因组大小
      要想估计基因组的大小，也就是整条基因组的长度，我们把这个值设为大 G。那么测序的所有碱基数可以计算出来，将所有 reads 的碱基加起来就可以，为大 S。用所有碱基数除以每个碱基的平均覆盖深度 D，碱基总数除以测序深度，那么就可以得到基因组的长度了。所以，要想估计基因组大小，必须计算出每个位点被覆盖的平均深度，因为我们已经有了总碱基数S。但是这个深度无法直接计算出来，所以，我们通过 kmer 的深度，来推测测序的深度，进而求出基因组大小。那么就是要推测出 kmer 深度与测序深度之间的关系，下面我们来看一下如何通过 kmer 的深度来计算测序的深度。
      看下面的公式。
      我们设 G 为实际基因组大小，k 为 kmer 长度，l 为 reads 长度，n_k为所有 kmer 的个数，n_b为所有碱基的个数，d_k为 kmer 期望深度，d_b为碱基期望深度，其中n_k和d_k皆可以统计出来。因为 Kmer 深度符合泊松分布，所以 d_k 即为 Kmer 深度的平 均值，而整个基因组约可产生 G 个 Kmer，这里面忽略边界效应。最终我们就得到了如下的公式。

G = n_k / d_k = n_b / d_b ;

      其中，只有 dk 一个未知量，我们只要得到 kmer 的深度就可以推测出基因组的大小。kmer 的深度，通过统计就可以得到，这样就得到了基因组的大小。


二、实战
      参考文章PMID:31483244,细菌基因组
      文章20个样本，每个样本分别使用二代、pacbio和纳米孔，总共60个sra。
2.1 准备数据

1. mkdir data;cd data;
   
2. #1 下载参考序列
   
3. #基因组下载
   
4. #Klebsiella pneumoniae MGH78578
   
5. #基因组： NC_009648.1  https://www.ncbi.nlm.nih.gov/nuccore/NC_009648.1/
   
6. #质粒： NC_009649.1 https://www.ncbi.nlm.nih.gov/nuccore/NC_009649
   
7. efetch -db nuccore -format fasta -id NC_009648.1 > MGH78578.fasta
   
8. efetch -db nuccore -format fasta -id NC_009649 >>MGH78578.fasta
   
9. #Escherichia coli CFT073
   
10. efetch -db nuccore -format fasta -id NC_004431.1 >CFT073.fasta
    
11. #下载测序数据#https://www.ncbi.nlm.nih.gov/bioproject/PRJNA422511
    
12. #下载illumina数据prefetch SRR8482586 -O ./ -o illumina.sra
    
13. #wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-21/SRR8482586/SRR8482586.1 -O illumina.sra
    
14. #下载不到，我去本地迅雷下载了，传回服务器改名
    
15. mv SRR8482586.1 illumina.sra
    
16. chmod u+x illumina.sra
    
17. fastq-dump --split-files --gzip illumina.sra
    
18. #下载pacbio数据prefetch SRR8494912 -O ./ -o pacbio.sra
    
19. mwget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-21/SRR8494912/SRR8494912.1
    
20. mv SRR8494912.1 pacbio.sra
    
21. chmod u+x pacbio.sra
    
22. fasterq-dump pacbio.sra
    
23. gzip pacbio.fastq
    
24. #下载nanopore数据prefetch SRR8494939 -O ./ -o nanopore.sra
    
25. wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-20/SRR8494939/SRR8494939.1
    
26. mv SRR8494939.1 nanopore.sra
    
27. chmod u+x nanopore.sra
    
28. fasterq-dump nanopore.sra
    
29. gzip nanopore.fastq
    
30. # 查看数据
    
31. less -S illumina_1.fastq.gz
    
32. less -S pacbio.fastq.gz
    
33. less -S nanopore.fastq.gz
    
34. rm -f *.sra

复制代码

2.2 质控

1. #fastqc质控
   
2. mkdir qc
   
3. fastqc -f fastq -o qc illumina_1.fastq.gz illumina_2.fastq.gz
   
4. 
   
5. #利用fastp进行数据过滤  
   
6. fastp -i illumina_1.fastq.gz -I illumina_2.fastq.gz -o clean.1.fq.gz -O clean.2.fq.gz -D -z 4 -q 20 -u 30 -n 10 -f 20 -t 10 -F 20 -T 10 -h clean.html
   
7. #过滤完在质控一遍
   
8. fastqc -f fastq -o qc clean.1.fq.gz clean.2.fq.gz

复制代码

2.3 jellyfish 统计 kmer

1. 安装jellyfish和genomescope2
   
2. cd ~/biosoft
   
3. wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.11.tar.gz
   
4. tar zxvf jellyfish-1.1.11.tar.gz
   
5. mkdir jellyfish
   
6. ./configure --prefix=$PWD/../jellyfish
   
7. make -j 8
   
8. make install
   
9. /share/home/xiehs/biosoft/jellyfish/bin写入环境变量
   
10. 
    
11. 安装genomescope2
    
12. git clone https://github.com/tbenavi1/genomescope2.0.git
    
13. cd genomescope2.0/
    
14. mkdir ~/R_libs
    
15. echo "R_LIBS=~/R_libs/" >> ~/.Renviron #有写入权限
    
16. Rscript install.R

复制代码

1. #  kmer估计基因组大小
   
2. mkdir 34.kmer;cd 34.kmer/
   
3. ln -s ../data/illumina_1.fastq.gz .
   
4. ln -s ../data/illumina_2.fastq.gz .
   
5. #jellyfish统计kmer
   
6. #不支持压缩格式
   
7. jellyfish count -m 15 -s 2G -o kmer15.count  -s 2G -C <(zcat clean.1.fq.gz )
   
8. jellyfish stats kmer15.count_0 -o kmer15.stats
   
9. jellyfish histo kmer15.count_0 > kmer15.histo
   
10. 
    
11. /share/home/xiehs/biosoft/genomescope2.0写入环境变量
    
12. 
    
13. #genomescope2估计基因组大小
    
14. genomescope2 -i kmer15.histo -o gscope -p 1 -k 15

复制代码

      gscope文件夹中linear_plot.png就可以看到基因组大小约5.989m
      
      transformed_linear_plot.png调整后的，还有取对数值的。
      以上对illumina1号测序，取kmer15bp的参数比较出。MGH78578.fasta的大小是5.49m
      我们可以更改kmer为17，用illumina2号再跑一遍流程，能对回5.4最好。一般用二代的估算。也可以用我们下载的pacbio和纳米孔。