以下の記事の中で、sortコマンドについて触れていました。
多くのコマンドはsorted BEDを要求しますので、事前に以下のコマンドで全てのBEDをソートしておくとストレスがないかと思います。
$ sort -k1,1 -k2,2n in.bed > in.sorted.bed
今日はこのコマンドが何をしているのかについて解説します。
sortについて
sortコマンドの使い方については以下のサイトがとてもわかりやすいです。
sortコマンドについて詳しくまとめました 【Linuxコマンド集】
ここではMACS2のピークファイルを使ってソートを実際に試してみましょう。
narrowPeakをsort
まずはサンプルとなるピークリストのダウンロードと解凍
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz $ gunzip wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz # 内容の表示 $ head wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak chr10 80917404 80917527 . 957 . 149.808548515806 -1 1.78532983501005 57 chr22 39570815 39571059 . 910 . 142.528642559654 -1 1.78532983501005 122 chr18 9643685 9643754 . 849 . 133.01812851997 -1 1.78532983501005 55 chr6 20811494 20811738 . 847 . 132.721884288658 -1 1.78532983501005 122 chr8 103918617 103918861 . 784 . 122.843985314485 -1 1.78532983501005 122 chr3 13692098 13692342 . 768 . 120.292626372352 -1 1.78532983501005 122 chr1 205263766 205264010 . 746 . 116.906226943541 -1 1.78532983501005 122 chr3 58035007 58035251 . 701 . 109.733323372861 -1 1.78532983501005 122 chr17 73781085 73781329 . 694 . 108.721733687985 -1 1.78532983501005 122 chr19 3968693 3968937 . 680 . 106.565969826384 -1 1.78532983501005 122
narrowPeakフォーマットは、染色体、start, end, ピーク名（ID）, score, strand, signalValue, p値、 q値、peak値(0-based offset from chromStart) からなります。
参考： Peak calling with MACS2 | Introduction to ChIP-Seq using high-performance computing
このサンプルは最初からピーク強度が強い順にソートされているようですね。
sort (デフォルト)
何もつけずにsortした結果が以下です。
$ sort wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr1 103716001 103716245 . 113 . 17.757293306914 -1 1.20462049496422 122 chr1 10739728 10739972 . 126 . 19.8727530712812 -1 1.23198670220658 122 chr1 108479111 108479355 . 139 . 21.9096721061029 -1 1.26200103555381 122 chr1 10856653 10856897 . 113 . 17.787081338271 -1 1.20987632444834 122 chr1 109359046 109359290 . 178 . 27.9341585774513 -1 1.37497420775533 122 chr1 110527136 110527380 . 121 . 19.033753770053 -1 1.23745570657267 122 chr1 111470122 111470366 . 120 . 18.8432010665144 -1 1.23273408572432 122 chr1 112946999 112947243 . 378 . 59.2558076715072 -1 1.53329986093226 122 chr1 113243465 113243709 . 114 . 17.9824166190626 -1 1.21252694920448 122 chr1 119623228 119623472 . 338 . 52.9747880171335 -1 1.74485156130976 122
左の列から優先で昇順でソートされるので、chr1が最初に来るようになります。ちなみにchr1の次はchr11です。
特定の列を基準にsort
仮に2列目の"start" を基準にソートしたい場合は、以下のように -k2,2 をつけます。
(-k2 だけだと2列目以降全ての列が評価されてしまいます。-k2,2 は「2列目から2列目まで」という指定です。)
$ sort -k2,2 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr14 100003838 100004082 . 162 . 25.3896151955265 -1 1.32856874236087 122 chr8 10009854 10010098 . 140 . 22.0324635890778 -1 1.26598847143077 122 chr15 100106095 100106339 . 165 . 25.8864092879034 -1 1.33959339080326 122 chr13 100150800 100151044 . 144 . 22.5409997536205 -1 1.27478892238628 122 chr15 100243871 100244115 . 287 . 45.0594193501982 -1 1.56124172764321 122 chr13 100340050 100340294 . 431 . 67.4827532683693 -1 1.68424674749618 122 chr13 100370580 100370824 . 134 . 20.9796948179283 -1 1.25368035557817 122 chr9 100394971 100395215 . 297 . 46.6277802709549 -1 1.52770570059048 122 chr7 100434051 100434295 . 208 . 32.5983076837242 -1 1.4658668221046 122 chr11 100449778 100450022 . 144 . 22.5396554016715 -1 1.27512069817163 122
なにか変ですね。"100003838"が最初に来ています。これは数値が文字列として評価されているためです。
数値としてソートしたい場合は -n をつけましょう。
$ sort -k2,2 -n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr16 122590 122834 . 156 . 24.5519243359907 -1 1.31865948396838 122 chr16 441938 442182 . 113 . 17.8207002619443 -1 1.20806676148905 122 chr9 444599 444843 . 191 . 30.0193503293631 -1 1.4152662953559 122 chr4 668097 668341 . 138 . 21.6323169668381 -1 1.25420898826128 122 chr16 790836 791080 . 122 . 19.1922172413508 -1 1.24286158989763 122 chr17 855058 855302 . 190 . 29.8288083711511 -1 1.41876633698804 122 chr1 936188 936432 . 278 . 43.6111912894822 -1 1.58313206719137 122 chr18 939233 939477 . 259 . 40.5622471525867 -1 1.48685535527511 122 chr7 985313 985557 . 133 . 20.8488735208324 -1 1.25369923175399 122 chr19 996416 996660 . 330 . 51.7245477629659 -1 1.65959931242927 122
無事、数値としてソートされました。
ちなみに -k2,2 -n は -k2,2n と書いても同じ意味になります。
$ sort -k2,2n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr16 122590 122834 . 156 . 24.5519243359907 -1 1.31865948396838 122 chr16 441938 442182 . 113 . 17.8207002619443 -1 1.20806676148905 122 chr9 444599 444843 . 191 . 30.0193503293631 -1 1.4152662953559 122 chr4 668097 668341 . 138 . 21.6323169668381 -1 1.25420898826128 122 chr16 790836 791080 . 122 . 19.1922172413508 -1 1.24286158989763 122 chr17 855058 855302 . 190 . 29.8288083711511 -1 1.41876633698804 122 chr1 936188 936432 . 278 . 43.6111912894822 -1 1.58313206719137 122 chr18 939233 939477 . 259 . 40.5622471525867 -1 1.48685535527511 122 chr7 985313 985557 . 133 . 20.8488735208324 -1 1.25369923175399 122 chr19 996416 996660 . 330 . 51.7245477629659 -1 1.65959931242927 122
降順でソートしたい場合は -r オプションをつけましょう。
$ sort -k2,2 -n -r wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr1 244436392 244436636 . 116 . 18.2779430757098 -1 1.22768405980467 122 chr1 244327430 244327674 . 236 . 37.0352452557144 -1 1.49136169383242 122 chr1 244090184 244090428 . 198 . 31.1366875436721 -1 1.45062610008222 122 chr1 242097512 242097756 . 219 . 34.351615073636 -1 1.46008791542184 122 chr1 239979463 239979707 . 144 . 22.5476334083494 -1 1.27445689294964 122 chr1 238811236 238811480 . 187 . 29.3108873795414 -1 1.42587671098019 122 chr2 238522973 238523217 . 206 . 32.3223759856553 -1 1.47135994807226 122 chr1 237360863 237361107 . 129 . 20.2738825008399 -1 1.24189631661435 122 chr2 237118258 237118502 . 183 . 28.6979797945465 -1 1.3746686662588 122 chr2 236688138 236688382 . 154 . 24.1282500315595 -1 1.31204005182904 122
-k2,2 -n -r は -k2,2nr と書くこともできます。 -k2,2n -r とは書けないようです。
染色体配置にもとづくソート
染色体上での並び順にもとづいてピークをソートしたい場合、まず1列目と2列目でソート（1列目優先）となります。これは以下のように書けます。
$ sort -k1,1 -k2,2n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr1 936188 936432 . 278 . 43.6111912894822 -1 1.58313206719137 122 chr1 1790095 1790339 . 147 . 23.0760153956621 -1 1.27356690138484 122 chr1 2143910 2144154 . 328 . 51.3985090086356 -1 1.67414642060339 122 chr1 3530010 3530254 . 195 . 30.5372043198762 -1 1.41709407864103 122 chr1 6333445 6333689 . 116 . 18.1876867422188 -1 1.22125052352078 122 chr1 6844488 6844732 . 167 . 26.1507228133338 -1 1.3422008190137 122 chr1 7669033 7669277 . 132 . 20.7725171056509 -1 1.24639610451395 122 chr1 7764587 7764831 . 125 . 19.6798804866856 -1 1.24285094687436 122 chr1 8933714 8933958 . 209 . 32.8047056184361 -1 1.45819850201426 122 chr1 9488912 9489156 . 300 . 47.0204903548675 -1 1.50753840495668 122
-k オプションを複数指定すると、左から優先でファイルがソートされます。１列目は文字列で２列目が数値なので、-k2,2の方だけに-n オプションがついています。
これが冒頭で示した、BEDtoolsで要求されるソートになります。
ピーク強度でソート
ピークが強い順にtop 1000ピークを抽出したい、というような場合は、以下のようになります。
MACS2の場合はq値がピーク強度を表しているので、9列目でソートします。
$ sort -k9,9nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak |head chr1 156643538 156643782 . 584 . 91.4369629380883 -1 1.78532983501005 122 chr1 205263766 205264010 . 746 . 116.906226943541 -1 1.78532983501005 122 chr1 214489937 214490181 . 556 . 87.1055294122947 -1 1.78532983501005 122 chr1 30447209 30447453 . 623 . 97.5499934901461 -1 1.78532983501005 122 chr10 65389313 65389557 . 649 . 101.626596581112 -1 1.78532983501005 122 chr10 80917404 80917527 . 957 . 149.808548515806 -1 1.78532983501005 57 chr10 92079336 92079580 . 569 . 89.1329212789991 -1 1.78532983501005 122 chr11 65264905 65265149 . 547 . 85.6244766911001 -1 1.78532983501005 122 chr11 86451724 86451968 . 565 . 88.4367426823928 -1 1.78532983501005 122 chr11 9634736 9634980 . 586 . 91.8571732244745 -1 1.78532983501005 122
同じq値をもつピークが（何故か）多くあり、同値でうまくソートできませんでした。
では、signalValueを表す7列目を追加で指定してみましょう。
$ sort -k9,9nr -k7,7nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head chr10 80917404 80917527 . 957 . 149.808548515806 -1 1.78532983501005 57 chr22 39570815 39571059 . 910 . 142.528642559654 -1 1.78532983501005 122 chr18 9643685 9643754 . 849 . 133.01812851997 -1 1.78532983501005 55 chr6 20811494 20811738 . 847 . 132.721884288658 -1 1.78532983501005 122 chr8 103918617 103918861 . 784 . 122.843985314485 -1 1.78532983501005 122 chr3 13692098 13692342 . 768 . 120.292626372352 -1 1.78532983501005 122 chr1 205263766 205264010 . 746 . 116.906226943541 -1 1.78532983501005 122 chr3 58035007 58035251 . 701 . 109.733323372861 -1 1.78532983501005 122 chr17 73781085 73781329 . 694 . 108.721733687985 -1 1.78532983501005 122 chr19 3968693 3968937 . 680 . 106.565969826384 -1 1.78532983501005 122
無事、思い通りのソートになりました。
sort -k9,9nr -k7,7nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head -n1000 とすれば、top1000ピークが手に入ります（ヘッダ行が含まれない場合。）
注意点
ソートする文字列の中に記号や日本語など複雑な文字列が含まれている場合のソート順は環境に依存します。複雑なファイルのソートの場合、場合によっては使ったOSなどでsortの結果が違うかもしれないので、ファイル比較などをする際は統一環境でソートしておくほうが無難です。
