基因组survey

物种的选择需要有研究价值，无论是经济价值还是其他的研究价值。二倍体、低杂合、低重复物种相对好组装。但是高杂合有利于进行单倍型分型。大基因组意味着更高的测序成本，以及更大的计算资源。

通过NCBI查询待测物种已发表的各类数据，比如是否有参考基因组、收录了多少蛋白序列、是否有转录组数据等。

查询PubMed等文献数据库，了解待测物种的研究情况。

C-value是单倍体细胞DNA含量或二倍体细胞中的DNA含量的一半，以重量计算，单位通常是皮克（pg），1pg约为978 Mb。部分植物物种可以通过 https://cvalues.science.kew.org/ 查询，动物可以通过 http://www.genomesize.com/ 查询。

先通过文献查找是否有测序物种的核型信息，如果没有可以通过实验获取。

基本原理是DNA含量越高，荧光强度越大。

对于没有参考基因组的物种需要先测低深度二代数据对基因组进行评估，一般基于流式评估测30-50X，包括：

评估基因组大小，决定测序量
重复序列比例，现在基本都用HIFI或者超长，一般问题不大
杂合度，对于二倍体或者多倍体，每套染色体会有差异
GC含量，过高或过低GC含量的区域二代测序可能难以被测到，现在基本都用HIFI或者超长，一般问题不大

之后一般按k=17进行kmer分析，注意一般选择奇数长度，也可以使用k=19，21，23等，k越大对内存的需求更大。可以几个值都分别进行评估，然后综合比较。一般kmer变大，主峰左移。

评估的前提假设是，基因组的测序是随机的、均匀的。碱基深度和kmer深度不一样，可以用换算公式计算。 Kmer出现的频数服从泊松分布。根据曲线获得Kmer深度估计值c，估计基因组大小。Kmer在depth=1时认为是错误情况，计算错误率，可以用于修正基因组大小。

二倍体中杂合度为1%时，杂合峰和纯合峰可能就基本持平了。杂合度变大时，杂合峰上升、纯合峰下降更剧烈。所以要分别以各峰作为主峰进行评估，然后看那个评估结果与流式细胞的结果更接近。

有污染的数据不适合做评估。

进行基因组 survey 分析的测序数据，需要满足高质量要求，因此一般使用二代测序数据进行，不建议使用 ONT 及 Pacbio CLR 等测序错误率高的数据进行 kmer 分析。。错误率高的数据不适合做评估，比如ONT。错误率变大，主峰左移，错误率达到4%时，kmer频率分布曲线偏离实际曲线已经十分严重了。所以用illumina数据进行评估，三代的Pacbio HIFI数据也可以。Pacbio CLR数据自矫正后的数据也可以用于评估。

可以使用 fastqc对数据质量进行检查，然后用fastp过滤。fastqc和fastp的安装及使用请在http://www.genomes.cn/主页搜索框检索“fastqc”和“fastp”。

项目地址：https://github.com/gmarcais/Jellyfish

安装

conda install -c conda-forge -c bioconda jellyfish

pre=Kmer_21

# 生成reads 解压脚本
ls R*.fastq.gz | \
awk '{print "gzip -dc "$0 }' > generate.file


# 进行kmer计算
jellyfish count \
-t 8 \ # 线程
-C \ # canonical k-mers
-m 21 \ # kmer大小
-s 1G \ # 哈希大小
-g generate.file \ # 数据
-G 2 \ # 数据生成并行任务数
-o $pre # 输出文件

# 统计kmer总数
jellyfish stats -o $pre.stat $pre

# 生成kmer频数统计表
jellyfish histo \
-v \
-t 8 \
-h 10000000 \ # 输出深度阈值，大于stat中的 Max_count 值
-o $pre.histo \ # 输出文件 
$pre # 输入文件

pre=HiFi_21mer

# 生成reads 解压脚本
ls hifi.fastq.gz | \
awk '{print "gzip -dc "$0 }' > generate.hifi

# 进行kmer计算
jellyfish count \
-t 8 -C \
-m 21 \
-s 1G \
-g generate.hifi -G 2 -o $pre

# 统计kmer总数
jellyfish stats -o $pre.stat $pre

# 生成kmer频数统计表
jellyfish histo -v -t 4 -h 10000000 -o $pre.histo $pre

GCE对于分析原理有较详细的说明。该软件适用于单倍体和二倍体的物种，并且GCE只适用于杂合峰低于纯合峰的情况。

gce \
-f Kmer_21.histo \ # 输入，kmer频数表
-c 72 \ # 期望深度估计值
-H 0\ # 纯合模式
-g 353761430 \ # 总kmer个数，在stat文件查看
-M 500 >gce.table 2>gce.log # 频率最大对应横坐标的3倍即可

# 杂合模式运行
gce \
-f Kmer_21.histo \ 
-c 72 \
-H 1 \ #杂合模式 
-g 353761430 \
-M 500 >gce.table 2>gce.log

结果存储在log文件

如果是二倍体，推荐findGSE，findGSE对于低测序深度（10X-20X）、高测序错误率（5%）、高杂合（2%以上）都有较高的容忍度。该软件适用于单倍体和二倍体的物种

library("  ")

# 自动获取kmer cov，纯合模式 
findGSE(histo="Kmer_21.histo", sizek=21, outdir="Kmer_21.histo.findGSEhomo")

# 输入kmer cov, 杂合模式
findGSE(histo="Kmer_21.histo", sizek=21, exp_hom = 140, outdir="Kmer_21.findGSEhet)

# 运行
Rscript step2.run_findGSE.R

如果是三倍体到六倍体，推荐GenomeScope2，该软件适用于单倍体到六倍体的物种

在线版：http://qb.cshl.edu/genomescope/genomescope2.0/

genomescope.R \ 
-i Kmer_21.histo \ # 输入文件
-o genomescope_out \ # 输出前缀
-p 1 \ # 倍型
-k 21 \ # kmer大小
-m 10000 # kmer最大深度

2023-10-02

基因组survey

公共数据库查询

核型信息

流式细胞仪估计基因组大小

低深度的二代数据

背景介绍

使用fastqc、fastp进行二代测序数据评估、过滤

使用jellyfish进行Kmer计数

基于二代测序数据运行

基于 HiFi 数据运行

对计数结果进行分析

GCE

findGSE（推荐）

GenomeScope2

本页目录