下準備

テスト用に使うデータをダウンロードしておきましょう．

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のワンライナーを使うとさらに強力な処理が可能
- ただ，私はあまり使わない．必要に迫られたときだけ使う