首先是用gffread提取cds序列,蛋白序列,转录本序列
gffread genome.gff3 -g genome.fa -x cds.fa
gffread genome.gff3 -g genome.fa -y protein.fa
gffread genome.gff3 -g genome.fa -w transcripts.fa
接下来我们利用组合工具来提取mRNA,和gene序列
python -m jcvi.formats.gff bed --type=mRNA --key=ID a.gff3 > mRNA.bed
这个需要借助Python 包jcvi包
抽取 GFF文件中mRNA, key是ID
基因也是如此
python -m jcvi.formats.gff bed --type=gene --key=ID a.chr.gff3 > gene.bed
image.png
注意用bedtools提取序列，要加-s 参数（能区分正反链）
bedtools getfasta -fi ./a.fasta -bed gene.bed -nameOnly -s -fo a.gene.fa
a.gene.fa就是我们需要的文件，同理也可以得到a.mRNA.fa。
再利用提取最长转录本脚本，获得基因中最长可变剪切的序列。
网友评论