基因组Survey 是对基因组的大小、杂合情况进行评估。我们de novo组装一个物种的基因组之前，一般要对所测物种基因组大小有一个估计，以便后续测三代数据的时候，确定大概需要测多少数据量（深度）以满足自身的组装要求。
公司一般会采用两种方式来综合估计基因组的大小：一是用流式细胞仪测DNA含量，根据一个已知基因组大小的物种为参考，分别测两个物种的DNA含量，来推算另一个物种的基因组大小；二是通过k-mer分析来估计基因组大小，通过k-mer估计基因组大小的原理这篇文章说得蛮通俗清楚(使用K-mer估计基因组大小 )，大致的原理就是:比如一个基因组大小为1000bp，如果用K=19bp的大小去逐碱基的去生成kmer，可以生成1000-19+1=981个唯一的19kmer序列，这个时候我们可以估计这个基因组大小为981bp，可想而知这种估计方式会随着基因组越大估计的误差越小。K的大小在20左右就能唯一对应到基因组的特定位置，和pcr引物大小一样。
##这里有几个实际情况需要考虑：
1#我们实际用于估计的测序数据并不是只测了1次，基因组每个位置他都测了多次（就是我们说的测序深度），我们用测序数据里的所有reads拆分得到kmer总数除以我们的测序深度就有得到了我们的基因组大小，这和上面说的1X的原理是一致的；
2我这里说的测序深度指的是绝大部分kmer的所测到/拆分到的kmer数量，事实上有序测序的偏向性会使得基因组各个的测序深度并不完全相同，所以我们的我们kmer频率分布图应该是一个泊松分布而非集中在一个深度，故我们取kmer数量最多的那个位置所对应的深度作为我们的测序深度来计算即可；
3#此外我们测序的那个个体，其重复序列，杂合度，倍性都会影响我们估计结果。
4#测序数据的准确性也会影响结果，一般选取二代测序数据或者pacbio的HIFI（CCS）数据等高质量测序数据，pacbio的CLR数据以及ont数据估计出的结果可能与实际大小想去甚远。
进行k-mer分析估计基因组大小大致分为两步：先用jellyfish生成kmer频数表，基于生成的kmer频数表使用GenomeScope2计算基因组大小。
#jellyfish生成kmer频数表
1.安装jellyfish
###手动安装
conda install -c bioconda jellyfish
wget https://github.com/gmarcais/Jellyfish/releases/download/v2.3.0/jellyfish-2.3.0.tar.gz
tar -zxvf jellyfish-2.3.0.tar.gz ###解压
mkdir jellyfish ###新建一个目录用于安装软件
cd jellyfish-2.3.0
./configure --prefix=$PWD/../jellyfish ###配置安装路径，不要在当前路目录下安装，和win下安装软件指定路径一样
make -j 4 ###编译
make install ###安装
###自动安装
conda install -c bioconda jellyfish
2.解压用于生成kmer的过滤后的.fastq.gz文件或者生成解压命令文件
###生成.fastq.gz解压脚本文件
mkdir jellyfish.out
ls ../fastp/liver.clean.CLEAN.R*.fastq.gz | awk '{print "gzip -dc "$0 }' > generate.file
###解压.fastq.gz
gunzip ../fastp/liver.clean.CLEAN.R*.fastq.gz
3.用测序文件.fastq进行kmer计算 jellyfish count
/newlustre/home/jfgui/pengfang/wt/2.3/jellyfish/bin/jellyfish count \
-t 8 \ ###线程数
-C \ ###help里注释“Count both strand, canonical representation”，大致意思是测序reads可能是方向互补的，划分的kmer要合并计算
-m 21 \ ### kmer大小
-s 1G \ # help里注释“*Initial hash size”，哈希大小，不太懂 \笑哭
-g generate.file \ ### 解压脚本文件，用于指定输入的测序文件。这里可以直接用.fstaq文件名代替
-G 2 \ ###解压并行任务数
-o kmer21 ###输出kmer计算文件
4.kmer统计 jellyfish stats
/newlustre/home/jfgui/pengfang/wt/2.3/jellyfish/bin/jellyfish stats -o kmer21.stat kmer21 ### -0输出统计文件名+kmer文件
##查看结果
cat kmer21.stat
Unique: 1218936212 ###kmer只出现一次的kmer数
Distinct: 1955488355 ###kmer种数
Total: 40783077884 ###总kmer数(包括重复)
Max_count: 123674750 ###重复次数（depth）的kmer其次数
5.生成kmer频数表 jellyfish histo
/newlustre/home/jfgui/pengfang/wt/2.3/jellyfish/bin/jellyfish histo \
-v \ ###Verbose，日志信息
-t 4 \ ###线程数
-h 123674750 \ ### High count value of histogram，要≥ kmer21.stat里的Max_count值，设置小了确实会使基因组估计偏小
-o kmer21.histo \ ### 输出文件名
kmer21 ### kmer文件
##查看结果
$head -30 kmer21.histo
1 1218936212
2 99999484
3 20935249
4 8552091
5 4228422
6 2768580
7 2005243
8 1628162
9 1406755
10 1287171
11 1230652
12 1221860
13 1237472
14 1283545
15 1350410
16 1439661
17 1543166
18 1657552
19 1794506
20 1937804
21 2095848
22 2271979
23 2467633
24 2689394
25 2941122
26 3224341
27 3555004
28 3933669
29 4373807
30 4876443 ###前面深度较低的一些kmer一般是由于测序错误导致，一些碱基测错导致这些kmer只出现了一次或几次
$tail kmer21.histo
2440713 1
2519432 1
2530141 1
2547146 1
8864574 1
8886826 1
8890472 1
9021202 1
9138838 1
10000001 4 ###后面深度较高的一些kmer可能是一些重复序列
总之大致趋势就是先极速下降后面就是一个单峰
#GenomeScope2估计基因组大小
1#安装GenomeScope2
git clone https://github.com/tbenavi1/genomescope2.0.git
cd genomescope2.0/
Rscript install.R
2#运行程序
mkdir genomescope.k21.out ###建个目录用于放输出文件
singularity exec ../../container/Assembly202306.sif genomescope.R \ ###用的基因课genek的服务器,实验室的超算没装上
-i kmer21.histo \ ### 输入kmer频数表文件
-o genomescope.k21.out \ ### 输出文件目录
-p 2 \ ### 所测物种倍型
-k 21 \ ### kmer大小
-m 123674750 ### MAX_KMERCOV，kmer最大深度，即kmer21.stat里统计的Max_count值
3#查看估计结果
输出结果有4个png文件以及3个txt文档，这里主要看linear_plot.png和summary.txt这两个够了，
###输出结果有4个png文件以及3个txt文档，这里主要看linear_plot.png和summary.txt这两个够了，
$ll genomescope.k21.out/
-rw-r--r-- 1 ug1603 VIP 197622 Jul 28 22:40 linear_plot.png
-rw-r--r-- 1 ug1603 VIP 196220 Jul 28 22:40 log_plot.png
-rw-r--r-- 1 ug1603 VIP 565 Jul 28 22:40 model.txt
-rw-r--r-- 1 ug1603 VIP 349 Jul 28 22:40 progress.txt
-rw-r--r-- 1 ug1603 VIP 592 Jul 28 22:40 summary.txt
-rw-r--r-- 1 ug1603 VIP 212302 Jul 28 22:40 transformed_linear_plot.png
-rw-r--r-- 1 ug1603 VIP 204044 Jul 28 22:40 transformed_log_plot.png
###查看summary.txt, 基因组大小估计771Mb
cat summary.txt
GenomeScope version 2.0
input file = kmer21.histo
output directory = genomescope.k21.out
p = 2
k = 21
max_kmercov = 123674750
property min max
Homozygous (aa) 99.7436% 99.762%
Heterozygous (ab) 0.237972% 0.256351%
Genome Haploid Length 770,470,462 bp 771,346,780 bp
Genome Repeat Length 227,617,245 bp 227,876,133 bp
Genome Unique Length 542,853,217 bp 543,470,648 bp
Model Fit 71.4851% 97.4163%
Read Error Rate 0.195431% 0.195431%
###后面设置了kmer大小为17和19，19结果大差不差，k=17时只有384Mb，不知道时那里出问题了明明一样的操作
cat summary.txt
GenomeScope version 2.0
input file = kmer19.histo
output directory = genomescope.k19.out
p = 2
k = 19
max_kmercov = 134518629
property min max
Homozygous (aa) 99.7264% 99.7508%
Heterozygous (ab) 0.249247% 0.273628%
Genome Haploid Length 774,637,484 bp 775,663,725 bp
Genome Repeat Length 266,158,351 bp 266,510,958 bp
Genome Unique Length 508,479,133 bp 509,152,767 bp
Model Fit 68.5259% 97.4533%
Read Error Rate 0.18645% 0.18645%
########################################################################
cat summary.txt
GenomeScope version 2.0
input file = kmer17.histo
output directory = genomescopek17_out
p = 2
k = 17
max_kmercov = 146388861
property min max
Homozygous (aa) 85.5287% 91.7584%
Heterozygous (ab) 8.24157% 14.4713%
Genome Haploid Length 383,488,328 bp 384,340,999 bp
Genome Repeat Length 197,557,658 bp 197,996,920 bp
Genome Unique Length 185,930,670 bp 186,344,079 bp
Model Fit 61.2391% 90.5093%
Read Error Rate 0.237057% 0.237057%
输出的png图片结果：
linear_plot.png
