muddy notes & thoughts
核酸医薬ヌシネルセンが相補鎖を構成し得る核酸を、ヒトゲノム塩基配列の中から検索することを試みます。
( 注 2019.01.12 bowtie_queryEditに、とんでもない初歩的な間違いが含まれていた旨の指摘をいただきました。いかに残念な間違いだったか、念のために(?)残しておきます )。
但し書き:
リファレンス・ゲノム内で見つかったヒットの数と、PMDAの審査資料 17ページで言及されている遺伝子の領域にヒットしたかをまとめると、下記の表のようになります。
BNC2 blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
EFCAB6 blast - bowtie2 bowtie_queryEdit gggenome last vmatch
FIGLA blast bowtie - bowtie_queryEdit gggenome last vmatch
FOXP1 blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
FREM2 - bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
MECOM blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
MSR1 blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
RPGRIP1L blast bowtie - bowtie_queryEdit gggenome last vmatch
RSF1 blast - bowtie2 bowtie_queryEdit gggenome last vmatch
STAT4 blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
TRHDE blast bowtie bowtie2 bowtie_queryEdit gggenome last vmatch
ULK2 blast - - bowtie_queryEdit gggenome last vmatch
hg38ヒット数 453 231 394 774 774 660 774
それぞれのヒット位置については、ゲノム軸上で閲覧できるようにしてあります ( BNC2 / EFCAB6 / FIGLA / FOXP1 / FREM2 / MECOM / MSR1 / RPGRIP1L / RSF1 / STAT4 / TRHDE / ULK2。テーブル・ブラウザ を経由して、ヒット位置をダウンロードすることもできます )。
blast, bowtie, gggenome, last, vmatchはいずれもアラインメント・ツールの名称です(詳細などは次節参照)。bowtie_queryEditは、クエリ配列に1塩基、2塩基の挿入・欠失を あらかじめ入れた上で bowtieを用い検索した結果です(bowtieはギャップを含むアラインメントを見つけられるようには設計されていないため)。
但し書き:
curl https://gggenome.dbcls.jp/hg38/2/TCACTTTCATAATGCTGG.bed \
| grep -v '_' \
| grep -v track
さすがに便利です。そして、早い。
blastn -db hg38.fa -query query.fa -outfmt "17 SR SQ" -word_size 4 -strand both -evalue 100000 \
| samtools calmd - hg38.fa \
| awk 'BEGIN{OFS="\t"}{
if (match($0,"^@")){
print
}else{
edit_dist = 0;
if ( match($0,"NM:i:[0-9]+") ) { edit_dist += substr( $0, RSTART + 5, RLENGTH - 5) }
if ( match($6, "^[0-9]+H" ) ) { edit_dist += substr( $6, RSTART, RLENGTH - 1) }
if ( match($6, "[0-9]+H$" ) ) { edit_dist += substr( $6, RSTART, RLENGTH - 1) }
$5 = edit_dist # put edit_dist as score
$1 = $1 ";edit_dist:"edit_dist ";CIGAR:" $6
if ($5 <= 2){ print }
} }' \
| samtools view -hu - \
| bamToBed
やっていることは、次のとおりです。
（久しぶりに使ってみたところ、なんとsamフォーマットでアラインメントが出力できるようになっていて、びっくりしました。blastもがんばっていますね。ただしmakeblastdb の際に -parse_seqids をつけないと、染色体名がきちんと認識されないので注意を。）
bowtie --maxbts 100000 -f --all -v 2 --sam hg38 query.fa \
| samtools view -hb - \
| bamToBed
ミスマッチを2個まで許し(-v 2)、みつかったものを全部出力(–all)しています (ins/delが探せない以外は、悪くない)。
bowtie2 --all -x hg38 -D 40 -R 3 -N 0 -L 6 -i S,1,0 --gbar 1 --score-min G,-5,-5 -f query.fa \
| grep -e @ -e NM:i:0 -e NM:i:1 -e NM:i:2 \
| samtools view -hb - \
| bamToBed
seed長を6(-L 6)にして、1bpごとにスライディングさせています(-i S,1,0)。低いアラインメントスコアでもとにかく出力するように (–score-min G,-5,-5)。出力結果から、NMタグを使ってedit distance2以下のものだけを選択しています (期待した程は、見つけてくれなかったですね)。
query=">query\nTCACTTTCATAATGCTGG"
delOne ()
{
awk '{
if ( $1 ~ "^>" ) { orig_name = substr($1,2); next}
len = length($1);
for (i=0;i<len;i++)
{
name = orig_name ";" substr($1,1,i) "_" substr($1,i+2)
str = substr($1,1,i) substr($1,i+2)
print ">" name "\n" str
}
}'
}
insOne ()
{
awk '{
if ( $1 ~ "^>" ) { orig_name = substr($1,2); next}
split( "A T G C" , nuc, " " )
len = length($1);
for (n in nuc)
{
for (i=0;i<=len;i++)
{
name = orig_name ";" substr($1,1,i) nuc[n] substr($1,i+1)
str = substr($1,1,i) nuc[n] substr($1,i+1)
print ">" name "\n" str
}
}
}'
}
indelOne ()
{
echo $1 | insOne
echo $1 | delOne
}
indelTwo ()
{
echo $1 | insOne | insOne
echo $1 | delOne | delOne
echo $1 | delOne | insOne
}
echo $query | bowtie --maxbts 100000 -f --all -v 2 --sam hg38 - > zero.sam
indelOne $query | bowtie --maxbts 100000 -f --all -v 1 --sam hg38 - > one.sam
indelTwo $query | bowtie --maxbts 100000 -f --all -v 0 --sam hg38 - > two.sam
samtools merge - zero.sam one.sam two.sam \
| bamToBed
無理矢理やなぁ…という力技(brute force)ですが、この程度の配列長とins/del数であれば、あっという間に終わります。下手に凝ったツールやパラメータを探したりするよりもよっぽど直球かもしれません(配列長がもっと長くなれば、bowtie2やlastとかが真価を発揮するのでしょうし)。ins/del無しのときは編集距離2以内の結果すべてを出力 (–all -v 2)、ins/delが1個のときは編集距離1以内(–all -v 1)、ins/delが2個のときは編集距離0の結果すべてを出力 (–all -v 0)しています。
lastal -L6 -m0 -d0 -y0 -a0 -x 2 -e 14 hg38-m1R00 query.fa \
| maf-convert sam /dev/stdin \
| tr '=' 'M' \
| samtools view -h -T hg38.fa - \
| samtools calmd - hg38.fa \
| awk 'BEGIN{OFS="\t"}{
if (match($0,"^@")){
print
}else{
edit_dist = 0;
if ( match($0,"NM:i:[0-9]+") ) { edit_dist += substr( $0, RSTART + 5, RLENGTH - 5) }
if ( match($6, "^[0-9]+H" ) ) { edit_dist += substr( $6, RSTART, RLENGTH - 1) }
if ( match($6, "[0-9]+H$" ) ) { edit_dist += substr( $6, RSTART, RLENGTH - 1) }
$5 = edit_dist # put edit_dist as score
$1 = $1 ";edit_dist:"edit_dist ";CIGAR:" $6
if ($5 <= 2){ print }
} }' \
| samtools view -hu - \
| bamToBed
lastはとてもいろいろできるツールですが、こちらのドキュメントを参考にしてパラメータを設定しました(“3. Guaranteeing to align reads with up to N mismatches or indels”)。インデックスにはlastdbを用い、可能な限り多くのseedをとり(-m1)、繰り返し配列をマスクしない(-R00)ように構築しました。編集距離については上記blastと同様、検索終了後にsamtoolsを使って計算しました。（でも、もっと設定を工夫すればきちんと見つかるかも。。。時間のあるときに、開発者の方に聞いてみたいと思います）
vmatch -d -p -e 2 -showdesc 10 -q query.fa -complete hg38.fa \
| grep -v ^# \
| awk 'BEGIN{OFS="\t"}{
left_len=$1; left_seq=$2; left_start=$3;
strand = "+"; if ( $4 == "P") {strand = "-"}
right_len=$5; right_seq=$6; right_start=$7;
dist=$8; eval=$9; score=$10;
name = right_seq ",dist:"dist ",eval:" eval
print left_seq, left_start, left_start + left_len, name, score, strand
}'
編集距離が2以下のもの(-e 2)、というオプションを指定しています (さすがvmatch、の結果でしたね)。
tags: japanese