有個師弟說他想要將對某參考基因組的各條染色體生成相同的bin區(qū)間，以bed格式儲存結(jié)果。
即參考序列的GFF格式文件為：
gff.jpg
希望根據(jù)參考序列進行操作，輸出的結(jié)果為以fasta格式儲存結(jié)果：
Chr bin_start bin_end
1 1 2000
1 2001 4000
1 4001 6000
...
2 1 2000
2 2001 4000
（Bed格式是一種0 based的文件，這種寫法是不對的...）
其實這個文件，寫腳本非常簡單，但是也有現(xiàn)成的輪子，即熟練使用samtools/bedtools解決這種需求，輕而易舉。
比如我們現(xiàn)在有個fake參考序列，fake_genome.fa，內(nèi)容如下：
>Chr1 57bp
GTTTGGTTTGTGCGTGATGTTAAGATCGGAAGAGCACACGTCTGGAGCACACGTCTG
>Chr2 50bp
NCCTCTGCAAACGGGTCTGATAGTATTTCAGATCGGAAGAGCACACGTCT
現(xiàn)在我們想將其切割成長度為20bp的一條一條的bins區(qū)間，并且用Bed格式保存結(jié)果用于后續(xù)分析。
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
最后輸出結(jié)果如下：
Chr1 0 20
Chr1 20 40
Chr1 40 57
Chr2 0 20
Chr2 20 40
Chr2 40 50
實現(xiàn)了師弟最初的想法。
也可以將某條染色體上的bins區(qū)間分別輸出成bed， 代碼如下：
bedtools makewindows -g fake_genome.chrome.size -w 20 | awk '{if ($1=="Chr1"){print $0 >> $1":"$2"-"$3".bed"}}'
