有个师弟说他想要将对某参考基因组的各条染色体生成相同的bin区间,以bed格式储存结果。
即参考序列的GFF格式文件为:希望根据参考序列进行操作,输出的结果为以fasta格式储存结果:
Chr bin_start bin_end
1 1 2000
1 2001 4000
1 4001 6000
...
2 1 2000
2 2001 4000
(Bed格式是一种0 based的文件,这种写法是不对的...)
其实这个文件,写脚本非常简单,但是也有现成的轮子,即熟练使用samtools/bedtools解决这种需求,轻而易举。
比如我们现在有个fake参考序列,fake_genome.fa,内容如下:
>Chr1 57bp
GTTTGGTTTGTGCGTGATGTTAAGATCGGAAGAGCACACGTCTGGAGCACACGTCTG
>Chr2 50bp
NCCTCTGCAAACGGGTCTGATAGTATTTCAGATCGGAAGAGCACACGTCT
现在我们想将其切割成长度为20bp的一条一条的bins区间,并且用Bed格式保存结果用于后续分析。
samtools faidx fake_genome.fa
cut -f 1,2 fake_genome.fa.fai > fake_genome.chrome.size
bedtools makewindows -g fake_genome.chrome.size -w 20 > result.bed
最后输出结果如下:
Chr1 0 20
Chr1 20 40
Chr1 40 57
Chr2 0 20
Chr2 20 40
Chr2 40 50
实现了师弟最初的想法。
也可以将某条染色体上的bins区间分别输出成bed, 代码如下:
bedtools makewindows -g fake_genome.chrome.size -w 20 | awk '{if ($1=="Chr1"){print $0 >> $1":"$2"-"$3".bed"}}'