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, ref_gene.gff, qry.fa, qry_gene.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的一部分。
参考资料
- https://github.com/tangerzhang/ALLHiC/issues/16
- https://github.com/tangerzhang/ALLHiC/wiki/ALLHiC:-identify-allelic-contigs
版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。
