本帖最后由 生信喵 于 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.1 https://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
- #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
- 安装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和纳米孔。
| 
发表于 2022-4-5 15:22:38
