ALLHiC續(xù): 如何構(gòu)建Allele.ctg.table
Allele.ctg.table是ALLHiC用于過濾多倍體基因組中因等位序列相似引起的HiC噪音的必要輸入。
構(gòu)建的方法有很多種，這里列舉兩個方法。
方法一: 基于BLAST
基于BLAST的方法需要你提供一個染色體級別組裝的近緣物種。需要提供4個文件，近緣物種的CDS和GFF文件，多倍體物種的CDS和GFF文件。最重要的是，你得保證CDS序列中的編號和GFF文件的基因編號一致，這也是后續(xù)處理至關重要的一步。
問題一: 我們?nèi)绾慰焖佾@取待組裝物種的CDS和GFF文件？最快速的方法就是利用轉(zhuǎn)錄組序列基于STAR + StringTie的組合進行拼接，并利用gffread從中提取出cdna序列。
假設你的多倍體文件名為contig.fa, 轉(zhuǎn)錄組數(shù)據(jù)為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
# 拼接轉(zhuǎn)錄本
stringtie ${prefix}Aligned.sortedByCoord.out.bam -p $THREADS -o qry.gtf
# GTF 轉(zhuǎn)成 GFF3
gffread qry.gtf -o qry.gff
# 提取cdna
gffread qry.gff -g $ref -w qry.fa
問題二: 后續(xù)的處理腳本要求輸入序列的編號和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文件，考慮到后續(xù)的腳本只認gene所在的行，而cds里面存放的是轉(zhuǎn)錄本序列，我們可以提取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
第三步: 過濾低質(zhì)量的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結(jié)果鑒定等位序列
classify.pl -i Eblast.out -p 4 -r ref_gene.gff -g qry_gene.gff
# -p 是等位基因的數(shù)目
最終輸出的Allele.ctg.table就是后續(xù)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
版權(quán)聲明：本博客所有文章除特別聲明外，均采用 知識共享署名-非商業(yè)性使用-禁止演繹 4.0 國際許可協(xié)議 (CC BY-NC-ND 4.0) 進行許可。
