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

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的一部分。

参考资料

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×