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.gzread_R2.fq.gz

mkdir -p STAR
ref=contig.fa
THREADS=40
fq1=read_R1.fq.gz
fq2=read_R2.fq.gz
prefix=sample

问题二: 后续的处理脚本要求输入序列的编号和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

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

方法二: 基于GMAP

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

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

qry=target.genome

第二步: 获取 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

阅读原文

ALLHiC续: 如何构建Allele.ctg.table》来自互联网,仅为收藏学习,如侵权请联系删除。本文URL:http://www.hashtobe.com/522753.html