kmer估计基因组大小
本帖最后由 生信喵 于 2022-4-5 15:59 编辑一、利用 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 准备数据
mkdir data;cd data;
#1 下载参考序列
#基因组下载
#Klebsiella pneumoniae MGH78578
#基因组： NC_009648.1https://www.ncbi.nlm.nih.gov/nuccore/NC_009648.1/
#质粒： NC_009649.1 https://www.ncbi.nlm.nih.gov/nuccore/NC_009649
efetch -db nuccore -format fasta -id NC_009648.1 > MGH78578.fasta
efetch -db nuccore -format fasta -id NC_009649 >>MGH78578.fasta
#Escherichia coli CFT073
efetch -db nuccore -format fasta -id NC_004431.1 >CFT073.fasta
#下载测序数据#https://www.ncbi.nlm.nih.gov/bioproject/PRJNA422511
#下载illumina数据prefetch SRR8482586 -O ./ -o illumina.sra
#wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-21/SRR8482586/SRR8482586.1 -O illumina.sra
#下载不到，我去本地迅雷下载了，传回服务器改名
mv SRR8482586.1 illumina.sra
chmod u+x illumina.sra
fastq-dump --split-files --gzip illumina.sra
#下载pacbio数据prefetch SRR8494912 -O ./ -o pacbio.sra
mwget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-21/SRR8494912/SRR8494912.1
mv SRR8494912.1 pacbio.sra
chmod u+x pacbio.sra
fasterq-dump pacbio.sra
gzip pacbio.fastq
#下载nanopore数据prefetch SRR8494939 -O ./ -o nanopore.sra
wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-20/SRR8494939/SRR8494939.1
mv SRR8494939.1 nanopore.sra
chmod u+x nanopore.sra
fasterq-dump nanopore.sra
gzip nanopore.fastq
# 查看数据
less -S illumina_1.fastq.gz
less -S pacbio.fastq.gz
less -S nanopore.fastq.gz
rm -f *.sra
2.2 质控
#fastqc质控
mkdir qc
fastqc -f fastq -o qc illumina_1.fastq.gz illumina_2.fastq.gz
#利用fastp进行数据过滤
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
#过滤完在质控一遍
fastqc -f fastq -o qc clean.1.fq.gz clean.2.fq.gz
2.3 jellyfish 统计 kmer
安装jellyfish和genomescope2
cd ~/biosoft
wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.11.tar.gz
tar zxvf jellyfish-1.1.11.tar.gz
mkdir jellyfish
./configure --prefix=$PWD/../jellyfish
make -j 8
make install
/share/home/xiehs/biosoft/jellyfish/bin写入环境变量
安装genomescope2
git clone https://github.com/tbenavi1/genomescope2.0.git
cd genomescope2.0/
mkdir ~/R_libs
echo "R_LIBS=~/R_libs/" >> ~/.Renviron #有写入权限
Rscript install.R#kmer估计基因组大小
mkdir 34.kmer;cd 34.kmer/
ln -s ../data/illumina_1.fastq.gz .
ln -s ../data/illumina_2.fastq.gz .
#jellyfish统计kmer
#不支持压缩格式
jellyfish count -m 15 -s 2G -o kmer15.count-s 2G -C <(zcat clean.1.fq.gz )
jellyfish stats kmer15.count_0 -o kmer15.stats
jellyfish histo kmer15.count_0 > kmer15.histo
/share/home/xiehs/biosoft/genomescope2.0写入环境变量
#genomescope2估计基因组大小
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和纳米孔。
页:
[1]