Bashでシーケンシングデータを扱う

このトピックを終えたら

ターミナルでFASTA/FASTQファイルを開いて、配列数を数えて、特定の配列を検索できるようになります。マウスなしでキーボードだけでデータを扱う初めての体験です。

ターミナルとは？

研究室でピペットを使って試薬を精密に扱うように、ターミナルはコンピュータのデータを精密に扱うツールです。クリックの代わりにコマンドを入力します。

Google Colabでは、セルの先頭に!を付けるとBashコマンドを実行できます。

ファイルを作成して確認する

まず、練習用のFASTAファイルを作りましょう。

bash

# 練習用FASTAファイルを作成
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

# ファイルの内容を確認
cat genes.fasta

bash

# 行数を数える
wc -l genes.fasta

wcはWord Countの略です。-lオプションは行数のみを数えます。

grep: 配列を検索する

grepはファイルから特定のパターンを見つけるコマンドです。論文でキーワードを検索するのと同じです。

bash

# FASTAヘッダー(>)の行のみを出力
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta

bash

# 配列数を数える（ヘッダー行 = 配列数）
grep -c ">" genes.fasta

-cオプションはマッチした行の数だけを出力します。

bash

# 特定の遺伝子を検索
grep "BRCA1" genes.fasta

FASTQ形式を理解する

実際のシーケンシング装置から出てくるデータはFASTQ形式です。FASTAと違い、品質スコア（quality score）が含まれています。

bash

# FASTQファイルを作成（4行で1セット）
echo "@read_001
ATGCGATCGATCGATCGATCG
+
IIIIIIIIIIIIIIIIIIIII
@read_002
GCTAGCTAGCTAGCTAGCTAG
+
IIIIIIIIIIIIIIIIIIIII
@read_003
NNNNNATCGATCGATCGATCG
+
!!!!!IIIIIIIIIIIIIIIII" > reads.fastq

cat reads.fastq

FASTQの4行構造：

@で始まるリードID
塩基配列
+（区切り線）
品質スコア（高いほど良い、I=最高、!=最低）

awk: データを加工する

awkはテキストデータを列（column）単位で処理するツールです。スプレッドシートで特定の列だけを選択するのに似ています。

bash

# タブ区切りデータを作成
echo "BRCA1	chr17	43044295	43170245
TP53	chr17	7661779	7687538
EGFR	chr7	55019017	55211628" > gene_locations.tsv

# 遺伝子名（1列目）のみを出力
awk '{print $1}' gene_locations.tsv

bash

# 遺伝子の長さを計算（4列目 - 3列目）
awk '{print $1, $4 - $3, "bp"}' gene_locations.tsv

パイプ(|): コマンドをつなげる

パイプ（|）はあるコマンドの出力を次のコマンドの入力に接続します。実験で試薬Aの生成物を試薬Bに入れるのと同じです。

bash

# FASTAからヘッダーだけを抽出して数を数える
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta | wc -l

やってみよう（Faded Example）

空欄を埋めて、FASTQファイルのリード数を数えるコマンドを完成させてください。

FASTQではリードIDは@で始まります。

穴埋め問題bash

# @で始まる行の数を数えるには？
grep  "" test.fastq

よくあるエラーと解決法

Q: command not foundが出ます

Colabでは!プレフィックスが必要です：!grep ">" genes.fasta

Q: 文字化けします

ファイルのエンコーディングの問題です。file genes.fastaでエンコーディングを確認し、UTF-8でなければiconv -f EUC-KR -t UTF-8 input.txt > output.txtで変換します。

Q: Permission deniedが出ます

ファイルに実行権限がない場合に発生します。chmod +x script.shで権限を付与してください。

次の記事では、Pythonでこのデータをより精密に分析する方法を学びます。

BioPlayground

Bashでシーケンシングデータを扱う

Bashでシーケンシングデータを扱う

このトピックを終えたら

ターミナルとは？

ファイルを作成して確認する

grep: 配列を検索する

FASTQ形式を理解する

awk: データを加工する

パイプ(|): コマンドをつなげる

やってみよう（Faded Example）

よくあるエラーと解決法

💬 質問・コメント