用k-mer分析进行基因组调查：（四）用GenomeScope评估基因组特征-锐单电子商城

(全文约3500字)

**【推荐】用Smudgeplot在评估物种倍性后，使用组合jellyfish GenomeScope1.对二倍体物种进行基因组调查，使用组合KMC GenomeScope2.对多倍体物种进行基因组调查。**

# 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.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。

建议根据物种情况进行修改，如1万，统计更准确。这个参数太小，可能导致过滤过多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深度。建议根据物种情况进行修改，如1万，统计更准确。这个参数太小，可能导致过滤过多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比例为杂合度。.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.isto：频数分布直方表，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评估基因组特征

详细介绍电流互感器功能区别3CT SR ZCT

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

详细介绍电流互感器功能区别3CT SR ZCT

最近热搜

历史搜索 清除历史记录

历史搜索清除历史记录