构建allelic.ctg.table

ALLHiC简介

ALLHiC算法包括pruning，partition，rescue，optimization，building5个步骤，通过修剪同源染色体之间的交联信号，将等位基因和同源序列分隔在各自的单倍型内独立组装，从而减少了大量拼接错误，通过优化算法改进了contig的排序和定向，尤其是连续性较低的contig，成功解决了染色体水平同源多倍体组装困难的问题。

在多倍体拼接中，比较重要的应该是pruning，该步骤根据提供的Allele.ctg.table过滤HiC—BAM文件中等位基因间坍缩区域和未坍缩区域的HiC信号，并保留最佳链接信号用于后续步骤。

构建allelic.ctg.table

1. based on BLAST

所需数据：拼接contigs(target.genome)，近缘物种CDS序列(reference.cds.fasta)，近缘物种注释信息(reference.gff3)，RNAseq reads(RNAseq_R1.fq.gz，RNAseq_R2.fq.gz)

# STAR alignment

## build index
STAR\
  --runThreadN 10 \
  --runMode genomeGenerate \
  --genomeDir STAR \
  --genomeFastaFiles draft.asm.fa
  
## alignment
mkdir RNA_based && cd RNA_based
STAR\
    --genomeDir path/to/STAR_out \
    --runThreadN 10 \
    --readFilesIn  RNAseq_R1.fq.gz  RNAseq_R2.fq.gz \
    --readFilesCommand zcat \
    --outFileNamePrefix sample\
    --outSAMtype BAM SortedByCoordinate \
    --outBAMsortingThreadN 10 \
    --outSAMstrandField intronMotif \
    --outFilterIntronMotifs RemoveNoncanonical
    
## Spliced transcript
stringtie \
sampleAligned.sortedByCoord.out.bam -p 10 -o qry.gtf

## GTF to GFF3
cufflinks/bin/gffread \
qry.gtf -o qry.gff

## get cdna
cufflinks/bin/gffread \
qry.gff -g draft.asm.fa -w qry.fa

# uniform cds name and gene name
## classify.pl 脚本中是以‘gene’和‘Name=’为锚点识别目标行，这里可以修改代码中的关键词或者修改输入文件
sed -e 's/transcript/gene/' -e 's/ID/Name/' qry.gff > qry_gene.gff

## 同样的，近缘物种也需要修改
## Notice：Ensembl下载的reference.gff3中的gene_name与qry_gene.gff中的gene_Name不是同一个元素。

## 一般情况
ref=reference.gff3
grep '[[:blank:]]mRNA[[:blank:]]' $ref | sed -e 's/mRNA/gene/' -e 's/ID/Name/' > ref_gene.gff

## 例外：
## reference.gff3 最后一列中已经有了Parent_name(Name=)
## 一般不影响classify.pl 脚本识别，但还是建议修改一下
ref=reference.gff3
grep '[[:blank:]]mRNA[[:blank:]]' $ref | sed -e 's/mRNA/gene/' -e 's/Name/geneID/' -e's/ID=transcript:/Name=/' > ref_gene.gff

# blast
## reference.cds.fasta, ref_gene.gff, qry.fa, qry_gene.gff

## build index
blast/bin/makeblastdb \
-in reference.cds.fasta -dbtype nucl

## alignment
blast/bin/blastn \
-query qry.fa -db reference.cds.fasta \
-out qry_vs_ref.blast.out -evalue 0.001 -outfmt 6 -num_threads 4 -num_alignments 1

## filter 
path_to/ALLHiC/scripts/blastn_parse.pl \
-i qry_vs_ref.blast.out -o Eblast.out -q qry.fa -b 1 -c 0.6 -d 0.8 

## Classify alleles 
path_to/ALLHiC/scripts/classify.pl \
-i Eblast.out -p 4 -r ref_gene.gff -g qry_gene.gff

2. Gmap-based

所需数据：拼接contigs(draft.fasta)，近缘物种CDS序列(reference.cds.fasta)，近缘物种注释信息(reference.gff3)

其中近缘物种信息一般可以在ensembl中下载，在这里同样要注意gmap.gff3和reference.gff3中的geneID和Name关键词锚点信息的统一性，否则会获得空白table。

# run gmap get a gff3 file
gmap_build -D . -d DB draft.fasta

gmap -D . -d DB -t 12 -f 2 -n 2 reference.cds.fasta > gmap.gff3

# generate the allelic.ctg.table
path_to/ALLHiC/scripts/gmap2AlleleTable.pl reference.gff3

3. Bed file based

Ensembl下载的近缘物种reference.gff3可能包含比较复杂的成分，我们可以通过通过提取信息的脚本构建一个包含ChrID，Start_Position， End_Position和GeneID的bed文件来增加可检索性。值得注意的是，bed文件中的起始位置是从0开始，并且等于gff文件中起始位置-1。

同样的，gma.gff3中gene和Name所锚定的ID与reference.gff3中gene和Name所锚定的ID并不对应，需要按照上面的方式进行修改。

# generate the allelic.ctg.table
perl gmap2AlleleTableBED.pl ref.bed

gmap2AlleleTableBED.pl.txt

参考

1. Based on BLAST: https://github.com/tangerzhang/ALLHiC/wiki/ALLHiC:-identify-allelic-contigs

2. Gmap-based: https://github.com/tangerzhang/ALLHiC/issues/16

3. ALLHiC续: 如何构建Allele.ctg.table: https://blog.csdn.net/u012110870/article/details/102829550