ALLHiC续: 如何构建Allele.ctg.table

Allele.ctg.table是ALLHiC用于过滤多倍体基因组中因等位序列相似引起的HiC噪音的必要输入。

构建的方法有很多种，这里列举两个方法。

方法一: 基于BLAST

基于BLAST的方法需要你提供一个染色体级别组装的近缘物种。需要提供4个文件，近缘物种的CDS和GFF文件，多倍体物种的CDS和GFF文件。最重要的是，你得保证CDS序列中的编号和GFF文件的基因编号一致，这也是后续处理至关重要的一步。

问题一: 我们如何快速获取待组装物种的CDS和GFF文件？最快速的方法就是利用转录组序列基于STAR StringTie的组合进行拼接，并利用gffread从中提取出cdna序列。

假设你的多倍体文件名为contig.fa, 转录组数据为read_R1.fq.gz, read_R2.fq.gz

mkdir -p STAR
ref=contig.fa
THREADS=40
fq1=read_R1.fq.gz
fq2=read_R2.fq.gz
prefix=sample
# 建立索引
STAR \--runThreadN $THREADS \--runMode genomeGenerate \--genomeDir STAR \--genomeFastaFiles $ref
#比对
STAR \--genomeDir STAR \--runThreadN $THREADS \--readFilesIn $fq1 $fq2 \--readFilesCommand zcat \--outFileNamePrefix $prefix\--outSAMtype BAM SortedByCoordinate \--outBAMsortingThreadN $THREADS \--outSAMstrandField intronMotif \--outFilterIntronMotifs RemoveNoncanonical
# 拼接转录本
stringtie ${prefix}Aligned.sortedByCoord.out.bam -p $THREADS -o qry.gtf
# GTF 转成 GFF3
gffread qry.gtf -o qry.gff
# 提取cdna
gffread qry.gff -g $ref -w qry.fa

问题二: 后续的处理脚本要求输入序列的编号和GFF的基因编号相同，并且编号在GFF文件中的格式为NAME=xxx。

对于我们上一步处理的qry.gff和qry.fa，我们只需要简单粗暴地替换qry.gff里的信息即可。

sed -e 's/transcript/gene/' -e 's/ID/Name/' qry.gff > qry_gene.gff

对于近缘物种的cds序列和gff文件，考虑到后续的脚本只认gene所在的行，而cds里面存放的是转录本序列，我们可以提取mRNA所在行，然后将其改名成gene。假设参考序列是ref.fa, 参考的GFF文件名是ref.gff

ref=ref.gff
grep '[[:blank:]]mRNA[[:blank:]]' $ref | sed -e 's/mRNA/gene/' -e 's/ID/Name/' > ref_gene.gff

综上你得到了如下四个文件, ref.fa, refgene.gff, qry.fa, qrygene.gff

下面的操作，参考自官方教程完成，

第一步：建立BLAST索引

makeblastdb -in ref.fa -dbtype nucl

第二步: 将我们物种序列比对到近缘物种

blastn -query qry.fa -db ref.fa -out qry_vs_ref.blast.out -evalue 0.001 -outfmt 6 -num_threads 4 -num_alignments 1

第三步: 过滤低质量的HIT，此处标准是相似度低于60%, 覆盖度低于80%

blastn_parse.pl -i qry_vs_ref.blast.out -o Eblast.out -q qry.fa -b 1 -c 0.6 -d 0.8

第四步: 基于BLAST结果鉴定等位序列

classify.pl -i Eblast.out -p 4 -r ref_gene.gff -g qry_gene.gff
# -p 是等位基因的数目

最终输出的Allele.ctg.table就是后续ALLHiC的prune输入。

方法二: 基于GMAP

第二个方法比较简单，只需要提供多倍体基因组的基因组序列和近缘物种的cds序列即可

第一步: 运行GMAP获取GFF3文件

qry=target.genome  # 你的多倍体序列
ref_cds=reference.cds.fasta # 同源物种的cds
N=4 #你的基因组倍性, 4表示的是4倍体
gmap_build -D . -d DB $qry
gmap -D . -d DB -t 12 -f 2 -n $N  $ref_cds > gmap.gff3

第二步: 获取 allelic.ctg.table

perl gmap2AlleleTable.pl referenece.gff3

gmap2AlleleTable.pl脚本是ALLHiC的一部分。

参考资料

----