用scaffold10x基于10x genomics数据做gen

背景介绍

在得到初步的组装结果之后，如果手上有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会发现它做的工作是去自动下载bwa、smalt和pigz这几个依赖软件，可以直接用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. 最近不知道选什么图片作为头图比较好，于是突发奇想就用自己拍的照片做头图好了~这样也减少了使用有版权照片的法律/商业风险。昨天的头图是我养的小兔子图图，今天的图是前段时间热气球节上拍的照片，希望你喜欢。