为什么记录这个过程是因为我发现我老了，许多之前做过的分析，写过的脚本，再看一眼都回忆不起来了，遇到同样的问题又要重新想一遍，解决完就又忘掉了。所以还是多记录一下，帮助自己回忆
问题来源：
hifiasm 的HIC 模式可以分出两套hap1和hap2，而且质量较高。我用同一套参数组装挂载了两个物种，其中杂合度高的物种分单倍型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
网友评论