BED、VCFをスッキリと染色体番号順にソートする方法

こんにちは、あんドーナツ好きで有名な久保（kubor）です。
ブログのネタを日頃から探しているのですが、なかには、過去に先輩が取り上げていそうで、意外にも紹介していない話題が結構あります。

例えば、sortコマンド。
このコマンドは、言わずと知れた基本的なUnixコマンドですが、
社内でもあまり知られていない、素敵なオプションがあります。

このオプションがあればBEDファイルやVCFファイルを綺麗にソートすることが出来ますよ。

では、実際にBEDファイルをソートしてみます。

sort -k 1,1 target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr2 37869024 37899342 CDC42EP3
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 161581163 161583014 AGPAT4-IT1
chr6 39282473 39290330 KCNK16
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2

-kオプションで一列目を指定していますが並び順が変ですね。
ポジションもソートされていませんね。
まずはポジションをソートしましょう。

sort -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2

-k 2,2n オプションで2列目を数値扱いでソートしました。
ポジションはバッチリですね。
あとは、染色体です。「-V」オプションを追加します。

sort -V -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2

綺麗にソートされましたね。
文句なしです。

本来「-V」オプションは、ソフトウェアのバージョン番号をソートするためのオプションですが、染色体番号にも使用可能です。

※ 残念ながらOS Xの場合は、付属のBSD版のsortコマンドに、「-V」オプションがありませんので、GNU版のsortが入っているCoreutils（http://www.gnu.org/software/coreutils/coreutils.html）をインストールする必要があります。

ちなみに、VCFファイルも1列目が染色体番号、2列目がポジションと、BEDファイルの例と同じですので、まったく同じオプションでソート可能です。
ですが、オプションを毎回書いてると面倒なので、エイリアスしておくといいですよ。

sortBedはbedtoolsですでにあるコマンド（http://bedtools.readthedocs.org/en/latest/content/tools/sort.html）ですし、sortFineで、登録しちゃいましょう。

bashをお使いの方は、~/.bashrc
zshをお使いの方は、~/.zshrc
などへ以下を記述します。

alias sortFine="sort -V -k 1,1 -k 2,2n"

必要であれば設定ファイルを再読み込みします。
source ~/.bashrc
or
source ~/.zshrc

それでは、使ってみましょう。

sortFine target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2

うーん、これは相当素敵です。相当ファイン、sortFineですね。