卖萌哥

V1

2022/10/17阅读:15主题:橙心

用scaffold10x基于10x genomics数据做genome scaffolding

背景介绍

在得到初步的组装结果之后,如果手上有10x genomics的基因组测序数据的话,除了可以用supernova基于10x数据独立组装出一个新版本的基因组外(详见我的往期推送10x基因组数据的组装),还可以借助10x的数据把contig序列给连成更长的片段,做一个初步的scaffolding。

scaff10x的开发者来自于大名鼎鼎的桑格研究所的高性能算法团队(High Performance Algorithms Group),下面是它的github地址: https://github.com/wtsi-hpag/Scaff10X

软件原理

scaff10x的工作原理是:

  1. 把barcoded从10x的原始数据中提取出来,放到序列的id行里以备后续使用;
  2. 把10x的序列比对到基因组上,这里可以选择用bwa或者SMALT
  3. 把barcode根据contig和比对坐标(mapping coordinates)给排个序;
  4. 建立一个关系矩阵(relation matrix)用以记录潜在的可相互连接的contig之间的共享barcodes信息;
  5. 在找到最接近的contigs后把他们根据顺序和方向连接起来。

软件安装

git clone  https://github.com/wtsi-hpag/Scaff10X.git
cd Scaff10X
./install.sh

打开这个install.sh会发现它做的工作是去自动下载bwasmaltpigz这几个依赖软件,可以直接用conda安装。我这里为了方便还是让他默认安装吧。

一个小bug是由于pigz从2.6版本更新到了2.7,因此需要手动修改一下install.sh里的pigz的版本,从2.6修改到2.7,否则会报错的哦。

安装好后记得把软件加入到环境变量中,当然也可以写绝对路径调用。

软件运行

scaff10x \
-nodes 120 \ # 设置运行的线程数
-size 2.0 \ # 基因组的大致大小,单位是Gb,可以写0.5, 1.0, 2.0 (Gb)
-longread 1 \ # 基因组是用什么组装的?1代表三代数据,0代表二代数据。
-gap 100 \ # 设置gap的大小,默认是100
-matrix 5000 \ # 设置relation matrix的大小,默认是2000
-reads 10 \ # 上面原理中第一步和第二步的最小共享barcode的reads数目,默认是10
-link 8 \ # 上面原理中第一步和第二步的最小的被共享的barcode的数目,默认是8
-score 20 \ # 最小的平均比对质量,默认是20
-edge 50000 \ # scaffolding时边界的长度,默认是50000
-block 10000 \ # 决定最接近的相邻者的长度。默认是50000
-plot hap2_length.png \ # 打印出barcode的长度分布。
/path/to/test.hic.hap2.p_ctg.fasta \ # 用于scaffold的contig,即前期组装结果
/path/to/test_L001_R1_001.fastq.gz \ # 10x数据reads 1
/path/to/test_L001_R2_001.fastq.gz \ # 10x数据reads 2
test.hap2.scaff10x_block10000.fasta # 最终结果。

其实看起来设置了很多的数据,其实大多数都是默认值。只是修改了block的数值。

最终效果还不错,我的数据从2155条contig减少到了1037条,N50也从3.5 Mb提升到了37 Mb。直接翻了十倍。当然,我这里用的是hifiasm的单倍体的数据,本身N50较短。

萌哥碎碎念

  1. 就我的观察而言,10x scaffolding对于较长的contig的贡献比较一般,但是一些较短的contig确实有明显的提高,这也非常符合预期。因为10x的数据是基于illumina平台的二代短序列,本身较短,即使有共享的barcode辅助延伸,对于长片段的作用也非常有限。
  2. 最近不知道选什么图片作为头图比较好,于是突发奇想就用自己拍的照片做头图好了~这样也减少了使用有版权照片的法律/商业风险。昨天的头图是我养的小兔子图图,今天的图是前段时间热气球节上拍的照片,希望你喜欢。

分类:

工具介绍

标签:

开源软件

作者介绍

卖萌哥
V1