複数のPCR産物をまとめてシークエンス
クイックスタート
解析データを用意
- シークエンスデータ(.fastq.gz)をraw_fastqディレクトリ直下に置く。
- リファレンス配列(.fa)をreferenceディレクトリ直下に置く。
- ファイルの配置はテストデータを参照してください。
Dockerを使用できる場合 (おすすめ)
$ docker run --rm -v ./:/workspace tsbioinfo/multilocus-seq demo # 動作確認用データダウンロード $ docker run --rm -v ./:/workspace tsbioinfo/multilocus-seq
- tsbioinfo/multilocus-seq というdocker imageを使用します。
- 解析結果を出力するフォルダを/workspace にマウントします。
- 引数なしで実行すると、raw_fastqフォルダにある全てのfastq.gzを解析します。
Dockerを使用できない場合
blast、R、PHPなどが動作する環境が必要です。 必要なパッケージはリポジトリに含まれるDockerfileを参照してください。
$ git clone https://kiku3.tsbio.info/git/multilocus-seq.git mls $ cd mls $ make demo # 動作確認用データダウンロード $ make
- gitリポジトリをクローンします。
- raw_fastq と reference ディレクトリをつくって、所定のファイルを置きます。
- make で解析します。
解析結果
ソフトウェア開発の背景 (これは何をするものか)
変異株をスクリーニングして、候補株が得られましたが、既知の変異株のアリルかもしれません。 そこで、候補遺伝子だけをシークエンスしたいと考えました。 1-2 kbpの領域に分けて、PCRして、シークエンス、となるわけですが、まあまあのコストがかかります。
そんなとき、Azenta Life Sciencesという会社のPCR-EZというサービスを知りました。 「単一なクローンで直鎖状」PCR産物を、(おそらく)Nanoporeを使ってシークエンスしてくれるものです。 1サンプルあたり4000円(2026年10月現在)と、サンガーシークエンスの10回分にも満たない価格です。
このサービスに、いくつかのPCR産物を混ぜたものを1サンプルとして送付したところ、シークエンスはちゃんとできるのですが、参照配列との比較はしてもらえませんでした(当たり前です)。
そこで、生のリード(fastq)を解析して、複数の遺伝子の変異の有無を調べることにしました。 うまくできたので、公開することにしました。
この方法の良いところは、ヘテロの挿入・欠失が同定できるということです。 サンガー法のシークエンスではぐちゃぐちゃになって読めないところが、次世代シーケンサーは1分子ごとに読むので、ヘテロもきれいに調べることができます。
Nanoporeシーケンサーについて
(簡単な説明をいずれ書く)
データ処理の流れ
- リードをリファレンスと比較して、アライメントを作成する。
- リファレンスの塩基ごとに、リードを数えてミスマッチを探す。
アライメントまで
azentaのデータはfastqファイル。nanoporeなので、かなり長い。平均で1000bpぐらい。リード数はサンプル当たり600弱。 bowtieは1000bpを超えるリードを扱えないので、BLASTを使うことにする。
BLASTは入力の形式がfastaなので、fastqからfastaに変換する必要がある。
BLASTの実行に必要なものは、リードとリファレンス。リファレンスは事前にmakeblastdbでインデックスをつくる必要がある。
数え上げ
BLASTが作成したアライメントからリードの重なりを数え上げる。 この部分は自分でスクリプトを書く必要がある。 処理がしやすいように、BLASTはXMLで出力することにする。
最後に図をつくることを考えるとRで数え上げをすることにする。 RでXMLのデータを扱うのは(自分には)難しいので、PHPで橋渡しをする。
BLAST (XML) -> PHP (タブ区切りテキスト) -> R という感じにする。
解析結果の提示
図をいれたレポートをつくりたいので、HTMLで出力することになる。 Rで描いた図と、数え上げの結果をPHPでHTMLの整える。
BLAST
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+
インデックスをつくる
上記のテストデータをダウンロードして、tar xvzfで展開する。 referenceディレクトリにSUDL.faがあり、これをリファレンスとする。
$ curl -O 'https://biochemistry.isc.chubu.ac.jp/labo/suzuki/wp-content/uploads/sysdev-01/multilocus-seq-demo.tar.gz' $ tar xvzf multilocus-seq-demo.tar.gz $ 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へ
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
<?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));
}
/* xpathをつかって必要な要素を抽出する
*/
// 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', 'hseq', 'qseq'] as $key){
$values[$key] = (string) $xhsp->{"Hsp_$key"};
}
/* BLASTはリード側の方向に合わせてリファレンスを相補鎖変換している。
* リファレンスを基準とするために、リードを相補鎖にして、位置を入れ替える
*/
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で行う。
BLASTのこういうデータから
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
これをリードごとに処理して、遺伝子ごとに数える。最後にRのオブジェクトのまま保存する(saveRDS)。
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行ごとに処理する。nrow(dr) は行数で、seq_lenで1:nrow(dr)と同じになる
d2 <- lapply(seq_len(nrow(dr)), function (i){
# hseqを空文字''で分割する。つまり1文字ごとに分けられる。
# strsplitは配列に対して作用し、返り値はリストになっている。
# ここではds$hseqは要素が1つの配列なので、[[1]]をつけてリストでなくしている。
hseq <- strsplit(dr$hseq[i], '')[[1]]
qseq <- strsplit(dr$qseq[i], '')[[1]]
# hseqが'-'でなければ、true(1)、-ならばfalse(0)にする。
# これを累積和(cumsum)することで、ギャップをスキップしつつ位置情報を得る
hpos <- dr$hfrom[i] - 1 + cumsum(hseq != '-')
# リードの配列を、位置情報毎に、paste（文字列連結)する。
# collapse = ''とすることで、配列を結合する ('A', 'C') -> 'AC' になる。
q <- tapply(qseq, hpos, paste, collapse = '')
hseq[hseq == '-'] <- ''
h <- tapply(hseq, hpos, paste, collapse = '')
data.frame(
pos = as.numeric(names(h)),
ref = as.vector(h),
read = as.vector(q),
stringsAsFactors = F
) # 返り値はデータフレーム
})
# 全データフレームを一つに結合
d3 <- do.call(rbind, d2)
# 位置ごとに数える
d4 <- as.data.frame(
table(d3$pos, d3$ref, d3$read),
stringsAsFactors = FALSE
)
names(d4) <- c("pos", "ref", "read", "count")
# 0を除いて、形式を整える
d4$count <- as.numeric(d4$count)
d4 <- d4[d4$count > 0, ]
d4$pos <- as.numeric(d4$pos)
d4 <- d4[order(d4$pos), ]
rownames(d4) <- NULL
return (d4)
})
# 保存する
saveRDS(counts, file = argv[2])
こんな感じで実行する。
$ Rscript scripts/count_read.R tmp/2p-7.aln tmp/2p-7.Robj
CSVで出力する
数え上げたデータをPHPで読めるように、CSVを出力する。
scripts/obj2csv.R
# 引数 Rojb reference.fa image_path (.gene_name.pngを加える)
argv <- commandArgs(trailingOnly = T)
counts <- readRDS(argv[1]) # データ読み込み
s <- Map(function (mat, gene){
data.frame(gene = gene, mat)
}, counts, names(counts))
d <- do.call(rbind, s)
write.csv(d, , file = "", row.names = FALSE, quote = FALSE)
図にする
上記で塩基ごとにリード数とミスマッチを数えることができた。これを図にする。 遺伝子毎に図にしつつ、縮尺をそろえたい。
縮尺を揃えるということは、遺伝子の長さの情報が必要で、それはリファレンスから読み込むことにする。
# リファレンス(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)
counts <- readRDS(argv[1]) # データ読み込み
fastaFile <- argv[2]
# 一致、不一致を数える
counts <- lapply(counts, function (count){
# refとreadが一致しているかどうかを判定判定する
# 全てTRUEのときでもFALSEが0となるように、factorに変換する
is_match <- factor(
count$ref == count$read,
levels = c(TRUE, FALSE),
labels = c('match', 'mismatch')
)
# 位置ごと、一致、不一致ごとに数える
tmp1 <- tapply(count$count, list(pos = count$pos, is_match), sum)
tmp1[is.na(tmp1)] <- 0
# データフレームにして返す
tmp2 <- data.frame(
pos = as.numeric(rownames(tmp1)),
mismatch = tmp1[,'mismatch'],
total = rowSums(tmp1),
stringsAsFactors = FALSE
)
tmp2[order(tmp2$pos), ]
})
# リファレンス(FASTA)を読み込む
readFasta <- function (file) {
buf <- scan(file, what = 'c')
# '>'で始まる行を先頭として、配列ごとにFASTAの行をまとめる
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(fastaFile)
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]]
x <- as.numeric(rownames(count))
# 描画範囲をつくる
plot(x, count[,'total'],
xlim = c(1, max_len), ylim = c(1, max_height),
type = 'n',
xlab = 'position', ylab = 'count', main = gene
)
# 全数で描画する
polygon(c(head(x, 1), rep(x, each = 2), tail(x, 1)),
c(0, rep(count[,'total'], each = 2), 0),
border = NA, col = '#9ACD32',
)
# ミスマッチを描画する
polygon(c(head(x, 1), rep(x, each = 2), tail(x, 1)),
c(0, rep(count[,'mismatch'], each = 2), 0),
border = NA, col = 'red',
)
dev.off()
})
次のように実行する。
$ Rscript scripts/pile.R tmp/reference.fa tmp/2p-7.Robj
レポートを作成
図とCSVをまとめてHTMLにする。
scripts/make_html.php
<?php
/* コマンドラインの引数は $argv に入っている。
* 一つ目の引数は画像のベースネーム
* 二つ目の引数は遺伝子、位置、ごとの塩基数
*/
if(! isset($argv[1])) die("引数がありません\n");
$images = glob(str_replace('.png', '.*.png', $argv[1]));
if(count($images) === 0) die("画像ファイルが見つかりません: {$argv[1]}\n");
$csvFile = $argv[2];
if(! file_exists($csvFile)) die("CSVファイルが見つかりません: {$csvFile}\n");
/* 関数
*/
function mismatchString($v){
return sprintf("%s(%d)", $v[3], $v[4]);
}
/* データ処理
*/
// CSVを読み込む
$lines = explode("\n", trim(file_get_contents($csvFile)));
if(count($lines) > 1){
array_shift($lines); // ヘッダ行を削除
$prev = "";
$combined = [];
foreach($lines as $line){
$v = explode(',', $line); // gene pos ref read count
$v[4] = intval($v[4]);
$mis = ($v[2] != $v[3]);
$gene = $v[0];
if($prev == $v[1]){
$cur = &$combined[$gene][count($combined[$gene])-1];
// カバー数を増やす
$cur['cover'] += $v[4];
// ミスマッチなら、ミスマッチの塩基を追加する
if($mis){
$cur['mismatches'][] = $v;
}else{
$cur['match'] += $v[4];
}
}else{
$v = array_merge($v, [
'cover' => $v[4],
'match' => $mis ? 0 : $v[4],
'mismatches' => $mis ? [$v] : []
]);
if(isset($combined[$gene])){
$combined[$gene][] = $v;
}else{
$combined[$gene] = [$v];
}
};
$prev = $v[1];
}
}
/* HTML生成
*/
?>
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>MultiLocus-Seq</title>
<style>
body { font-family: monospace; }
.blast-image { max-width: 100%; height: auto; }
table {
width: 100%;
border-collapse: collapse;
table-layout: fixed; /* 列幅を固定化 */
tr.high-mismatch {
background-color: #FFCCCC; /* 赤系の背景色 */
}
tr.low-mismatch {
background-color: white; /* 緑系の背景色 */
}
/* thead と tbody をブロック化 */
thead, tbody tr {
display: table;
width: 100%;
table-layout: fixed;
}
thead {
th:nth-child(1) {
width: 70px;
}
th:nth-child(2),
th:nth-child(3),
th:nth-child(4){
width: 50px;
}
}
tbody {
/* tbody のみスクロール領域を設定 */
display: block;
max-height: 200px;
overflow-y: auto;
td:nth-child(1) {
width: 70px;
}
td:nth-child(2),
td:nth-child(3),
td:nth-child(4) {
width: 50px;
}
td:nth-child(1),
td:nth-child(2),
td:nth-child(3) {
text-align: right;
}
}
th, td {
padding: 8px;
border: 1px solid #ddd;
text-align: left;
}
}
</style>
<script>
function filterTable(cont, tableId) {
const start = Number(cont.querySelector('.start').value);
const end = Number(cont.querySelector('.end').value);
const table = document.getElementById(tableId);
const rows = table.querySelectorAll('tbody tr');
rows.forEach(row => {
const bp = parseInt(row.cells[0].textContent, 10);
if (!isNaN(bp) && (bp < start || bp > end)) {
row.style.display = 'none';
} else {
row.style.display = '';
}
});
}
</script>
</head>
<body>
<?php
foreach($images as $i => $image){
$gene = explode('.', basename($image, '.png'));
$gene = array_pop($gene);
$tbody_html = '';
$bp_min = $bp_max = 0;
foreach($combined[$gene] as $line){
$tbody_html .= sprintf("<tr class='%s'>%s</tr>\n",
$line['match'] / $line['cover'] < 0.3 ? 'high-mismatch' : 'low-mismatch',
join("", [
sprintf("<td>%s</td>", $line[1]),
sprintf("<td>%s</td>", $line['cover']),
sprintf("<td>%s</td>", $line['match']),
sprintf("<td>%s</td>", $line[2]),
sprintf("<td>%s</td>", implode(', ', array_map(
fn($v) => sprintf("%s(%d)", $v[3], $v[4]),
$line['mismatches']
))),
])
);
$bp_min = $bp_min == 0 ? $line[1] : min($bp_min, $line[1]);
$bp_max = max($bp_max, $line[1]);
}
?>
<section>
<h2>Gene: <?php echo $gene; ?></h2>
<div>
<img src="<?php echo $image; ?>">
</div>
<details>
<summary>塩基一覧</summary>
<div class="table-wrap">
<form onsubmit="return false;"
oninput="filterTable(this, '<?php echo "table-$i"; ?>')"
>
<label>絞り込み:
<input type="number" class="start" step="1"
min="<?php echo $bp_min; ?>" max="<?php echo $bp_max; ?>"
value="<?php echo $bp_min; ?>" placeholder="<?php echo $bp_min; ?>"
> -
<input type="number" class="end" step="1"
min="<?php echo $bp_min; ?>" max="<?php echo $bp_max; ?>"
value="<?php echo $bp_max; ?>" placeholder="<?php echo $bp_max; ?>"
>
<button type="reset"
onclick="setTimeout(() => this.form.dispatchEvent(
new Event('input', { bubbles: true })
), 0)"
>リセット</button>
</label>
</form>
<table id="<?php echo "table-$i"; ?>">
<thead>
<tr><th>bp</th><th>Cover</th><th>Match</th><th>Ref</th><th>Mismatch</th></tr>
</thead>
<tbody>
<?php echo $tbody_html;?>
</tbody>
</table>
</div>
<details>
</section>
<?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からrds (カウントデータ)
- rdsからcsv
- rdsからpng
- pngとcsvから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) %.html: $(SCRIPT_DIR)/make_html.php $(WORK_DIR)/%.png $(WORK_DIR)/%.csv php $^ >$@ $(WORK_DIR)/%.csv: $(SCRIPT_DIR)/obj2csv.R $(WORK_DIR)/%.rds Rscript $^ >$@ $(WORK_DIR)/%.png: $(SCRIPT_DIR)/pile.R $(WORK_DIR)/%.rds $(REF).fa Rscript $^ $(basename $@) && touch $@ $(WORK_DIR)/%.rds: $(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 $^ >$@ DEMO_DATA = https://biochemistry.isc.chubu.ac.jp/labo/suzuki/wp-content/uploads/sysdev-01/multilocus-seq-demo.tar.gz demo: @echo "==> Downloading demo dataset..." curl -O "$(DEMO_DATA)" tar xzf $(notdir $(DEMO_DATA)) @echo "" DOCKER_USER = tsbioinfo IMAGE_NAME = multilocus-seq VERSION = 1.0 -include .env export build: docker build \ -t $(DOCKER_USER)/$(IMAGE_NAME):$(VERSION) \ -t $(DOCKER_USER)/$(IMAGE_NAME):latest . push: login docker push -a $(DOCKER_USER)/$(IMAGE_NAME) login: @echo "$(DOCKER_TOKEN)" | docker login -u $(DOCKER_USER) --password-stdin
こんな感じで実行する。
$ 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 \
libuv1-dev pkg-config pandoc \
&& rm -rf /var/lib/apt/lists/*
WORKDIR /app
COPY scripts/ /app/scripts/
COPY Makefile /app/Makefile
ENV SCRIPT_DIR=/app/scripts
# 作業ディレクトリの設定
WORKDIR /workspace
# コンテナ起動時に自動で make を実行する場合の設定
ENTRYPOINT ["make", "-f", "/app/Makefile"]
CMD []
イメージを作成
$ docker build -t tsbioinfo/multilocus-seq
実行
$ docker run --rm -v .:/workspace tsbioinfo/multilocus-seq