Biostarとbedtoolsの公式サイトに、bedtoolsを使ったpromoter配列の抽出の仕方がまとめられている（How To Use Bedtools To Extract Promoters From A Mouse Bed File）。試してみる。
このようなbedファイルの各featureの上流を抽出する。
user$ head -5 genes.bed
chr1 11873 12227 NR_046018_exon_0_0_chr1_11874_f
chr1 12612 12721 NR_046018_exon_1_0_chr1_12613_f
chr1 13220 14409 NR_046018_exon_2_0_chr1_13221_f
chr1 14361 14829 NR_024540_exon_0_0_chr1_14362_r
chr1 14969 15038 NR_024540_exon_1_0_chr1_14970_r
まずクロモソームのサイズを書いたファイルを準備する。単純なタブ区切りのテキストファイルであればよい。ここではsamtools faidxを使いfastaから作っている。
samtools faidx reference.fa
cut -f 1,2 reference.fa.fai > genome.txt
上流領域を抽出する。bedtoolsで隣接 featureを取れるコマンドは複数あるが、今回はflankコマンドを使う。
bedtools flank -i genes.bed -g genome.txt -l 1000 -r 0 -s > genes.1kb.promoters.bed
flankコマンドで隣接領域を抽出している。この時-l 1000 -r 0とすることで隣接する上流1000-bp、下流0-bpが抽出されるように工夫している。ちなみにクロモソームの末端でクロモソームからはみ出すfeatureがある場合、トリミングして出力される。
さらにbedtoolsのgetfastaで配列を抽出する。
bedtools getfasta -fi reference.fa -bed genes.1kb.promoters.bed |fold -w 60 > genes.1kb.promoters.bed.fa
- -fi Input FASTA file
- -bed BED/GFF/VCF file of ranges to extract from -fi
-fiでリファレンスゲノムを指定して使う。また、パイプを使いfoldで60文字ごとに改行している。少し重くなるので、不要なら外す。上記コマンドはヘッダー名が違うとワークしないので気をつける。
出力は、例えば grep -n "^>" input.fasta | wc -l で配列数をカウントできる。
使用したpromoter配列の塩基組成や偏りなどの基本情報をまとめたい場合、bedtoolsのnucを使う。
先ほど作ったbedを使い、以下のように打つ。
bedtools nuc -fi reference.fa -bed genes.1kb.promoters.bed > ATGC.bed
カラムの後半に以下のような情報がつく。
- 1) %AT content
- 2) %GC content
- 3) Number of As observed
- 4) Number of Cs observed
- 5) Number of Gs observed
- 6) Number of Ts observed
- 7) Number of Ns observed
- 8) Number of other bases observed
- 9) The length of the explored sequence/interval.
- 10) The seq. extracted from the FASTA file. (opt., if -seq is used)
- 11) The number of times a user's pattern was observed.
平均や最大/最小などの情報はRでまとめる。
上で使ったbedtoolsのコマンドはこちらでまとめています。
引用
How To Use Bedtools To Extract Promoters From A Mouse Bed File
