用k-mer分析进行基因组调查：（四）用GenomeScope评估基因组特征

（全文约3500字）

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

# 1. k-mer进行基因组调查的软件概况
k-mer进行基因组调查分为**k-mer频数统计**和**基因组特征评估**两步。
- GenomeScope可以实现第二步基因组特征评估。
- 需要在jellyfish/KMC等软件的第一步结果k-mer频数分布表的基础上，GenomeScope才可实现。

推荐第一步获取k-mer频数分布表的命令：
1. jellyfish
```
jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o sample.jf #计算k-mer频率，生成sample.jf
jellyfish histo -t 10 sample.jf > sample.histo #生成k-mer频数直方表sample.histo和k-mer直方图
```

2. KMC
```
mkdir tmp
ls *.fastq.gz > FILES
kmc -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

1. 上传第一步获得的k-mer频数分布表sample.histo文件；
2. 设置参数Kmer length为第一步选择的k-mer长度值，这里是17；
3. 参数Read length为序列读长，一般为150；
4. 参数Max kmer coverage默认是1000。

建议按照物种情况修改，比如10000，以统计更准确。

这个参数太小，可能造成过滤过多的Kmer，导致估计的基因组大小偏小的情况。

这个参数太大则可能把高拷贝数量的DNA，比如叶绿体DNA，包括进Kmer的统计，造成GenomeScope算法的误差，所以还是不推荐使用-1或太大的值。
5. 提交后几分钟就可以得到结果，保存结果图片可用于发表。

**Figure 1. GenomeScope1.0结果示例**

# 4. GenomeScope2.0 网页版【推荐】 —— 适用于多倍体物种
[GenomeScope2.0 网页版](http://qb.cshl.edu/genomescope/genomescope2.0)：http://qb.cshl.edu/genomescope/genomescope2.0

GenomeScope2.0版本相较于1.0，进行了许多改进，主要是增加了多倍体物种的基因组调查，并提出Smudgeplot方法来估计基因组的倍性和基因组结构。

## 4.1. GenomeScope 2.0 使用步骤
1. 上传第一步获得的k-mer频数分布表histo文件；
2. 设置参数Kmer length为第一步选择的k-mer长度值，这里是17；
3. 参数倍性Ploidy根据物种的倍性设定，默认是二倍体，设置成2；
4. 参数Max k-mer coverage默认是-1，即不限制最大k-mer深度。

建议按照物种情况修改，比如10000，以统计更准确。

这个参数太小，可能造成过滤过多的Kmer，导致估计的基因组大小偏小的情况。

这个参数太大则可能把高拷贝数量的DNA，比如叶绿体DNA，包括进Kmer的统计，造成GenomeScope算法的误差，所以还是不推荐使用-1或太大的值。
5. 参数Average k-mer coverage for polyploid genome默认是-1，即不进行筛选，可以根据情况调整。
6. 提交后几分钟就可以得到结果，保存结果图片可用于发表。

## 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**，可以评估多倍体的基因组特征。

1. 四倍体共有两种可能的拓扑结构，代表着同源四倍体和异源四倍体，每种拓扑包含三种杂合基因型和一种纯合基因型，共有五种基因型。(五倍体有五种可能的拓扑，六倍体有十六种)

2. 根据结果中杂合基因型的分布模式可以区分异源四倍体和同源四倍体。

**Figure 3. 四倍体的两种拓扑结构可能性 a异源四倍体，b同源四倍体。图片来源：[GenomeScope 2.0 paper](https://www.nature.com/articles/s41467-020-14998-3)**

3. GenomeScope2.0的**四倍体**结果

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

**Figure 4. GenomeScope2.0 多倍体结果示例 a异源四倍体，b同源四倍体。图片来源： [GenomeScope 2.0 paper](https://www.nature.com/articles/s41467-020-14998-3)**

# 5. GenomeScope1.0本地版 —— 适用于二倍体物种
GenomeScope1.0的本地版是用一个R脚本实现的，在[GenomeScope github](https://github.com/schatzlab/genomescope)可以下载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. 下载和安装

```shell
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实践经验
1. 实际使用中发现，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安装
1. 依赖

依赖是[tbenavi1/KMC](https://github.com/tbenavi1/KMC)和[GenomeScope2.0](https://github.com/tbenavi1/genomescope2.0)。
- 用于统计k-mers频数的软件。建议[tbenavi1/KMC](https://github.com/tbenavi1/KMC)，里面包括一个smudge_pairs程序，用来找杂合k-mer对。也可以用jellyfish代替KMC，参考manual of smudgeplot with jellyfish：https://github.com/KamilSJaron/smudgeplot/wiki/manual-of-smudgeplot-with-jellyfish。
- [GenomeScope2.0](https://github.com/tbenavi1/genomescope2.0)

2. 安装

`conda install -c bioconda smudgeplot` #conda安装

## 8.3. Smudgeplot使用
1. 用KMC计算k-mer频率，生成k-mer直方图
```
mkdir tmp
ls *.fastq.gz > FILES
kmc -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的最大值。

2. 选择覆盖阈值
- 可以目视检查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
```

3. 提取阈值范围的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".dump
smudgeplot.py hetkmers -o kmcdb_L"$L"_U"$U" < kmcdb_L"$L"_U"$U".dump
```

4. 生成污点图(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](https://github.com/KamilSJaron/smudgeplot)**

## 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
1. GenomeScope 1.0 github：https://github.com/schatzlab/genomescope
2. GenomeScope 2.0 github：https://github.com/tbenavi1/genomescope2.0
3. GenomeScope 1.0 paper：https://academic.oup.com/bioinformatics/article/33/14/2202/3089939
4. GenomeScope 2.0 + Smudgeplot paper：https://www.nature.com/articles/s41467-020-14998-3
5. Smudgeplot github：https://github.com/KamilSJaron/smudgeplot

用k-mer分析进行基因组调查：（四）用GenomeScope评估基因组特征相关推荐

生信步骤|kmc+genomescope进行基因组调查
在组装未知基因组时,往往需要利用重测序数据提前进行基因组调查,以获取其基因组规模,杂合率,重复序列比例,GC含量等信息.从而更好地拟定后继测序策略.基因组调查可以采用kmers方法.kmers基因组调 ...
2019微生物组——16S扩增子分析专题培训第四期
文章目录课程简介课程大纲一.生信基础知识和技巧二.图表解读和绘制三.扩增子基础和分析流程四.可重复计算和统计绘图五.功能预测和机器学习六.网络和环境因子分析往期精彩回顾主讲教师助 ...
新秀nginx源代码分析数据结构篇（四）红黑树ngx_rbtree_t
新秀nginx源代码分析数据结构篇(四)红黑树ngx_rbtree_t Author:Echo Chen(陈斌) Email:chenb19870707@gmail.com Blog:Blog.csd ...
linux 查看握手时间,实战：tcpdump抓包分析三次握手四次挥手
本文档以实战的形式介绍tcpdump抓包分析三次握手四次挥手的过程. 执行tcpdump命令 tcpdump -n -i ens32 host 192.168.10.10 and 42.186.113 ...
Week2 Teamework from Z.XML - 必应缤纷桌面助手 - 软件分析与用户需求调查
软件分析与用户需求调查(2013) from Z.XML 本次团队作业要求: 通过定性, 定量地分析, 总结和评定某软件是否满足了目标用户的需求,并把分析的过程和结果用博客表达出来. 选题:必应缤纷桌 ...
EEG巨型分析I：跨研究的频谱和振幅特征
导读通过汇集多项研究的统计结果(元分析),fMRI领域取得了重大成就.最近,fMRI标准化工作的重点是实现跨研究(巨型分析)的fMRI原始数据的联合分析,以期获得更详细的见解.然而,目前尚不清楚在E ...
iMeta｜湘雅医院刘庆组-泛癌分析揭示铜死亡调节子的临床和分子特征
点击蓝字关注我们泛癌分析揭示铜死亡调节子的临床和分子特征 https://doi.org/10.1002/imt2.68 RESEARCH ARTICLE ●2022年12月5日,中南大学湘雅医院 ...
『深度学习项目四』基于ResNet101人脸特征点检测
相关文章: [深度学习项目一]全连接神经网络实现mnist数字识别 [深度学习项目二]卷积神经网络LeNet实现minst数字识别 [深度学习项目三]ResNet50多分类任务[十二生肖分类] 『深度 ...
[统计学笔记] （四）数据分布的数字特征
(四)数据分布的数字特征数据的分布特征与使用的描述统计量数据集中趋势在统计研究中,需要搜集大量数据并对其进行加工整理,大多数情况下数据都会呈现出一种钟形分布,即各个变量值与中间位置的距离越近,出 ...
Nature子刊：涵盖20多万个基因组的人体肠道微生物参考基因组集
Nature子刊:涵盖20多万个人体肠道微生物基因组的参考基因组集 A unified catalog of 204,938 reference genomes from the human gut ...

用k-mer分析进行基因组调查：（四）用GenomeScope评估基因组特征

用k-mer分析进行基因组调查：（四）用GenomeScope评估基因组特征相关推荐

最新文章

热门文章