下準備
テスト用に使うデータをダウンロードしておきましょう．
curl -O http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refFlat.txt.gz
gunzip refFlat.txt.gz
このデータは人のゲノムの中でどこにどの遺伝子があるかをリストアップしたもの．
データの各カラムの説明
- geneName (遺伝子シンボル)
- name (RefSeq ID)
- chrom (染色体)
- strand (ストランド)
- txStart (転写開始点)
- txEnd (転写終了点)
- cdsStart (Coding DNA Sequence Start)
- cdsEnd (Coding DNA Sequence Start)
- exonCount (exonの数)
- exonStarts (exonの開始地点．カンマ区切り)
- exonEnds (exonの終了地点．カンマ区切り)
テキスト処理に使うコマンド一覧
-
grep: 特定の文字列を含む行を抜き出す -
cat: ファイルの連結 -
head: ファイルの先頭を抜き出す -
tail: ファイルの末尾を抜き出す -
sed: 置換やその他の様々な操作が行える -
tr: 1文字置換 -
cut: 特定の列を抜き出す -
paste: 列を結合する -
wc: Word Count -
sort: 並び替え -
uniq: 重複を排除する
今までに解説済みもあるが一部は改めて解説する
grep
特定の文字を含む行を抜き出す．正規表現が使える
$ grep chr1 refFlat.txt
CCHE1 NR_131782 chr10 - 59173962 59176464 59176464 59176464 1 59173962, 59176464,
JMJD1C NM_001282948 chr10 - 63167220 63269223 63168044 63219884 25 63167220,63168434,63176296,63177716,63183446,63184607,63185562,63186214,63189167,63190893,63192937,63193344,63194285,63197410,63198512,63200475,63206594,63209062,63213472,63215262,63215552,63217206,63219877,63264650,63268662, 63168134,63168566,63176473,63177856,63183569,63184738,63185653,63186383,63189446,63191108,63193151,63193472,63194375,63197563,63198727,63200677,63208801,63209235,63215151,63215455,63215696,63217331,63219983,63264764,63269223,
適当にchr1を含む行だけを抜き出そうとするとchr10を含む行も一緒に出てきてしまう
$ grep '[[:space:]]chr1[[:space:]]' refFlat.txt
HSD3B2 NM_000198 chr1 + 119415119 119423039 119415419 119422620 4 119415119,119415330,119419417,119421808, 119415202,119415561,119419582,119423039,
ANKRD65 NM_001243536 chr1 - 1418419 1421444 1419541 1421005 3 1418419,1420815,1421132, 1419549,1421005,1421444,
[[:space:]]はいかなるスペースにもマッチする．この例ではchr1の前後が空白である行にだけマッチしている
正規表現 1
- ゆらぎを許したパターンを判定させる
- 有限オートマトンで受理可能な範囲でしか書けない
- なので，括弧の対応関係は書けない（証明できますか？）
- （ただし，拡張の種類次第で書けたりする）
- 標準的な書き方がある
- 詳しくは
man re_formatもしくはman 7 regex
主な正規表現 1
-
.: 任意の1文字 -
*: 直前の文字の0回以上の繰り返し -
+: 直前の文字の1回以上の繰り返し -
?: 直前の文字の0回もしくは1回の繰り返し -
[HOGE]: H, O, G, Eのどれでも良い- 括弧の中の任意の1文字にマッチする
-
[a-z]とするとaからzまでの任意の1文字にマッチする - 文字を並べる代わりに
[:alpha:]と書けばアルファベットのみ -
[^HOGE]とするとH, O, G, E以外になる
主な正規表現 2
-
^: 行の先頭 -
$: 行の末尾 -
(aa|bb):aaもしくはbbにマッチ
拡張モードを利用しない場合は+, ?, (|)は\を前においてエスケープする必要がある．grepで拡張モードを利用するためには-Eをオプションとして与えるか，egrepを代わりに使う．正規表現を使いたくない場合には-Fを付けるか，fgrepを代わりに使う
正規表現の詳細はman re_formatもしくはman 7 regexを参照すること
課題
5番染色体と10番染色体だけを抜き出しなさい
解答
$ grep '[[:space:]]chr\(5\|10\)[[:space:]]' refFlat.txt
CCHE1 NR_131782 chr10 - 59173962 59176464 59176464 59176464 1 59173962, 59176464,
JMJD1C NM_001282948 chr10 - 63167220 63269223 63168044 63219884 25 63167220,63168434,63176296,63177716,63183446,63184607,63185562,63186214,63189167,63190893,63192937,63193344,63194285,63197410,63198512,63200475,63206594,63209062,63213472,63215262,63215552,63217206,63219877,63264650,63268662, 63168134,63168566,63176473,63177856,63183569,63184738,63185653,63186383,63189446,63191108,63193151,63193472,63194375,63197563,63198727,63200677,63208801,63209235,63215151,63215455,63215696,63217331,63219983,63264764,63269223,
...
$ grep -E '[[:space:]]chr(5|10)[[:space:]]' refFlat.txt
...
$ grep -E '\schr(5|10)\s' refFlat.txt
...
-wをつけることで正規表現の前後が[[:<:]]で囲まれたことと同じになり，単語ごとにマッチさせることが可能になる．
$ grep -w -E 'chr(5|10)' refFlat.txt
head / tail
- ファイルの先頭と末尾を表示する
-
-nオプションで表示する行数を指定できる- 省略したら 10 行
-
tailは-fを付けることでファイルを監視できる
$ tail -f /var/log/system.log
...
/var/log/system.logに行がが追加されるたびに出力される．終了はControl-C
tr (translate)
1文字置換を行う
以下は0をaに，1をbに置き換える例
$ cat refFlat.txt|tr 01 ab
CCT8L2 NM_ab44a6 chr22 - b659a757 b65928ba b659a876 b659255a b b659a757, b65928ba,
CCHEb NR_b3b782 chrba - 59b73962 59b76464 59b76464 59b76464 b 59b73962, 59b76464,
...
あんまり使わない
cut
- 特定の列を抜き出す
-
-fで抜き出す列番号を指定する - 複数指定するときは
,で区切る
-
- 以下は1列目と3列目だけを抜き出している
$ cut -f 1,3 refFlat.txt
CCT8L2 chr22
CCHE1 chr10
JMJD1C chr10
TIFA chr4
...
paste
- 列を結合する
$ cut -f 1 refFlat.txt > column-1.txt
$ cut -f 3 refFlat.txt > column-3.txt
$ paste column-1.txt column-3.txt|head
CCT8L2 chr22
CCHE1 chr10
JMJD1C chr10
TIFA chr4
...
あんまり使わないかも
sort
- 特定の列でソートする
sort refFlat.txt # 遺伝子シンボルでソート
A1BG NM_130786 chr19 - 58346805 58353499 58347021 58353437 8 58346805,58347352,58350369,58351390,58352282,58352927,58353291,58353403, 58347029,58347640,58350651,58351687,58352555,58353197,58353327,58353499,
A1BG-AS1 NR_015380 chr19 + 58351969 58355183 58355183 58355183 4 58351969,58353378,58353713,58354368, 58353044,58353474,58353857,58355183,
A1CF NM_001198818 chr10 - 50799408 50885675 50806728 50859940 14 50799408,50809893,50811039,50813856,50816005,50820551,50828130,50836073,50841861,50843987,50859841,50862888,50864032,50885580, 50806880,50810042,50811176,50814038,50816279,50820649,50828295,50836312,50841992,50844122,50859985,50862981,50864080,50885675,
...
- 列番号を指定するためには
-kを使う - そのままだと文字列でソートされるので
10が2の前に来てしまう- 数字としてソートするためには
-nをつける -
sort -n -k 5 refFlat.txtとsort -k 5 refFlat.txtの結果を比較してみる
- 数字としてソートするためには
uniq
- 重複のないユニークな行を出力できる
- 入力はソートされていることが前提となる
-
-cを付けると何回出現したかを表示できる
$ cut -f 1 refFlat.txt|sort|uniq|wc -l
26433
26433遺伝子のデータがあることが分かる
sort|uniqはsort -uでも代用可
課題
Exon の数の度数分布表を作りなさい．ただし，出現頻度順にソートされているものとする
- 度数分布表: 各値が出現した回数を並べた表
解答
$ cut -f 9 refFlat.txt|sort|uniq -c|sort -n
1 102
1 103
1 104
1 107
...中略...
4743 5
5270 4
5281 3
5721 1
sed (stream editor)
- テキスト編集が何でもできる
- 難解
- 特定のパターンだけ覚えればOK
- 一番よく使うのは置換
- 以下の例は
chrをChromosomeに置換する- 置換には正規表現が使える
$ sed -e 's/chr/Chromosome/g' refFlat.txt
CCT8L2 NM_014406 Chromosome22 - 16590757 16592810 16590876 16592550 1
16590757, 16592810,
CCHE1 NR_131782 Chromosome10 - 59173962 59176464 59176464 59176464 1
59173962, 59176464,
- 全部は解説できないので
man sedをどうぞ
まとめ
- テキスト処理はコマンドをパイプでつなぎ様々なコマンドを組み合わせることでプログラムを書くことなく色々な処理ができる
- テキスト処理には
grep,cat,head,tail,sed,tr,cut,paste,wc,sort,uniqなどを使う
おまけ
-
perlやawkのワンライナーを使うとさらに強力な処理が可能- ただ，私はあまり使わない．必要に迫られたときだけ使う