一文搞懂ALLHiC构建Allele.ctg.table

        本人接触生信4个多月了，发现与机器学习不同的是生信涉及的一些知识尤其是一些坑点很少有资料可以查询到，所以我决定写一些文章去记录一下平时遇到的一些问题，一来可以作为笔记供日后查询，二来可以给一下有需要的人提供一些思路，注意本人为生信方面的萌新宝宝，所写不一定全部正确，欢迎大家指正！！！！

        好接下来进入正题，ALLHiC作为一款2019年发布的scaffolding工具，网上很多人在教大家如何使用这款工具去组装自己的基因组，但是缺很少有人将构建Allele.ctg.table说清楚，接下来我将从下载近缘物种的cds序列以及gff3文件讲起，将大家一步步构建自己的Allele.ctg.table并避免其中的坑点。

        目前官方提供了俩种方法去构建Allele.ctg.table，一种是基于BLAST，另外一种是基于GMAP。基于BLAST的方法较为繁琐且需要的文件较多，这里不推荐。我平时使用的也是基于GMAP的方法。

qry=waxapple.asm.bp.p_utg.gfa.fa  # 你的多倍体序列
ref_cds=cds_from_grande_genomic.fna # 同源物种的cds

N=4 #你的基因组倍性, 4表示的是4倍体
/data/liyanchun/gmap-2023-12-01/util/gmap_build -D . -d DB $qry -t 16
/data/liyanchun/gmap-2023-12-01/src/gmap -D . -d DB -t 12 -f 2 -n $N  $ref_cds > gmap.gff3

perl /data/liyanchun/ALLHiC/scripts/gmap2AlleleTable.pl  oleosum_genomic.gff3

其中qry为自己的多倍体序列，ref_cds为同源物种的cds序列，N为你自己基因组的倍性，几倍体就填几，接下俩俩行代码需要修改的就是 -t 这是指代的是你的线程数，最后一行代码代码这里就有一个小小的坑点，这里的oleosum_genomic.gff3为你的近缘物种的gff文件，需将其修改为自己下载好的文件。还有一个坑点是倒数第二行代码中最后输出文件的名字只能是 gmap.gff3！！！！！不然代码会报错，报错的原因为最后一步中gmap2AlleleTable.pl 指定的输入文件名称为 gmap.gff3

        那么，这里有俩个我们需要且没有的文件，就说近缘物种的cds序列和gff文件，那怎么知道你自己的基因组的近缘物种有什么呢？一个简单且有效的方法GPT。以我的莲雾为例

然后就去NCBI上面搜索（https://www.ncbi.nlm.nih.gov/genome/）

点击红框即可，然后你会发现进去之后不一定会有提供对应的gff和cds文件，那么咱们就得换一种思路了，还是以莲雾为例，这几个物种的拉丁文前面一摸一样，那么我们直接搜索Syzygium，

经查找发现Syzygium grande有提供对应的文件，那么我们点进去，再往下滑找到Download并点击.

将cds序列和gff文件勾选，并再次下滑，找到Download并点击即可。
        那么我们现在就拥有了需要的全部文件，这是否意味着我们可以构建自己的table表了，理论上是这样的但是，注意这里出现了最大的坑点。我们直接上代码

#!/usr/bin/perl -w

die "Usage: perl $0 ref.gff3\n" if(!defined ($ARGV[0]));
my $refGFF = $ARGV[0];
open(IN, "grep 'gene' gmap.gff3 |") or die"";#可修改（1）
while(<IN>){
	chomp;
	my @data = split(/\s+/,$_);
	my $gene = $1 if(/Name=([^;\n]*)/);#可修改（2）
	#print $gene;
	$infordb{$gene} .= $data[0]."	";
	}
close IN;

open(OUT, "> Allele.ctg.table") or die"";
open(IN, "awk '\$3==\"gene\"' $refGFF |") or die"";#可修改（3）
while(<IN>){
	chomp;
	my @data = split(/\s+/,$_);
	my $gene = $1 if(/Name=(\S+)/);#可修改（4）
	   $gene =~ s/;.*//g;
	next if(!exists($infordb{$gene}));
	my @tdb = split(/\s+/,$infordb{$gene});
	my %tmpdb = ();
	map {$tmpdb{$_}++} @tdb;
	print OUT "$data[0]	$data[3]	";
	map {print OUT "$_	"} keys %tmpdb;
	print OUT "\n";
	}
close IN;
close OUT;

这是gmap2AlleleTable.pl里的代码，这也就解释了为什么前面我们说使用gmap时输出文件的名字必须是gmap.gff3文件。这段代码什么意思呢，一句话概括就是用gmap.gff3中代表gene这一行中Name标识符对应的内容去匹配近缘物种中gff3文件中代表gene这一行中Name标识符对应的内容从而构建table表。我帮大家说：不懂给我绕迷糊了QAQ ，那么好 我们直接上图

这是gmap.gff3中的内容，从这我们就能很什么叫做gene这一行中Name标识符对应的内容。大家也可以发现这里不止有代表gene的一行还有代表着mRNA和CDS的行，那么这有什么用呢，别急，接着往下看。

这是下载好的莲雾的近缘物种的gff表，大家还记得我们之前提到的gmap2AlleleTable.pl匹配原则吗？但是这里发现gene这行中Name对应的内容与gmap.gff3中代表gene这行Name对应的内容毫无相关，反正是CDS这行的Name有点相关，但也不完全一样，多了很多东西。怎么办？老规矩转化思路。多了给他删了就行。我们直接上代码。

sed -E 's/Name=lcl\|[A-Za-z0-9]+\.[0-9]+_cds_([A-Za-z0-9]+_[0-9]+\.[0-9]+)_([0-9]+)[^;]*/Name=\1/' gmap.gff3 > gene_g.gff3

这段代码什么意思呢？一句话就是用正则表达式去匹配Name=lcl|JAMXDD010000100.1_cds_KAI6668043.1_3;中的KAI6668043.1，并将其修改为Name=KAI6668043.1;，那么好，具体怎么匹配的呢？大家如果换成了自己的数据集又改如何修改呢？老规矩GPT，我也不是预言家，也不知道大家到时候的Name格式是啥样的，马上都2025年了大家一定要学会使用GPT，你不学会使用他你就会被他取代。这里给大家提供一个我匹配的思路，我是根据下划线的数量去匹配的。

修改好后的文件就变成了这样。但是我们上面对Name就行删减时将gmap.gff3 修改为gene_g.gff3怎么办？而且这里修改的是CDS对应的Name也不是gene对应的名字还是匹配不上呀。那么接下来我们就对gmap2AlleleTable.pl中的代码就行修改。

一共就是这五处地方可以修改，2这里我们就可以将其修改为自己用gmap构建的gff文件最终的名字，在此例中我们需要将4这里修改为CDS即可。那么当大家自己建table时遇到其他的匹配原则时再具体进行修改即可，如这里用的时Name进行的匹配，若大家匹配时发现使用ID进行匹配更加简单，将3，5改为ID即可。好了到这里为止，所有的坑点都已讲完需要本文对大家有所帮助，最后如果本文中有错误欢迎大家指正！！！

参考：

ALLHiC续: 如何构建Allele.ctg.table_allhic allele.ctg.table-CSDN博客

GMAP一款比对工具用于ALLHiC构建等位基因表_gmap比对用法-CSDN博客