需求:提取 1 号染色体的 21856782-21856982 和 43942666-43942866 区间的 fa 序列;
需要准备的test.gff3
文件如下所示,第一列是染色体号,第二列和第三列的内容随意,第四列和第五列是提取的开始位置和终止位置,第六列和第八列内容随意,第七列是基因组的正负链,第九列序列的信息,第九列必须有Parent=
, ;ID=
, ;gene_id=
信息,不然序列提取不出来。分隔符为tab:
chr1 gene transcript 21856782 21856982 . + . Parent=ENST0000001234.1;ID=ENST0000001234.1;gene_id=ENST0000001234.1
chr1 gene transcript 43942666 43942866 . + . Parent=ENST00000567890.1;ID=ENST00000567890.1;gene_id=ENST00000567890.1
另外还需要基因组文件GRCh38.primary_assembly.genome.fa
, 如下所示:
准备好以上两个文件后,用命令gffread -w test.fa -g GRCh38.primary_assembly.genome.fa test.gff3
即可获取需要的序列。序列存储在test.fa
文件中。