关于经验:用kmer分析进行基因组调查四用GenomeScope评估基因组特征

（全文约3500字）

【举荐】用Smudgeplot评估物种倍性后，用组合jellyfish+GenomeScope1.0做二倍体物种的基因组考察，用组合KMC+GenomeScope2.0做多倍体物种的基因组考察。

1. k-mer进行基因组考察的软件详情

k-mer进行基因组考察分为k-mer频数统计和基因组特色评估两步。

GenomeScope能够实现第二步基因组特色评估。
须要在jellyfish/KMC等软件的第一步后果k-mer频数散布表的根底上，GenomeScope才可实现。

举荐第一步获取k-mer频数散布表的命令：

jellyfish

jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o sample.jf #计算k-mer频率，生成sample.jfjellyfish histo -t 10 sample.jf > sample.histo #生成k-mer频数直方表sample.histo和k-mer直方图

KMC

mkdir tmpls *.fastq.gz > FILESkmc -k21 -t16 -m64 -ci1 -cs10000 @FILES kmcdb tmp #计算k-mer频率kmc_tools transform kmcdb histogram sample.histo -cx10000 #生成k-mer直方图

2. GenomeScope详情

GenomeScope能够利用第一步jellyfish或KMC等其他软件剖析失去k-mer频数散布表(sample.histo文件)实现第二步基因组特色评估。
GenomeScope1.0在2017年发表，用于二倍体物种的基因组考察；2020年又发表了GenomeScope 2.0版本，用于多倍体物种的基因组考察，并公布了用于判断物种倍性的Smudgeplot。
GenomeScope1.0的基因组特色后果包含基因组大小(genome size)，杂合度(heterozygosity)，反复序列比例，GC含量等；
GenomeScope2.0的基因组特色后果包含基因组大小(genome size)，杂合度(heterozygosity)，反复序列比例，GC含量，基因型比例，和基因组构造(同源/异源多倍体)等。
GenomeScope有网页版和Linux本地版，性能一样；举荐网页版，免去装置的麻烦。

3. GenomeScope1.0 网页版【举荐】—— 实用于二倍体物种

GenomeScope1.0 网页版：http://qb.cshl.edu/genomescope

上传第一步取得的k-mer频数散布表sample.histo文件；
设置参数Kmer length为第一步抉择的k-mer长度值，这里是17；
参数Read length为序列读长，个别为150；
参数Max kmer coverage默认是1000。
倡议依照物种状况批改，比方10000，以统计更精确。
这个参数太小，可能造成过滤过多的Kmer，导致预计的基因组大小偏小的状况。
这个参数太大则可能把高拷贝数量的DNA，比方叶绿体DNA，包含进Kmer的统计，造成GenomeScope算法的误差，所以还是不举荐应用-1或太大的值。
提交后几分钟就能够失去后果，保留后果图片可用于发表。

Figure 1. GenomeScope1.0后果示例

4. GenomeScope2.0 网页版【举荐】 —— 实用于多倍体物种

GenomeScope2.0 网页版：http://qb.cshl.edu/genomescop...

GenomeScope2.0版本相较于1.0，进行了许多改良，次要是减少了多倍体物种的基因组考察，并提出Smudgeplot办法来预计基因组的倍性和基因组构造。

4.1. GenomeScope 2.0 应用步骤

上传第一步取得的k-mer频数散布表histo文件；
设置参数Kmer length为第一步抉择的k-mer长度值，这里是17；
参数倍性Ploidy依据物种的倍性设定，默认是二倍体，设置成2；
参数Max k-mer coverage默认是-1，即不限度最大k-mer深度。
倡议依照物种状况批改，比方10000，以统计更精确。
这个参数太小，可能造成过滤过多的Kmer，导致预计的基因组大小偏小的状况。
这个参数太大则可能把高拷贝数量的DNA，比方叶绿体DNA，包含进Kmer的统计，造成GenomeScope算法的误差，所以还是不举荐应用-1或太大的值。
参数Average k-mer coverage for polyploid genome默认是-1，即不进行筛选，能够依据状况调整。
提交后几分钟就能够失去后果，保留后果图片可用于发表。

4.2. GenomeScope 2.0 后果

4.2.1. 二倍体后果

二倍体的GenomeScope 2.0 后果与GenomeScope 1.0 后果的次要不同之处在于杂合度后果(het)变成了2.0版本的代表基因型的aa和ab的比例，其中杂合基因型ab的比例即为杂合度。2.0后果中的p值代表设置的物种倍性。

Figure 2. GenomeScope2.0 二倍体后果示例

4.2.2. 辨别异源四倍体和同源四倍体

GenomeScope2.0增加了参数倍性Ploidy，能够评估多倍体的基因组特色。

四倍体共有两种可能的拓扑构造，代表着同源四倍体和异源四倍体，每种拓扑蕴含三种杂合基因型和一种纯合基因型，共有五种基因型。(五倍体有五种可能的拓扑，六倍体有十六种)
依据后果中杂合基因型的分布模式能够辨别异源四倍体和同源四倍体。

Figure 3. 四倍体的两种拓扑构造可能性 a异源四倍体，b同源四倍体。图片起源：GenomeScope 2.0 paper

GenomeScope2.0的四倍体后果

在GenomeScope 2.0 的后果中，如果杂合基因型aaab的比例大于aabb，则认为该物种是异源四倍体；如果杂合基因型aaab的比例小于aabb，则认为该物种是同源四倍体。

Figure 4. GenomeScope2.0 多倍体后果示例 a异源四倍体，b同源四倍体。图片起源： GenomeScope 2.0 paper

5. GenomeScope1.0本地版 —— 实用于二倍体物种

GenomeScope1.0的本地版是用一个R脚本实现的，在GenomeScope github能够下载genomescope.R脚本，下载后把genomescope.R文件退出环境变量即可应用。。

Rscript genomescope.R sample.histo k-mer_length read_length output_dir [kmer_max] [verbose]

必须参数：

sample.histo：频数散布直方表，jellyfish的后果。
k-mer_length：k-mer长度，通常是17，21，与jellyfish统一。
read_length：reads长度，这里是150bp的PE reads，所以是150。
output_dir：输入目录，后果图和文本都输入到这个目录。

6. GenomeScope2.0本地版 —— 实用于多倍体物种

6.1. 下载和装置

git clone https://github.com/tbenavi1/genomescope2.0.git #下载cd genomescope2.0/mkdir ~/R_libs #创立主目录下的R_libs文件夹用于装置本地R库echo "R_LIBS=~/R_libs/" >> ~/.Renviron #创立/编辑.Renviron文件，使得R在创立的R_libs文件夹加载库Rscript install.R #装置

装置后把目录下的genomescope.R文件退出环境变量即可应用。

6.2. 应用

genomescope.R -i histogram_file -o output_dir -k k-mer_length

参数：

-i histogram_file：频数散布直方表，jellyfish或KMC的后果。
-k k-mer_length：k-mer长度，通常是17，21，与jellyfish/KMC的设置统一。
-o output_dir：输入目录，后果图和文本都输入到这个目录。
-p ploidy：设置倍性。
-l lambda：设置测序的均匀k-mer覆盖率的初始猜想。
-n name_prefix：设置输入文件的前缀。
-m max_kmercov：设置从剖析中排除高频k-mers的截止值，依据物种状况确定，举荐1000或10000。

7. GenomeScope实践经验

理论应用中发现，GenomeScope1.0和2.0经常估算差别较大。倡议二倍体还是应用GenomeScope1.0。
在估算一个约300Mb的二倍体基因组时，GenomeScope1.0估算进去267Mb，GenomeScope2.0估算进去149Mb。
在估算一个约6Gb的四倍体基因组时发现，GenomeScope1.0估算进去5.5Gb，GenomeScope2.0估算进去2.7Gb。

8. Smudgeplot

Smudgeplot是2020年与GenomeScope2.0一起发表的用于预计物种的倍性的软件。开发者打算接下来把Smudgeplot整合进GenomeScope。

8.1. Smudgeplot原理

Smudgeplot从k-mer数据库中提取杂合k-mer对，而后训练杂合k-mer对。

通过比拟k-mer对覆盖度的总数(CovA + CovB)和绝对覆盖度(CovB / (CovA + CovB))，统计杂合k-mers对的数量，Smudgeplot能够解析基因组构造。

8.2. Smudgeplot装置

依赖

依赖是tbenavi1/KMC和GenomeScope2.0。

用于统计k-mers频数的软件。倡议tbenavi1/KMC，外面包含一个smudge_pairs程序，用来找杂合k-mer对。也能够用jellyfish代替KMC，参考manual of smudgeplot with jellyfish：https://github.com/KamilSJaro...。
GenomeScope2.0

装置

conda install -c bioconda smudgeplot #conda装置

8.3. Smudgeplot应用

用KMC计算k-mer频率，生成k-mer直方图

mkdir tmpls *.fastq.gz > FILESkmc -k21 -t16 -m64 -ci1 -cs10000 @FILES kmcdb tmp #计算k-mer频率kmc_tools transform kmcdb histogram sample.histo -cx10000 #生成k-mer频数直方表sample.histo和k-mer直方图

kmc命令参数：

-k21：k-mer长度设置为21
-t16：线程16
-m64：内存64G，设置应用RAM的大抵数量，范畴1-1024。
-ci1 -cs10000：统计k-mer coverages覆盖度范畴在[1-10000]的。
@FILES：保留了输出文件列表的文件名为FILES
kmcdb：KMC数据库的输入文件名前缀
tmp：长期目录

kmc_tools命令参数：

-cx10000：贮存在直方图文件中counter的最大值。

抉择笼罩阈值
能够目视查看k-mer直方图，抉择笼罩阈值上(U)下(L)限。
也能够用命令预计笼罩阈值上(U)下(L)限。L的取值范畴是[20-200]，U的取值范畴是[500-3000]。

L=$(smudgeplot.py cutoff kmcdb_k21.hist L)U=$(smudgeplot.py cutoff kmcdb_k21.hist U)echo $L $U # these need to be sane values

提取阈值范畴的k-mers，计算k-mer pairs
用kmc_tools提取k-mers，而后用KMC的smudge_pairs计算k-mer pairs。

smudge_pairs比smudgeplot.py hetkmers应用更少内存，速度更快地寻找杂合k-mer pairs。

kmc_tools transform kmcdb -ci"$L" -cx"$U" reduce kmcdb_L"$L"_U"$U"smudge_pairs kmcdb_L"$L"_U"$U" kmcdb_L"$L"_U"$U"_coverages.tsv kmcdb_L"$L"_U"$U"_pairs.tsv > kmcdb_L"$L"_U"$U"_familysizes.tsv

如果没有KMC，能够用kmc_dump提取k-mers，而后用smudgeplot.py hetkmers计算k-mer pairs。

kmc_tools transform kmcdb -ci"$L" -cx"$U" dump -s kmcdb_L"$L"_U"$U".dumpsmudgeplot.py hetkmers -o kmcdb_L"$L"_U"$U" < kmcdb_L"$L"_U"$U".dump

生成污点图(smudgeplot)

smudgeplot.py plot kmcdb_L"$L"_U"$U"_coverages.tsv

生成两个根底的污点图，一个log尺度，一个线性尺度。

8.4. Smudgeplot后果

热度图，横坐标是绝对覆盖度 (CovB / (CovA + CovB)) ，纵坐标是总覆盖度 (CovA + CovB) ，色彩是k-mer对的频率。
每个单倍型构造都在图上出现一个"污点(smudge)"，污点的热度示意单倍型构造在基因组中呈现的频率，频率最高的单倍型构造即为预测的物种倍性后果。(比方这个图提供了三倍体的证据，AAB的频率最高)

**Figure 5. Smudgeplot污点图。
图片起源： Smudgeplot github**

8.5. Smudgeplot的KMC后果用于GenomeScope进行基因组考察

通过KMC取得的频数散布表后果可用于GenomeScope进行基因组考察

Rscript genomescope.R kmcdb_k21.hist <k-mer_length> <read_length> <output_dir> [kmer_max] [verbose]

9. references

GenomeScope 1.0 github：https://github.com/schatzlab/...
GenomeScope 2.0 github：https://github.com/tbenavi1/g...
GenomeScope 1.0 paper：https://academic.oup.com/bioi...
GenomeScope 2.0 + Smudgeplot paper：https://www.nature.com/articl...
Smudgeplot github：https://github.com/KamilSJaro...