记录下如何手动解决hifiasm 装单倍型失败的例子， 许多之前做过的分析，写过的脚本，没整理都有点模糊了，遇到同样的问题又要重新想一遍。所以还是对典型问题多记录一下，帮助自己回忆

问题来源：

hifiasm 的HIC 模式可以分出两套hap1和hap2，而且质量较高，当然如果没有父母本的二代测序，仅用HIC mode分单倍型，在某些区域是有可能由于物种杂合度过低而分不出来的，我用同一套参数组装挂载了两个物种，其中杂合度高的物种分单倍型hap的效果非常好，但另一个杂合度低的物种在调图时遇到了这样的情况

我是将两套hap cat到一起挂载的，这样可以识别出一些潜在的组装错误。 在这两张hic图中，可以发现有些contig的部分片段和两个单倍型都有着很强的互作，在juicer调图的时候，我把他当作污染或者杂质将它从scaffold上切下来放到contig里面了。

e93ca7cb744ed8ace0baf68e4b4e1c6.png

a1c8ee0f7b7c96101229b28d72affc4.png

在挂载结束后，我将这两个物种的两个hap与已发表的近缘物种做共线性发现这是个组装错误，是hifiasm在面对杂合度很低的片段时，将本应分成两个单倍型的片段仅仅组装出了一个单倍型，这样在挂载的时候，这个片段和两套单倍型都有很强的互作，下图所示

image.png

问题

因此这里应该把这个剪下来的片段*2，手动插入到这两个scaffold中间。

实现方法

提取这三个scaffold， 将他们根据500个N拆成contig，做比对，看看这个嵌合的片段应该插入到哪里，代码如下

#!/usr/bin/envs perl -w

use strict;

my $usage=<<USAGE;
###usage:

perl $0 target_scaffold_name_list.txt scaffolding.fasta > result.fa


###Explanation:

target_scaffold_name_list.txt:
---Start---
HiC_scaffold_19
HiC_scaffold_20
---End---

scaffold.fasta
---Start---
>HiC_scaffold_1
...
...
>HiC_scaffold_19
...
>HiC_scaffold_20
...
---End---


USAGE

die "$usage\n" unless @ARGV == 2;

my $list=shift;
my $fa=shift;

my %name;

open IN, "$list" or die "";

while (<IN>){
        chomp;
        $name{$_}=1;
}

close IN;

my ($seqID, %seq);
open FA, "$fa" or die "";

while (<FA>){
        chomp;
        if (m/^>(\S+)/){$seqID=$1;}
        else{$seq{$seqID}.=$_;}
}

close FA;

foreach my $tar_ID (sort keys %name){
        die "target_id can't be found in scaffolding.fa, please check Names of target_list\n" unless my $seq=$seq{$tar_ID};
        my $N=("N") x 500;
        $seq=~ s/$N/1/g;
        my @fragment=split /1/, $seq;
        my $num=1;
        foreach my $contig (@fragment){
                next if $contig eq "";
                print ">${tar_ID}_${num}\n$contig\n";
                $num++;
        }
}

第一个输入文件为所需提取的scaffold名称如下图，第二个为总的scaffold.fasta，结果文件是如下图

image.png

image.png

再做共线性

Ptri_Chr06-21_22_1202.PNG

观察共线性图可知，我们可以将这些contig 重新按照正确的顺序用N连接，需要反向的contig则在其名字后面标注":R"即可，文件如下

image.png

代码如下

#!/usr/bin/envs perl -w

use strict;

my $usage=<<USAGE;

#usage:

perl $0 contig_name_for_scaffolding.list total_contig.fa > after_joint.fa

#explanation

contig_name_for_scaffolding.list
---start---
HiC_scaffold_21_1 HiC_scaffold_21_2
HiC_scaffold_22_2 HiC_scaffold_22_4 HiC_scaffold_1202_7
---end---


USAGE

die "$usage\n" unless @ARGV == 2;

my %hash;
my $ad=\%hash;
my @name;

open ID, "$ARGV[0]" or die "";

while (<ID>){
        chomp;
        my @line = split/\s+/, $_;
        my $name= shift @line;
        push @name, $name;
        push @{$ad -> {$name}}, @line;


}

close ID;


my ($seqID, %seq);
open FA, "$ARGV[1]" or die "";

while (<FA>){
        chomp;
        if (m/>(\S+)/){$seqID=$1;}
        else{$seq{$seqID}.=$_;}
}

close FA;

my $N=("N") x 500;
#$seq=~ s/$N/1/g;

foreach my $id (keys %hash){
        my @line = @{$hash{$id}};
        my @seq;
        foreach my $contig ( @line ){
                if ($contig=~ m/:R$/){
                        $contig =~s/:R//;
                        my $seq = $seq{$contig};
                        $seq=~ tr/ATCGagtc/TAGCtcag/;
                        $seq=reverse ($seq);
                        push @seq, $seq;
                }else{
                        my $seq = $seq{$contig};
                        push @seq, $seq;
                }
        }
        my $scaffold=join "$N", @seq;
        print ">$id\n$scaffold\n";
}

再做共线性：问题解决。这里我们注意到HiC_scaffold_22中间是有一块空白的，这是因为我用的是last比对软件，由于我们是直接将一个片段复制成两个片段，因此这个片段被last软件过滤掉了，所以图中没有任何共线性，但其实他和HiC_scaffold_21的片段是完全一样的。这里用last的话，一定要把过滤1e-5去掉，不然我们复制进去的两个片段是没有任何映射的。

image.png

问题解决

hifiasm这个软件速度快，效果好，还能自动装出两套单倍型，是现在hifi数据的热门选择之一，但这个模式也不是完美的，他在某些纯合度片段也会分不开两套单倍型，当然如果有父本母本的二代reads，用hifiasm的另一个模式的话效果应该比hic 模式好很多。这里主要记录了怎么手动将这种大片段的嵌合错误纠正过来。