VCFファイルから周辺配列（flanking sequence）を取得したいとき、様々な方法が考えられると思います。今回はその中で bedtools を用いた方法を考えていきたいと思います。
bedtools flank というコマンドが周辺配列をとってくる際に用いるコマンドですが、このコマンドはリファレンスの配列ごとの塩基長が記されたファイルを必要とします。このファイルは以下のサイトを参考に以下のコマンドで作成することができます。
https://www.danielecook.com/generate-fasta-sequence-lengths/
cat file.fa | awk '$0 ~ ">" {if (NR > 1) {print c;} c=0;printf substr($0,2,100) "\t"; } $0 !~ ">" {c+=length($0);} END { print c; }' > my.genome
また、今回はSNP情報が記載されたVCFファイルをスタートに進めていこうと思います。そこで最初にVCFをBED形式に変換してあげます。今回は bcftools query を用いました。コマンドは以下のようになります。
bcftools query -f'%CHROM\t%POS0\t%END\t%ID%REF\n' variant.vcf >variant.bed
そしてBEDファイルができたら早速 bedtools flank で周辺配列の座標情報を取得しBEDファイルに出力します。
bedtools flank -i variant.bed -g my.genome -b 100 > flank.bed
奇数行と偶数行に分けることでSNPの左側の情報と右側の情報とに分けます。
awk 'NR%2!=0' flank.bed > left.bed awk 'NR%2==0' flank.bed > right.bed
最後に bedtools getfasta で配列情報を出力させて完成です。
bedtools getfasta -fi genome.fa -bed left.bed bedtools getfasta -fi genome.fa -bed right.bed
