图泛基因组的组装

什么是图泛基因组#

因为真核生物的基因组包含大量的间隔区域, 同时存在较大的个体差异, 因此泛基因组可以表示为该物种内的所有DNA序列集合, 因此泛基因组组装的核心包含一个核心基因组和表示其他个体的可选基因组.
图泛基因组和以往的线性泛基因组组装不同, 线性泛基因组只是将所有的可选基因组中识别的PAVs组装为线性基因组, 但是该方法无法描述个体间变异来源或者对基因组中存在的PAVs进行精准定位; 图泛基因组则将所有线性基因组之间的变异以图形的格式储存.

图泛基因组的组装#

图形泛基因组的组装主要依赖minigraph-cactus进行:

准备基因组配置文件#

需要准备一个配置文件包含所有的基因组信息(seqFile), 格式为:

1
GenomeA /path/to/genomeA.fa
2
GenomeB /path/to/genomeB.fa
3
GenomeC /path/to/genomeC.fa

使用cactus进行graph pangenome的组装:#

在该步骤中使用了singularity 已经cactus的镜像构建以避免源代码构建cactus产生的错误:
singulairty pull cactis.sif docker://quay.io/comparative-genomics-toolkit/cactus:latest
cactus是一个组装泛基因组的复杂工具, 具有大量的可调参数与可选软件.
最简单的运行方式为: cactus-pangenome ./jobStore ./seqFile --reference GenomeA --maxCores 32 --outDir ./ --outName Out

NOTE

解读cactus 的参数设置:
jobStore: cactus的运行结果存放位置.
seqFile: 为步骤1 准备的基因组配置文件, 默认的图泛组装核心(—reference)为配置文件中的第一个基因组.
—maxCores (int): 运行一个job时的最大核心数.
—outDir: 运行结果存放的目录.
—outName: 运行结果存放的文件.
—gbz —gfa —vcf —hal: 结果输出的格式, 可以多选 可选参数
cactus具有丰富的可调参数, 以下给出一些我自己使用的时候觉得可以调整, 以提高检测效率的参数.
—reference: 默认图泛基因组核心, 这个不用过多解释了.
—minLen default=100k: 只有大于该值的片段才会被识别为该图需要的结构变异(SV), 如果构建图泛的物种亲缘关系接近可以将该值调小.
—filter default=2: 过滤掉出现在少于--filter个单倍型中的节点, 如果需要所有节点可以设置为0.
—clip default=1000: 剪切掉未对齐且长度小于--clip值的末端序列.

使用vg进行自动索引#

由于后续的比对与计算SNP的位点等功能需要使用到vgtools, 而比对的核心是使用vg giraffe, 因此需要先自动索引使得生成min文件等.
在该步骤中, vg会读取由cactus组装得到的gfa文件, 在读取该步骤前, 我们需要先使用--workflow giraffe 确定vg的工作流是使用giraffe 进行后续reads的比对, 因此他会自动运行一个pipeline, 包含construct, index, gbwt(如果有的话), snarls, index, 分别表示: 构建图, 构建图的索引, 构建单倍体索引(如果有的话), 构建气泡图, 寻找基于气泡图(snarls)的索引距离.

NOTE

.snarls是文件拓扑索引, 其中记录了所有参与构建基因组的哪些部分是等位基因变异或剪接分支.

如果有多个群体和多个部位的RNA-seq部位, 还可以进行RNA-seq的比对, 用于观察到底是什么基因变异, 而RNA-seq的比对需要使用vg index中计算的gcsa文件, 因此在设计的pipeline中可以重新跑一次vg autoindex, 但是workflow改为mpmap以生成gcsa文件, 否则无法进行图转录组比对.

NOTE

gcsa文件是什么?
gcsa是一种全文索引, 可以用于处理graph的结构, 包含的信息为: K-mer路径信息, 记录了图中所有长度达256bp的可能路径.

图泛基因组的注释#

图泛基因组与DNA测序数据#

由于图泛基因组以graph的方式储存了多个基因组相对于reference基因组的差异, 因此也可以使用多个群体(物种)的全DNA测序数据去mapping到图泛基因组上, 以展示DNA测序层面上, 测序的群体与图泛基因组所有包含的基因组之间的结构差异关系, 并且总结得到不同群体中所包含的变异关系.

1. DNA测序数据mapping到图泛基因组
   由于已经先使用vg autoindex提前构建好了基于giraffe所需的索引, 因此可以直接运行vg giraffe, 并且输入测序得到的fq文件即可.

   1
   vg giraffe -Z graph.gbz [-d graph.dist [-m graph.withzip.min -z graph.zipcodes]] <input options> [other options] > output.gam
   2
   # .dist, .min, .zipcodes 都是通过autoindex可以自动生成的
   

由于gam文件很大, 如果是输入的是群体测序的数据, 那么会得到非常多个gam文件, 可能会对硬盘储存的压力很大, 因此可以使用vg pack压缩gam文件.

2. 将比对得到的图数据转换为线性坐标
   通过图泛基因组比对得到的变异数据(.gam)事实上是记录了reads比对上图泛基因组的详细信息. 其中储存了图中每个node, 每个edge被reads覆盖的信息, 其中存在了极多的冗余的信息, 因此需要使用vg call, 因此我们可以判断在图中某一个分支处，Read 的覆盖程度是否足以证明该变异在样本中真实存在, 亦或者得到图泛基因组中的结构差异。但是由于其是基于泛基因组的路径比对得到的, 所以他不能基于未知的群体去判断未知群体中拥有的特异序列, 此时可能需要新的方法 — vg augment