1.介绍
Merqury是一款基于k-mer数据集评估无参考基因组的方法。对于trios，Merqury还可以评估特定单倍型的准确性、完整性、相位块连续性以及相位切换错位(switch errors,一个单倍型上的某个碱基被错误组装到另一个单倍型上)。
2.安装
conda create -n merqury -c conda-forge -c bioconda merqury openjdk=11
conda activate merqury
Rscript $MERQURY/plot/plot_spectra_cn.R --help
3.输入文件
ln -s $MERQURY/merqury.sh
###Tips1：Merqury默认所有meryl数据库均为.meryl结尾
###Tips2：<> 必须输入文件 [] 可选输入文件
Usage: merqury.sh <read-db.meryl> [<mat.meryl> <pat.meryl>] <asm1.fasta> [asm2.fasta] <out>
<read-db.meryl> : k-mer counts of the read set
<mat.meryl> : k-mer counts of the maternal haplotype (ex. mat.hapmer.meryl)
<pat.meryl> : k-mer counts of the paternal haplotype (ex. pat.hapmer.meryl)
<asm1.fasta> : Assembly fasta file (ex. pri.fasta, hap1.fasta or maternal.fasta)
[asm2.fasta] : Additional fasta file (ex. alt.fasta, hap2.fasta or paternal.fasta)
*asm1.meryl and asm2.meryl will be generated. Avoid using the same names as the hap-mer dbs
<out> : Output prefix
Arang Rhie, 2020-01-29. arrhie@gmail.com
4.构建数据库
4.1 计算合适的k-mer值
###当不知道使用什么k-mer值合适时，运行下面命令
###人无论是haploid(3.1G)还是diploid(6.2G)，k-mer均为21
sh $MERQURY/best_k.sh <genome_size>
4.2 单个reads建库
meryl k=$k count *.fastq.gz output $genome.meryl
4.3 多个reads建库
###有多少个reads列表，就有多少个库
meryl k=$k count output read$i.meryl read$i.fastq.gz
###将上面的库合并为一个库，union-sum取的是各个数据库的并集
meryl union-sum output $genome.meryl read*.meryl
###作者设置了很多参数便于筛选k-mer库，你可以自定义使用哪些库和哪些k-mer
less-than N return kmers that occur fewer than N times in the input. accepts exactly one input.
greater-than N return kmers that occur more than N times in the input. accepts exactly one input.
equal-to N return kmers that occur exactly N times in the input. accepts exactly one input.
not-equal-to N return kmers that do not occur exactly N times in the input. accepts exactly one input.
increase X add X to the count of each kmer.
decrease X subtract X from the count of each kmer.
multiply X multiply the count of each kmer by X.
divide X divide the count of each kmer by X.
divide-round X divide the count of each kmer by X and round results. count < X will become 1.
modulo X set the count of each kmer to the remainder of the count divided by X.
union return kmers that occur in any input, set the count to the number of inputs with this kmer.
union-min return kmers that occur in any input, set the count to the minimum count
union-max return kmers that occur in any input, set the count to the maximum count
union-sum return kmers that occur in any input, set the count to the sum of the counts
intersect return kmers that occur in all inputs, set the count to the count in the first input.
intersect-min return kmers that occur in all inputs, set the count to the minimum count.
intersect-max return kmers that occur in all inputs, set the count to the maximum count.
intersect-sum return kmers that occur in all inputs, set the count to the sum of the counts.
subtract return kmers that occur in the first input, subtracting counts from the other inputs
difference return kmers that occur in the first input, but none of the other inputs
symmetric-difference return kmers that occur in exactly one input
###举个例子，你想选出在db1和db2中均出现，且该k-mer在db2中还是unique的，将筛选出来的k-mer保存到db3
meryl output db3.meryl intersect db1.meryl [ equal-to 1 db2.meryl ]
4.4 构建hap-mer dbs
###本地运行-选一种运行就行
sh $MERQURY/trio/hapmers.sh maternal.meryl paternal.meryl child.meryl
###提交任务-选一种运行就行
sh $MERQURY/_submit_hapmers.sh maternal.meryl paternal.meryl child.meryl
###输出结果
* parental specific dbs: `mat.only.meryl` and `pat.only.meryl`
* inherited dbs: `mat.inherited.meryl` and `pat.inherited.meryl`
* inherited hap-mer dbs (which will be used for evaluation): `mat.hapmers.meryl` and `pat.hapmers.meryl`
* inherited_hapmers.png: k-mer distribution of the inherited dbs and cutoffs used to generate hap-mer dbs
5.评估基因组
5.1 一个基因组(混合单倍型)
# I don't have the hap-mers
$MERQURY/merqury.sh read-db.meryl asm1.fasta out_prefix
# Using the example above
$MERQURY/merqury.sh F1.k18.meryl athal_COL.fasta test-1
# I have the hap-mers
$MERQURY/merqury.sh read-db.meryl mat.meryl pat.meryl asm1.fasta out_prefix
# Using the example above
$MERQURY/merqury.sh F1.k18.meryl col0.hapmer.meryl cvi0.hapmer.meryl athal_COL.fasta test-1
5.2 两个基因组(分型单倍型)
# I don't have the hap-mers
$MERQURY/merqury.sh read-db.meryl asm1.fasta asm2.fasta out_prefix
# Using the example above
$MERQURY/merqury.sh F1.k18.meryl athal_COL.fasta athal_CVI.fasta test-2
# I have the hap-mers
$MERQURY/merqury.sh read-db.meryl mat.meryl pat.meryl asm1.fasta asm2.fasta out_prefix
# Using the example above
$MERQURY/merqury.sh F1.k18.meryl col0.hapmer.meryl cvi0.hapmer.meryl athal_COL.fasta athal_CVI.fasta test-2
6.结果结果
6.1 QV文件
col 682359 123511626 35.1183 0.000307729
cvi 506306 122349407 36.3761 0.00023035
Both 1188665 245861033 35.6991 0.00026921
第一列是两个分型单倍型的名字；第二列为仅存在于组装结果中的k-mer；第三列为组装结果中的所有k-mer；第四列为QV值；第五列为错误率
Tips：QV方法描述用于理解
image.png
6.2 完整度文件
#完整度计算方法：第三列/第四列
#
col all 104975080 125303808 83.7764
cvi all 104809523 125303808 83.6443
both all 123134729 125303808 98.2689
第一列为两个分型单倍型的名字；第二列表示哪些k-mer用于计算完整度；第三列为组装中存在的k-mer；第四列为reads数据库中存在的k-mer；第五列为完整度
6.3 其它文件
asm_only.bed和asm_only.wig两个文件一般用于追踪错误，
### 主要是看个深度，看看有没有位点的深度比较低的
igvtools count ${asm}_only.bed ${asm}_only.tdf $ref.fai
7. 如何提取仅存在于reads-only中的reads
Question: a way to extract read-only k-mers · Issue #9 · marbl/merqury
8.参考
GitHub - marbl/merqury: k-mer based assembly evaluation
https://zerobio.github.io/archives/1959810678.html
