複数のPCR産物をNanoporeでシークエンス
目的
変異株のスクリーニングを行い、いくつか候補を得ることができました。 次は原因遺伝子の同定に進みたいのですが、既知の遺伝子のアリルかもしれません。 その場合、候補遺伝子をサンガー法でシークエンスして、新規であることを確認しておきたいところです。 ですが、候補が複数あったり、大きかったりすると、いくつもの領域に分けてPCRをして、それぞれをシークエンスをする必要があります。 まあまあのコストがかかります。
Azenta Life Sciencesという会社は、PCR-EZというサービスを提供しています。 「単一なクローンで直鎖状」PCR産物を、(おそらく)Nanoporeを使ってシークエンスしてくれるものです。 1サンプルあたり4000円(2026年10月現在)と、サンガーシークエンスの10回分にも満たない価格です。
このサービスに、いくつかのPCR産物を混ぜたものを1サンプルとして送付したところ、リードがアセンブルされて、もっとも長い領域についてのみ解析がされました。 そこで、生のリード(fastq)を解析して、複数の遺伝子の変異の有無を調べることにしました。
この方法の一つよいところは、ヘテロの挿入・欠失が同定できるということです。 サンガー法のシークエンスではぐちゃぐちゃになって読めないところが、次世代シーケンサーは1分子ごとに読むので、ヘテロもきれいに調べることができます。
Nanopore
(簡単な説明をいずれ書く)
考え方
多数のリードをリファレンスと比較して、アライメントを作成する。この部分は自分でプログラムを書けないので、既存のツールを使う。
- こういうアライメントがつくれれば、リファレンスの各塩基ごとにリードを数えることができる。ミスマッチを見つけることができる。
read: ACGTAGAAA-CAGCTATTT
ref : ACGTACAAAACAGGTATTT
- 解析システムの前半はシーケンサーの出力からアライメントをつくる部分。後半はアライメントからリードを数える部分に分けることができる。
アライメントまで
azentaのデータはfastqファイル。nanoporeなので、かなり長い。平均で1000bpぐらい。リード数はサンプル当たり600弱。 bowtieでやってみると1000bpを超えるリードはあてられない。 BLASTを使うことにする。
BLASTは入力の形式がfasta。fastqからfastaに変換する必要がある。
前半手順
- BLASTのインデックスを作成する。makeblastdb -in -out -dbtype nucl リファレンスの起き場所と、インデックスの起き場所を決める必要がある。referenceディレクトリを使うことにする。
- fastq.gzを展開する zcatならパイプでつぎの処理へ。中間ファイルをつくるならgunzipで展開する。
- fastaに変換する。1行目の先頭@を>に置換し、3行目と4行目を削除する。sedを使う。そうするとzcatの出力からパイプでsedにつないで、そのままblastnにもっていくことができる。
- blastでアライメントを作成する。リファレンスとリードを比較する必要があるので、-outfmt 5のオプションをつけて、XMLで出力する。
後半手順
- XMLからアライメント部分を抽出する。こういうのをパース(parse)という。自分でプログラムを書くにせよ、ツールを使うにせよ、xpathを知っておく必要がある。この後で、リードを数え上げるのに都合のよい形式に変換する。
- 数え上げる。プログラムを書く。この後図にすることを考えるとRを使うのがよいか。
- Rで図をつくる。
- 図、その他を一覧できるようにHTMLをつくる。
アライメント作成までの手順
BLASTのインデックスをつくる
とりあえずこのシステム開発を行うディレクトリをつくり、そこに移動する。以後はそのディレクトリで行う。
リファレンスをダウンロードする。ウェブからファイルをダウンロードするのにcurlを使う。-OをつけるとURLからファイル名を決めてくれる。
$ curl -O https://kiku3.tsbio.info/~takamasa/260917Honda/ref/SUDL.fa $ mkdir reference $ mv SUDL.fa reference
blastをインストールする。 パッケージは apt update; apt search ncbi-blast; で探す。 apt updateはパッケージのデータベースを更新するコマンドなので、1週間に1回でもやれば十分。
見つかったら apt install ncbi-blast+ でインストールする。 管理者権限が必要なので適宜sudoをつけて実行する。
$ sudo apt update $ sudo apt search ncbi-blast $ sudo apt install ncbi-blast+
リファレンスを作成する。
$ makeblastdb -in reference/SUDL.fa -out reference/SUDL -dbtype nucl $ ls reference SUDL.fa SUDL.ndb SUDL.nhr SUDL.nin SUDL.not SUDL.nsq SUDL.ntf SUDL.nto
SUDL.ninなどができていればよい。
fastqからfastaへ
サンプルデータ(azentaが納品するデータ)をダウンロードし、展開する。
$ curl -O https://kiku3.tsbio.info/~takamasa/260917Honda/60-1380693045_FASTQ.zip $ unzip 60-1380693045_FASTQ.zip
sedでfastqをfastaに加工するコマンドをつくる。
- 1行目から4行ごとに(1~4)、先頭の@を>に置換(s/^@/>/)する。-> 1~4s/^@/>/
- 3行目から4行ごとに(3~4)、行を削除(d)する。-> 3~4d
- 4行目から4行ごとに(4~4)、行を削除(d)する。-> 4~4d
これらのコマンドを;でつなぐ。
$ mkdir tmp $ zcat raw_fastq/2p-7.fastq.gz | sed -e "1~4s/^@/>/;3~4d;4~4d" >tmp/2p-7.fasta
blast検索
塩基配列同士の検索なので、blastnを使う。オプションの意味は下記の通り。
- -out 出力ファイル名
- -outfmt 5 出力形式をXMLにする
- -max_target_seqs 1 最も一致度が高い配列だけを出力する。
- -max_hsps 1 アライメントを一つだけ出力する。リファレンスが環状(プラスミド)の場合、max_hspsは2にする必要があるかもしれない。
$ blastn -query tmp/2p-7.fasta -db reference/SUDL -out tmp/2p-7.xml -outfmt 5 -max_target_seqs 1 -max_hsps 1
検索結果はこんな感じ
<Iteration>
<Iteration_iter-num>1</Iteration_iter-num>
<Iteration_query-ID>Query_1</Iteration_query-ID>
<Iteration_query-def>97b33a7e-4236-4a9c-9ceb-ff9b41540e2c runid=cde8f4c2-9722-4978-885d-aca7c6634fc5 ch=990 start_time=2026-09-16T21:44:48.358392+09:00 flow_cell_id=PBI88494 basecall_gpu=NVIDIA_A800_80GB_PCIe protocol_group_id=20260916-RUN1-PBI88494-Japan_KK sample_id= barcode=barcode89 barcode_alias=barcode89 parent_read_id=97b33a7e-4236-4a9c-9ceb-ff9b41540e2c basecall_model_version_id=dna_r10.4.1_e8.2_400bps_sup@v5.2.0</Iteration_query-def>
<Iteration_query-len>937</Iteration_query-len>
<Iteration_hits>
<Hit>
<Hit_num>1</Hit_num>
<Hit_id>gnl|BL_ORD_ID|1</Hit_id>
<Hit_def>SUDL2</Hit_def>
<Hit_accession>1</Hit_accession>
<Hit_len>5430</Hit_len>
<Hit_hsps>
<Hsp>
<Hsp_num>1</Hsp_num>
<Hsp_bit-score>1583.7</Hsp_bit-score>
<Hsp_score>857</Hsp_score>
<Hsp_evalue>0</Hsp_evalue>
<Hsp_query-from>1</Hsp_query-from>
<Hsp_query-to>937</Hsp_query-to>
<Hsp_hit-from>4761</Hsp_hit-from>
<Hsp_hit-to>3822</Hsp_hit-to>
<Hsp_query-frame>1</Hsp_query-frame>
<Hsp_hit-frame>-1</Hsp_hit-frame>
<Hsp_identity>918</Hsp_identity>
<Hsp_positive>918</Hsp_positive>
<Hsp_gaps>13</Hsp_gaps>
<Hsp_align-len>945</Hsp_align-len>
<Hsp_qseq>GAGAATTCTACATGTATCTGGTGAAACACTCCTGTAAATCAATATCAGGTCTCGCTACTGAAGAAAAACATGAATGGTTAAAAGTTGAGCCTGCAAAGAGAAAACAAGAACTAAGGTGTAATGTGCTATGAGCAGAGCTTGATCATGAAAAACATAAGAAAATTGATACTCTAAAGCAGAATTCTCTTCCTTGCTCAACACATTCACAACTCTCTTCAAGATAACTTCAATAGGCTGGTCGAGCATTTTCAACCGCTTTGGAGATGGCTTGCCATTTCTCGCCATGCTTCGGCTCACAAGCTACACACTTCGCCACAACTTCCTTCCTGTCCTCATCAGAGCCATGTTGAAGTTGAAGATTGCTTAAACAACGCCCAGAAATCGCCAATATCTGGGCGAGCGGTAACAGCACGTTCAAACCATGCTCTGGCTTTCTCCACCTTCTTGTCTTGCCAAAAGAGCTTGGCAACAGCTATGGTGACATAGA-GGTCGCGATCACACTTCTTCATAGCATCAATACTCTTTGTTTTCCGT-GAGGGACGTGGCGCCATCTCAATGTCCGCAGCCCAGAGGATACCACTCTTGGGACAGTCTTGCAGAGCCTTTGACATCAAGTGCTCTGCTTCTCTCTTGTTGTCATGCCTCAGTTCAGCGCGAATAGCAGCTAACCATAGCTCAGCCCCCC-AGGATTCTTCTTCCTGGCCGTAGT-AGA-TAGCC-GAGCTTTGTT-AGCCCCATTCACTTTCTCTTCAAGATCAGCGAGTGAGAGCCACAGTGGTATGCAGTGGGGACAGTGCTTCAAACCAGTGTCATAAGCTTTCCTGGCCTGTTCTAGATGCTTGAACCATTTGTGGAGCT-CCCAAGCATTAACCAAAGCTTGAAGAATGTTGGGAATTGCTTCAAACCTTCATTAAGCAATCTCCTCTCCTCCTCTACGTTG</Hsp_qseq>
<Hsp_hseq>GAGAATTCTACATGTATCTGGTGAAACACTCCTGTAAATCAATATCAGGTCTCGCTACTGAAGAAAAACATGAATGGTTAAAAGTTGAGCCTGCAAAGAGAAAACAAGAACTAAGGTGTAATGTGCTATGAGCAGAGCTTGATCATGAAAAACATAAGAAAATTGATACTCTCAAGCAGAATTCTCTTCCTTGCTCAACGCATTCACAACTCTCTTCAAGATAACTTCAATAGGCTGGT-GAGCATTTTCAACCGCTTTGGAGATGGCTTGCCATTTCTCGCCATGCTTCGGCTCACAAGCTACACACTTCGCCACAACTTCCTTCCTGTCCTCATCAGAGCCATGTTGAAGTTCAA-ATTTGTAGAACAACGCCCAGAAATCGCCAATATCTGGGCCAACGGTAACAGCACGTTCAAACCATGCTCTGGCTTTCTCCACCTTCTTGTCTTGCCAAAAGAGCTTGGCAACAGCTATGGTGACAT-GAGGGTCGCGATCACACTTCTTCATAGCATCAATACTCTTTGTTTTCCGTCGA-GGACGTGGCGCCATCTCAATGTCCGCAGCCCAGAGGATACCACTCTTGGGACAGTCTTGCAGAGCCTTTGACATCAAGTGCTCTGCTTCTCTCTTGTTGTCATGCCTCAGTTCAGCGCGAATAGCAGCTAACCATAGCTCAGCCCCCCCAGGATTCTTCTTCCTGGCCGTGGTGAGAATAGCCCGAGCTTTGTTTAG-CCCATTCACTTTCTCTTCAAGATCAGCGAGTGAGAGCCACAGTGGTATGCAGTGGGGACAGTGCTTCAAACCAGTGTCATAAGCTTTCCTGGCCTGTTCTAGATGCTTGAACCGTTCCTCGAGCTGCCCAAGCATTAACCAAAGCTTGAAGAATGTTGGGAATTGCTTCAAACCTTCATTAAGCAATCTCCTCTCCTCCTCTACGTTG</Hsp_hseq>
<Hsp_midline>|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||| ||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| || ||| | ||||||||||||||||||||||||||||||| | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| || ||||||||||||||||||||||||||||||||||||||||||||||| || |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ||||||||||||||||||||| || ||| ||||| |||||||||| || ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| || | ||||| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||</Hsp_midline>
</Hsp>
</Hit_hsps>
</Hit>
</Iteration_hits>
<Iteration_stat>
<Statistics>
<Statistics_db-num>4</Statistics_db-num>
<Statistics_db-len>24161</Statistics_db-len>
<Statistics_hsp-len>16</Statistics_hsp-len>
<Statistics_eff-space>22193337</Statistics_eff-space>
<Statistics_kappa>0.46</Statistics_kappa>
<Statistics_lambda>1.28</Statistics_lambda>
<Statistics_entropy>0.85</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
後半
アライメントと位置情報を取り出す
BLASTが出力したXMLの中から必要な情報を抜き出す。
リードごとにIterationがつくられ、その中にHitがある。Hit_defがリファレンス側の塩基配列の名前で、Hit_hspsがアライメントを示す。 Hsp_query-fromにBLASTでいうところのクエリーの位置情報があり、これがリードの一致した部分の始まり位置。Hsp_hit-fromがリファレンス側の位置。 アライメントは Hsp_qseqとHsp_hseqにある。
これらの情報を抜き出し、次の処理に持っていきやすい形に整える必要があるので、プログラムを作成する。 どのような言語でつくってもよいが、今回は私の得意なPHPで行う。
XMLを簡単に扱えるようにするsimplexmlを使うので、php-xmlもインストールする必要がある。
$ sudo apt install php-cli php-xml
scripts/parse_blastxml.php
#!/usr/bin/php
<?php
/* コマンドラインの引数は $argv に入っている。
* 一つ目の引数が存在し、読込可能なファイルであればxmlとして読み込む。
* 無い場合は、標準入力を読み込む
*/
if(isset($argv[1]) && file_exists($argv[1])){
$xml = simplexml_load_file($argv[1]);
}else{
$xml = simplexml_load_string(stream_get_contents(STDIN));
}
// Hitというタグを探す
foreach($xml->xpath('.//Hit') as $xhit){
// Hit_defというタグの値を$valuesにいれる
$values = ['subject' => (string) $xhit->{"Hit_def"}];
// Hspというタグを探す
foreach($xhit->xpath('.//Hsp') as $xhsp){
// Hsp_hit-fromなどの値を取り出す。先頭のHsp_は共通なので、まとめてつける
foreach(['hit-from', 'hit-to', 'query-from', 'query-to', 'qseq', 'hseq'] as $key){
$values[$key] = (string) $xhsp->{"Hsp_$key"};
}
// 値をつなげて出力する
echo join("\t", array_values($values)) . "\n";
}
}
blastはクエリーの向きにサブジェクトをそろえる。今回はリファレンスを基準にリードを並べたいので、逆向きの場合は相補鎖に変換する必要がある。そのための処理を加える。
scripts/parse_blastxml.php
<?php
/* コマンドラインの引数は $argv に入っている。
* 一つ目の引数が存在し、読込可能なファイルであればxmlとして読み込む。
* 無い場合は、標準入力を読み込む
*/
if(isset($argv[1]) && file_exists($argv[1])){
$xml = simplexml_load_file($argv[1]);
}else{
$xml = simplexml_load_string(stream_get_contents(STDIN));
}
// Hitというタグを探す
foreach($xml->xpath('.//Hit') as $xhit){
// Hit_defというタグの値を$valuesにいれる
$values = ['subject' => (string) $xhit->{"Hit_def"}];
// Hspというタグを探す
foreach($xhit->xpath('.//Hsp') as $xhsp){
// Hsp_hit-fromなどの値を取り出す。先頭のHsp_は共通なので、まとめてつける
foreach(['hit-from', 'hit-to', 'query-from', 'query-to', 'qseq', 'hseq'] as $key){
$values[$key] = (string) $xhsp->{"Hsp_$key"};
}
if(intval($values['hit-to']) < intval($values['hit-from'])){
$values['hseq'] = complement($values['hseq']);
$values['qseq'] = complement($values['qseq']);
swap($values['hit-from'], $values['hit-to']);
swap($values['query-from'], $values['query-to']);
}
// 値をつなげて出力する
echo join("\t", array_values($values)) . "\n";
}
}
/* 相補鎖変換する関数 */
function complement($seq){
$replace = [
'AaCc' => '1234',
'TtGg' => 'AaCc',
'1234' => 'TtGg',
];
foreach($replace as $from => $to){
$seq = strtr($seq, $from, $to);
}
$rev = '';
for($i = 0 ; $i < strlen($seq) ; $i ++){
$rev = $seq[$i] . $rev;
}
return $rev;
}
/* 値を入れ替える */
function swap(&$a, &$b){
$buf = $b;
$b = $a;
$a = $buf;
}
こんな感じで実行する。
$ php scripts/parse_blastxml.php tmp/2p-7.xml >tmp/2p-7.aln
リードを数える
Rで行う。
read: 1 CGCCTAAA....GTCATACA 4224
ref : 2937 CGCCTAAA....GTCATACA 1286
こういうデータから
2937 C
2938 G
2939 C
...
というデータをつくる。
難しいのは
read: 937 CAA-CGTAGAGGGAG....ATTCTC 1
ref : 3822 CAAACGTAGAGG-AG....ATTCTC 4761
こういうデータで、ギャップがあるときにどう処理するのかということ。
リード側の欠失は – という塩基に置換されたことにする。 リード側の挿入は一つ前の塩基とつなげて、G -> GG という置換として処理する。
上記の例だと、
3822 C
3823 A
3824 A
3825 -
3826 C
3827 G
3828 T
3829 A
3830 G
3831 A
3832 G
3833 GG
3834 A
と書き出すことにする。
BLASTのアライメントが挿入をどういれるのかというところまで調整できていないので、相補鎖変換したときのギャップの位置がずれている可能性がある。 今後問題になったときに考える
まずはRで読み込んでみる。そのためにRをインストール。
$ sudo apt install r-base-core
シンプルなデータを使って、計算方法を考える
$ R # Rを起動
# 挿入と欠失とミスマッチが1つずつあるアライメントを用意。1行だけの表としてデータを作成
> d1 <- data.frame('SUDL2', 1, 10, 1, 10, 'AAATTCC-GG', 'AA-TTTCCGG')
# 列名をつける
> names(d1) <- c('subject', 'hfrom', 'hto', 'qfrom', 'qto', 'hseq', 'qseq')
# データを確認
> d1
subject hfrom hto qfrom qto hseq qseq
1 SUDL2 1 10 1 10 AAATTCC-GG AA-TTTCCGG
# 一文字ずつ分割する
> strsplit(d1$hseq, '')
[[1]]
[1] "A" "A" "A" "T" "T" "C" "C" "-" "G" "G"
# リストになっているので、それをとりだす
> hseq <- strsplit(d1$hseq, '')[[1]]
> hseq
[1] "A" "A" "A" "T" "T" "C" "C" "-" "G" "G"
# リード側も同じように処理
> qseq <- strsplit(d1$qseq, '')[[1]]
> qseq
[1] "A" "A" "-" "T" "T" "T" "C" "C" "G" "G"
# 塩基の位置をリファレンスを基準につくる
# リファレンス側が-の場合、位置を増やさない方法を考える
> hseq != '-'
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
# TRUEが1、FALSEが0として、累積和を求める
> cumsum(hseq != '-')
[1] 1 2 3 4 5 6 7 7 8 9
# リファレンス側のアライメント開始位置から、塩基位置を計算する
> d1$hfrom - 1 + cumsum(hseq != '-')
[1] 1 2 3 4 5 6 7 7 8 9
> hpos <- d1$hfrom - 1 + cumsum(hseq != '-')
# 塩基の位置ごとにリード側の文字列を結合する。
> tapply(qseq, hpos, paste, collapse = '')
1 2 3 4 5 6 7 8 9
"A" "A" "-" "T" "T" "T" "CC" "G" "G"
> q <- tapply(qseq, hpos, paste, collapse = '')
# リファレンス側はギャップなしになるように-を削除する。
> hseq[hseq == '-'] <- ''
> tapply(hseq, hpos, paste, collapse = '')
1 2 3 4 5 6 7 8 9
"A" "A" "A" "T" "T" "C" "C" "G" "G"
> h <- tapply(hseq, hpos, paste, collapse = '')
# リードがカバーするリファレンスをつくった。
> q == h
1 2 3 4 5 6 7 8 9
TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE TRUE
これをリードごとに処理して、遺伝子ごとに数える
scripts/count_read.R
# Rscript script.R alignment というように処理対象を引数で渡したい。
# それを取得するための関数 commandArgs。
argv <- commandArgs(trailingOnly = T)
# タブ区切りテキストを読み込む
d1 <- read.delim(argv[1], header = F, stringsAsFactor = F)
# 列に名前をつける。これは前のPHPの処理と対応させる
names(d1) <- c('subject', 'hfrom', 'hto', 'qfrom', 'qto', 'hseq', 'qseq')
# split(d1, d1$subject) でd1のデータをsubject列、つまり遺伝子ごとに分ける。
# lapply は最初の引数のリストを個別に、2番目の引数の関数に渡す
counts <- lapply(split(d1, d1$subject), function (dr){
# d1の一つ目の遺伝子のデータがdrとして渡されている
# 今度はそれを1行ごとに処理する。rownames(dr) は行名で、今回はread.delimした時点での行番号(つまり重複がない)
d2 <- lapply(split(dr, rownames(dr)), function (ds){
# hseqを空文字''で分割する。つまり1文字ごとに分けられる。strsplitは配列に対して作用し、返り値はリストになっている。ここではds$hseqは要素が1つの配列なので、[[1]]をつけてリストでなくしている。
hseq <- strsplit(ds$hseq, '')[[1]]
qseq <- strsplit(ds$qseq, '')[[1]]
# hseqが'-'でなければ、true(1)、-ならばfalse(0)にする。これを累積和(cumsum)することで、ギャップをスキップしつつ位置情報を得る
hpos <- ds$hfrom - 1 + cumsum(hseq != '-')
# リードの配列を、位置情報毎に、paste（文字列連結)する。collapse = ''とすることで、配列を結合する ('A', 'C') -> 'AC' になる。
q <- tapply(qseq, hpos, paste, collapse = '')
hseq[hseq == '-'] <- ''
h <- tapply(hseq, hpos, paste, collapse = '')
q == h
})
pos <- unlist(lapply(d2, names))
match <- unlist(d2, use.names = F)
res <- cbind(total = table(pos), match = tapply(match, pos, sum))
res[order(as.numeric(rownames(res))),]
})
このデータを次の処理にどうもっていくのかを考える。
- この処理の続きとしてグラフを書く。
- メリット: 中間ファイルの形式を考えなくてよい。
- デメリット: たぶんこれはそれなりに重い(時間がかかる)処理なので、プログラムの検証に時間がかかる。
- Rのオブジェクトとして書き出す。
- メリット: リストなど複雑な形式のデータがそのままもっていける。
- デメリット: 作られたファイルの中身を知ることができない。
- テキストとして書き出す。
- メリット: 計算結果を確認できる。
- デメリット: 複雑なデータを書き出すと、読込の処理が面倒。
今回はRのオブジェクトのまま保存することにする。 拡張子、Robjの保存先をスクリプトの2番目の引数として渡す。
scripts/count_read.R
# Rscript script.R alignment というように処理対象を引数で渡したい。
# それを取得するための関数 commandArgs。
argv <- commandArgs(trailingOnly = T)
# タブ区切りテキストを読み込む
d1 <- read.delim(argv[1], header = F, stringsAsFactor = F)
# 列に名前をつける。これは前のPHPの処理と対応させる
names(d1) <- c('subject', 'hfrom', 'hto', 'qfrom', 'qto', 'hseq', 'qseq')
# split(d1, d1$subject) でd1のデータをsubject列、つまり遺伝子ごとに分ける。
# lapply は最初の引数のリストを個別に、2番目の引数の関数に渡す
counts <- lapply(split(d1, d1$subject), function (dr){
# d1の一つ目の遺伝子のデータがdrとして渡されている
# 今度はそれを1行ごとに処理する。rownames(dr) は行名で、今回はread.delimした時点での行番号(つまり重複がない)
d2 <- lapply(split(dr, rownames(dr)), function (ds){
# hseqを空文字''で分割する。つまり1文字ごとに分けられる。strsplitは配列に対して作用し、返り値はリストになっている。ここではds$hseqは要素が1つの配列なので、[[1]]をつけてリストでなくしている。
hseq <- strsplit(ds$hseq, '')[[1]]
qseq <- strsplit(ds$qseq, '')[[1]]
# hseqが'-'でなければ、true(1)、-ならばfalse(0)にする。これを累積和(cumsum)することで、ギャップをスキップしつつ位置情報を得る
hpos <- ds$hfrom - 1 + cumsum(hseq != '-')
# リードの配列を、位置情報毎に、paste（文字列連結)する。collapse = ''とすることで、配列を結合する ('A', 'C') -> 'AC' になる。
q <- tapply(qseq, hpos, paste, collapse = '')
hseq[hseq == '-'] <- ''
h <- tapply(hseq, hpos, paste, collapse = '')
q == h
})
pos <- unlist(lapply(d2, names))
match <- unlist(d2, use.names = F)
res <- cbind(total = table(pos), match = tapply(match, pos, sum))
res[order(as.numeric(rownames(res))),]
})
# 保存する
save('counts', file = argv[2])
こんな感じで実行する。
$ Rscript scripts/count_read.R tmp/2p-7.aln tmp/2p-7.Robj
図にする
上記で塩基ごとにリード数とミスマッチを数えることができた。これを図にする。 遺伝子毎に図にしつつ、縮尺をそろえたい。
縮尺を揃えるということは、遺伝子の長さの情報が必要で、それはリファレンスから読み込むことにする。
# リファレンス(FASTA)を読み込む
readFasta <- function (file) {
buf <- scan(file, what = 'c')
ref <- sapply(split(buf, cumsum(substr(buf, 1, 1) == '>')), function (single) {
name <- strsplit(single[1], '[[:space:]]')[[1]][1]
name <- substring(name, 2)
seq <- paste(single[-1], collapse = '')
c('name' = name, 'seq' = seq)
})
data.frame(t(ref))
}
出力するファイル名を引数で与える必要がある。
lapply(names(counts), function (gene){
png(paste0(base_name, gene, ".png"))
})
というところを考えると、引数にはreference.faと、画像のファイル名が必要ということになる。
scripts/pile.R
# 引数 Rojb reference.fa image_path (.gene_name.pngを加える)
argv <- commandArgs(trailingOnly = T)
load(argv[1]) # countsが読み込まれる
# リファレンス(FASTA)を読み込む
readFasta <- function (file) {
buf <- scan(file, what = 'c')
ref <- sapply(split(buf, cumsum(substr(buf, 1, 1) == '>')), function (single) {
name <- strsplit(single[1], '[[:space:]]')[[1]][1]
name <- substring(name, 2)
seq <- paste(single[-1], collapse = '')
c('name' = name, 'seq' = seq)
})
data.frame(t(ref))
}
ref <- readFasta(argv[2])
max_len <- max(nchar(ref$seq))
max_height <- max(sapply(counts, function (count){
max(count[,'total'])
}))
lapply(names(counts), function (gene){
png(paste0(argv[3], '.', gene, '.png'), width = 1024, height = 300)
count <- counts[[gene]]
plot(as.numeric(rownames(count)), count[,'total'],
xlim = c(1, max_len), ylim = c(1, max_height),
type = 'h', col = 'green',
xlab = 'position', ylab = 'count', main = gene
)
points(as.numeric(rownames(count)), count[,'total'] - count[,'match'], type = 'h', col = 'red')
dev.off()
})
次のように実行する。
$ Rscript scripts/pile.R tmp/reference.fa tmp/2p-7.Robj
レポートを作成
図を一覧表示するindex.htmlを作成する。
scripts/make_html.php
<?php
/* コマンドラインの引数は $argv に入っている。
* 一つ目の引数が存在し、読込可能なファイルであればxmlとして読み込む。
* 無い場合は、標準入力を読み込む
*/
if(! isset($argv[1])) die("引数がありません\n");
$images = glob(str_replace('.png', '.*.png', $argv[1]));
if(count($images) === 0) die("画像ファイルが見つかりません: {$argv[1]}\n");
?>
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>BLAST結果</title>
<style>
body { font-family: monospace; }
.blast-image { max-width: 100%; height: auto; }
</style>
</head>
<body>
<?php
foreach($images as $image){
$gene = explode('.', basename($image, '.png'));
$gene = array_pop($gene);
?>
<div>
<h3>Gene: <?php echo htmlspecialchars($gene); ?></h3>
<img src="<?php echo htmlspecialchars($image); ?>">
</div>
<?php
}
?>
</body>
</html>
次のように実行する。
$ php scripts/make_html.php tmp/2p-7.png >2p-7.html
Dockerのコンテナにする
広く一般の利用者に使ってもらうためにDockerのイメージを作成する。 Dockerのイメージをつくることによって、blastなどをインストールする必要がなくなる。
手順を一般化する
固有のファイル名をなくして、連続的にプログラムが実行されるように make というしくみを使う。
今回行ったのは下記の手順で、前のファイルから次のファイルをつくるというようになっている。 この依存関係をMakefileに書き出す。
- リファレンスからmakeblastdb
- fastq.gzからfasta
- fastaからblast_xml
- blast_xmlからaln
- alnからRobj (カウントデータ)
- Robjからpng
- pngからindex.html
Makefile
SHELL = /bin/bash
WORK_DIR = tmp
SCRIPT_DIR ?= scripts/
.SECONDARY:
FASTQ = $(wildcard raw_fastq/*.fastq.gz)
HTML = $(patsubst raw_fastq/%.fastq.gz, %.html, $(FASTQ))
REF = $(WORK_DIR)/reference
all: $(HTML)
clean:
rm -rf $(WORK_DIR) $(HTML)
download-test:
@mkdir -p $(dir $@)
curl -o $@ https://kiku3.tsbio.info/~takamasa/260917Honda/ref/SUDL.fa
curl -O https://kiku3.tsbio.info/~takamasa/260917Honda/60-1380693045_FASTQ.zip
unzip 60-1380693045_FASTQ.zip
rm 60-1380693045_FASTQ.zip
%.html: $(SCRIPT_DIR)/make_html.php $(WORK_DIR)/%.png
php $^ >$@
$(WORK_DIR)/%.png: $(SCRIPT_DIR)/pile.R $(WORK_DIR)/%.Robj $(REF).fa
Rscript $^ $(basename $@) && touch $@
$(WORK_DIR)/%.Robj: $(SCRIPT_DIR)/count_read.R $(WORK_DIR)/%.aln
Rscript $^ $@
$(WORK_DIR)/%.aln: $(SCRIPT_DIR)/parse_blastxml.php $(WORK_DIR)/%.xml
php $^ >$@
$(WORK_DIR)/%.xml: $(WORK_DIR)/%.fasta $(REF).nin
blastn -query $< -db $(REF) -out $@ -outfmt 5 -max_target_seqs 1 -max_hsps 1
$(WORK_DIR)/%.fasta: raw_fastq/%.fastq.gz
@mkdir -p $(dir $@)
zcat $< | sed -e "1~4s/^@/>/;3~4d;4~4d" >$@
$(REF).nin: $(REF).fa
makeblastdb -in $< -out $(REF) -dbtype nucl
$(REF).fa: $(wildcard reference/*.fa) $(wildcard reference/*.fasta)
@mkdir -p $(dir $@)
cat $^ >$@
DOCKER_USER = takamasa
IMAGE_NAME = analyze-ez-pcr
VERSION = 1.0
docker-build:
docker build \
-t ${DOCKER_USER}/${IMAGE_NAME}:${VERSION} \
-t ${DOCKER_USER}/${IMAGE_NAME}:latest .
docker-push:
docker push ${DOCKER_USER}/${IMAGE_NAME}:${VERSION}
docker push ${DOCKER_USER}/${IMAGE_NAME}:latest
こんな感じで実行する。
$ make 2p-7.html
Dockerのイメージをつくる
Dockerfileをつくる。 これは既存のイメージからカスタマイズしたイメージをつくる手順を書いたもの。
一連の処理を行うのに必要なパッケージをインストールし、作成したスクリプトをコピーする。
Dockerfile
FROM ubuntu:22.04
# 環境変数の設定（対話型プロンプトを無効化）
ENV DEBIAN_FRONTEND=noninteractive
# 必須パッケージのインストール
RUN apt-get update && apt-get install -y --no-install-recommends \
make curl gzip unzip \
ncbi-blast+ \
php-cli php-xml libxml2-dev \
r-base-core r-base-dev \
libcurl4-openssl-dev libssl-dev \
&& rm -rf /var/lib/apt/lists/*
# Rの追加パッケージが必要な場合はここで指定（例: ggplot2 や Bioconductor等）
# RUN Rscript -e "install.packages(c('ggplot2', 'dplyr'), repos='https://cloud.r-project.org/')"
WORKDIR /app
COPY scripts/ /app/scripts/
COPY Makefile /app/Makefile
ENV SCRIPT_DIR=/app/scripts
# 作業ディレクトリの設定
WORKDIR /workspace
# コンテナ起動時に自動で make を実行する場合の設定
CMD ["make", "-f", "/app/Makefile"]
イメージを作成
$ docker build -t mypipeline
実行
$ docker run --rm -v .:/workspace mypipeline