Written by bonohu in misc on 日 28 5月 2017.
2017年国立遺伝学研究所 国際シンポジウムとして、DDBJ 30周年記念シンポジウムが三島にて。
Janet ThorntonさんはBritish Airwaysのシステム障害のため、飛行機が飛ばず来れずにキャンセルとのことで残念だったが、Peer Borkさんのメタゲノムの話は大変参考になった。また、ヒトゲノム計画を牽引されて来た榊佳之先生の話は、現在の研究においても示唆を含む内容だった。
[「初期の自動車が馬のように上手に走れないから価値がないと言っているようなものだ。このレフェリーの評価の間違いは歴史が示すであろう」](https://twitter.com/bonohu/status/868691804811440128)
[研究費を獲得するための応援の「外圧」と判断した](https://twitter.com/bonohu/status/868695024136892417)
[“作成されたデータについては作成から24時間を基本として全て公開して全ての研究者が自由に利用できるようにするという項目を含む、バミューダ原則”](https://twitter.com/bonohu/status/868695964801875968)
とくに、みんなと画一的である必要はない、いろんなタイプの人が居ていい、という高校生へのメッセージは響いて欲しいな、と。また、我らが統合DBの高木利久先生による話は
[「DBは研究のインフラでありフロンティア」](https …
Continue reading »
Written by bonohu in misc on 土 27 5月 2017.
生命科学データベースの構築・管理・運営を担う
の4つの組織が連携して、ここ数年合同講習会を開催してきたのが、今回は三島での開催。自分は、それの統合TV撮影部隊として業務参加。
第一部は、前日まで行われていたINSDCメンバーのうち、NCBIの3人とEBIの1人が今日まで残ってくれて素晴らしい話をしてくれた。その感動したフレーズをtwitterで流した。
[データベースは学問のインフラ=公共基盤であり、普段無意識に使って居るが、使えなくなって初めて重要性に気がつく](https://twitter.com/bonohu/status/868321448359636992)
[全世界の科学研究のためにオープンデータ=科学情報の公開・共有を](https://twitter.com/bonohu/status/868322604322398208)
[あなたのデータの最もよい使い道は、あなた以外の他の誰かが思いつくかもしれない](https://twitter.com/bonohu/status/868332772896555008)
とくに追加説明はいらないだろう …
Continue reading »
Written by bonohu in misc on 金 26 5月 2017.
また別のtranscript alignerのSpliced Transcripts Alignment to a Reference (STAR)によるgenome mappingを試してみた。まずは、STAR用のReference genomeのindex作成から。
[shell]
time STAR --runThreadN 12 --runMode genomeGenerate
--genomeDir genome4star --genomeFastaFiles hogenome.fa
[/shell]
メモリが足りないといわれたら、--limitGenomeGenerateRAM=3200000000のようなオプションを追加して。パラメータの渡し方が他とは違って=で指定しないといけないらしい。
そして、mapping本体は以下のように。
[shell]
time STAR --runThreadN 24 --genomeDir genome4star
--readFilesCommand pigz -dc --outFileNamePrefix hogenome
--readFilesIn fuga_1.fq …
Continue reading »
Written by bonohu in misc on 木 25 5月 2017.
ゲノム配列も読んである場合に使えるGenome-guided Trinity De novo Transcriptome Assembly。まずは、TopHatやSTARなどでリードをゲノム配列マッピング(RNA用)して、その結果のBAMファイル(以下の例ではhogenome.bam)でde novo transcriptome assemblyを実行できる。
[shell]
time Trinity --genome_guided_bam hogenome.bam
--genome_guided_max_intron 10000
--max_memory 64G --CPU 24
[/shell]
Continue reading »
Written by bonohu in misc on 水 24 5月 2017.
SAMファイルをmergeしてからBAMファイルに変換してsortしようとしたが。一度、BAMに変換してからsortするしかないのか?
まず1ファイルづつ、SAM->BAM変換。
[shell]
for f in .sam;
do g="${f%.}"
time samtools view -@ 4 -bS $f > $g.bam
done
[/shell]
それらをmergeして、そのままsort。
[shell]
time samtools merge -@ 4 - *.bam
| time samtools sort -@ 4 -o merged.bam -
[/shell]
sortする際のtemporary領域を/tmp以下に指定した方がなおよいかも。具体的には-T /tmp/hogeのオプションををsortのコマンドに足す。
Continue reading »
Written by bonohu in misc on 火 23 5月 2017.
RNA-seqのためのtranscript assemblyとquantificationのプログラムStringTie。またオプションが変わっているかもしれんが、前に動かした時のそれ。
[shell]
stringtie fuga.bam -p 4 -o fuga.gtf -G hogenome.gff -A fuga_abd.txt
[/shell]
-Gで指定しているhogenome.gffはGuideとなるreference annotation。-Aのそれがgene abundance estimation。-oでassembled transcriptsの出力を指定(GTFファイル)。
Continue reading »
Written by bonohu in misc on 月 22 5月 2017.
hisat2でreference genomeにmappingする場合。hisat2のウェブサイトにすでにindexずみのそれがある場合はしなくていいが、まずはindex作成。hisat2-buildコマンドにて。
[shell]
hisat2-build -p 4 hogenome.fa hoge
[/shell]
そして、実際のmapping。
[shell]
hisat2 -p 4 -x hoge -1 fuga_1.fastq -2 fuga_2.fastq -S fuga.sam
[/shell]
出力はSAM形式であることに注意。
Continue reading »
Written by bonohu in misc on 金 19 5月 2017.
pip updateコマンドはないので、pipで入れたパッケージのアップデートは以下のようにする。
[shell]
pip3 list --outdated
| awk '{print $1}'
| xargs pip3 install -U
[/shell]
そろそろ、python3をデフォルトのpythonにして、version3系のpipをpip3と打たないで済むようにしたいところ。
Continue reading »
Written by bonohu in misc on 木 18 5月 2017.
近頃はChromatin ImmunoPrecipitation(ChIP)データの再利用がしやすくなっている。それらのデータをアノテーションしているDBとして、UCSC Genome BrowserのTrackにも入っているthe Open Regulatory Annotation database (ORegAnno)に注目していたが。現状temporalなページになっていて、だがデータはダウンロードできる状況。そのうち、復活するのだろうか?
Continue reading »
Written by bonohu in misc on 火 16 5月 2017.
統合TVのコンテンツ(.movファイル)だけでなく、Togo Picture Galleryのコンテンツも生命科学系データベースアーカイブに移行。統合TV関連コンテンツの完全クラウド化が実現した。そしてついに、この日を持ってtdiary版統合TVウェブインターフェースが廃止となった。「togotv+キーワード」でググった時にそちらが優先的に出るのがずっと気になっていたのだが、今回の廃止でついに現在の三代目インターフェースに統一されることに。関係者の皆さん、おつかれさまでした!
AOEの更新は結局終わらず、継続課題に。
Continue reading »
Written by bonohu in misc on 月 15 5月 2017.
AOEの追加機能の計算をpythonでやろうとするが、データがでかすぎて思うようにできず。できる手段を探しつつ、データをまとめて減らすことも検討。
メタデータの仕様に変更があったようで、それの原因究明。やはり、これまで取れていたメタデータがなくなっている。どこからそれを取ってくるか、要検討。
Continue reading »
Written by bonohu in misc on 金 12 5月 2017.
きっとあるだろうと思ったら、やっぱりググって出て来たこのページを参考に。
transposeというコマンドもあるらしいが、ここは手堅くawkで。とおもったら、メモリ不足で途中でkillされた。以下のdatamashコマンドを使うやり方だと数万x数十万の行列の転置が(少々かかったが)できた。
[shell]
brew install -v datamash
datamash transpose < matrix.txt > transposed.txt
[/shell]
いよいよpythonでのコーディング覚醒の悪寒。
Continue reading »
Written by bonohu in misc on 水 10 5月 2017.
みんなのPython勉強会に参加してきた。今回、奇しくも第24回目。月1回なので、ちょうど丸2年ということだった。
内容的にはそれぞれの方の会社で開発されているツールやそれにまつわる周辺情報。この種の会への参加が久しぶりだったので、勉強になったし、新鮮だった。
Continue reading »
Written by bonohu in misc on 月 08 5月 2017.
Mishima.syk開催打ち合わせで、三島広小路のあの店へ。その際の写真はこちらに。個人的には、生姜まみれのホタルイカが圧巻だった。次回はついに10回目のMishima.syk。いろんな意味で期待。
Continue reading »
Written by bonohu in misc on 日 07 5月 2017.
BLATはThe BLAST Like Alignment Toolで、UCSC Genome Browser のサイトにあるリファレンスゲノム配列に特化した配列類似性というか配列マッピングツールである。ゲノムランディングツールとも呼ばれる。だが、商用利用にはライセンスが必要なためか、便利なのだが広まっていない。そういうツールだからHomebrewにはまさか入っていないだろう、と。
ぶらっと
[shell]
brew install -v blat
[/shell]
してみたら、インストールが始まった…。商用でなければアカデミア、非商用、個人利用はライセンスいらない模様。ちなみに
[shell]
blat refgenome.fa query.fa output.psl
[/shell]
という感じで使い、出力はPSL形式で<kbd.output.pslに。BWAやbowtie、BLASTのように実行前に特別なindexingは必要ない。
Continue reading »
Written by bonohu in misc on 土 06 5月 2017.
配列パターン検索はリファレンスゲノム中に探すのであれば、ゲゲゲのゲノムことGGGenomeでやるのが手っ取り早い。ゲノム中で一気に探しておいて、その領域がどういった場所であったかは後で絞り込むというやり方で。
AGGTCANNNTGACCTというパターンに一塩基ミスマッチを許してヒトリファレンスゲノム(hg38)中に探すのであれば、以下のURLで。
<a href="http://gggenome.dbcls.jp/hg38/1/AGGTCANNNTGACCT" target="_blank">http://gggenome.dbcls.jp/hg38/1/AGGTCANNNTGACCT</a>
また、この結果を大量取得するには。例えばGFF3で保存するときにはURLの最後に.gffをつければGFF3形式で保存できる。
[shell]
curl -O http://gggenome.dbcls.jp/hg38/1/AGGTCANNNTGACCT.gff
[/shell]
てな具合に。結果はカレントディレクトリにAGGTCANNNTGACCT.gffというファイル名で。GFF3形式の他に、単なるテキスト(txt …
Continue reading »
Written by bonohu in misc on 金 05 5月 2017.
配列パターン検索はUNIX的にはgrepでできるが、バイオな配列に対してだと、パターン中に改行が入った場合やヘッダ行中の「誤爆」を防ぎたい。さらには、いくつかのミスマッチも許容するには、EMBOSSパッケージのfuzznuc(塩基配列)やfuzzpro(タンパク質配列)を使えばよい。例えば、1塩基のミスマッチまで許して、AGGTCAというパターンをFASTA形式のファイルhoge.faに探す際には以下のようにする。
[shell]
fuzznuc -sequence hoge.fa -pattern AGGTCA -pmismatch 1 -outfile hoge.fuzznuc
[/shell]
-pmismatchというオプションがキモ。
Continue reading »
Written by bonohu in misc on 木 04 5月 2017.
もちろん、ChIP-seqのデータがあればそれを利用すればいいのだが、多くの場合そういったデータのない非モデル生物種では、転写因子の結合の有無を調べるのに、転写因子結合サイトを予測する。
TRANSFACがそのデータベースとして老舗だが、有料になっている。
しばらくフォローしてなかったけど、JASPARが良くなっている。JASPAR CORE databaseとしてVertebrata, Nematoda, Insecta, Plantae, Fungiと生物グループごとにセットが分けられていて便利になっている。それ自体は狭山茶やっていた10年前に比べて増えているが、予測法自体は変わってない模様で、やはり閾値は自分で決めないといけないのと、生物学的にはfalse positiveが多い。
Continue reading »
Written by bonohu in misc on 水 03 5月 2017.
モデル生物ではきっちりアノテーションがなされ、UTR(UnTranslated Region)の配列抽出とか、EnsemblのBiomartを使えば簡単にできる。Ensemblにない生物種でも、Ensembl Genomesのそれを使えば良いのであるが、こちらの場合生物種によってはUTRのアノテーションがなくて抽出できないことがある(あった)。アノテーションがきっちりなされていない非モデル生物のUTRの配列抽出は大変である。
しかしながら、簡単にやる方法があった。それなりにdeepなRNA-seqデータがある場合に、であるが。それはTrinityによるde novo transcript assemblyとその結果を元にOpen Reading Frame(ORF)を予測するTransdecoderによるアノテーションを利用するというものである。Transdecoderを実行(過去のブログエントリ参照)した後に出て来る結果のGFF3形式の出力をBEDファイルとして保存して、それを元に部分配列抽出する。
[shell]
grep UTR Trinity.fasta.transdecoder.gff3 > UTR.gff3
bedtools getfasta -fi Trinity.fasta -bed UTR.gff3 -fo UTR.fasta …
Continue reading »
Written by bonohu in misc on 火 02 5月 2017.
割とすぐに終わるからこれまで特に気にしていなかったが。調べてみたら、やはりsamtools indexもスレッドオプションがあった。他のコマンドと同じで-@のあとに上限スレッド数を指定する。この例の場合、4。
[shell]
for f in *.bam;
do samtools index -@ 4 $f
done
[/shell]
並列化の効果あって、結果が得られるのが早くなった。「indexがない！」と別のアプリケーション(例えば、IGV)で怒られてindexを作ることが多いので、早く返ってくるのは嬉しいかと。最初から作っとけよ、という話もないではないが…。
複数のファイルを引数指定できるといいのだが、上述のように書けば済む話なのでよしとする。
Continue reading »
Written by bonohu in misc on 月 01 5月 2017.
今日から5月。絶好のデータ解析日和ということで(^_^)。
TopHatの結果を処理するコマンドは以前はやっつけでinteractiveに処理していたが、バッチ化というか並列化というか。中間ファイルがかさばるのでパイプライン処理して一気にsortされたファイルだけを出力しようということで先日覚えた-(マイナス)オプションの練習がてら。
状況としては、TopHatの結果のBAMファイルがhoge1,hoge2,hoge3のようなディレクトリの中にaccepted_hits.bamというファイルで入っているのが前提で、そのディレクトリがあるところと同じ階層にhoge.bamというファイル名で新規のソートされたBAMファイルを作成する。
| #!/bin/sh
type=$1
p=4
tmp=/tmp
samtools merge -@ $p - ${type}*/accepted_hits.bam
| samtools sort -@ $p -T $tmp/$type.$$ -o $type.bam -
|
というスクリプトをsamtools-merge_sort.shという名前で保存して、
[shell]
sh samtools-merge_sort.sh …
Continue reading »
Written by bonohu in misc on 金 28 4月 2017.
ファイル変換weekになってしまったので、BAMからCRAMへのバッチスクリプトも紹介しておく。
SAMBAM変換とは異なり、リファレンスゲノム配列が必要で、それは各環境で違う場所にあると思うので、それは自分の環境のそれを指定しないといけないことに注意。
| #!/bin/sh
p=4
gpath="/somewhere/reference_genome.fa"
for f in *.bam;
do g="${f%.*}"
echo $g
time samtools view -@ $p -T $gpath -C -o $g.cram $g.bam
done
|
上記のスクリプト(bam2cram.shとする)を、BAMファイルの置いてあるディレクトリに移動してから、実行する …
Continue reading »
Written by bonohu in misc on 木 27 4月 2017.
SAMからBAMに変換して、そのまま中間ファイルを作らずにBAMをソートする。それをバッチで処理するには。
SAMファイルの置いてあるディレクトリに移動(cd)してから、以下のようなスクリプトを実行。
| #!/bin/sh
p=4
tmp=/tmp
for f in *.sam;
do g="${f%.*}"
echo $g
samtools view -@ $p -bS $g.sam | samtools sort -@ $p -T $tmp/$g.$$ -o $g.bam -
done
|
samtoolsの行の一番最後の-(マイナス記号)がポイント。これが標準入力から入力を受けることを示すコマンドで …
Continue reading »
Written by bonohu in misc on 水 26 4月 2017.
現在お世話になっている研究所でOxford Nanopore Technologies のセミナーがあり、参加させていただいた。
大変興味深く聞かせてもらったが、MinIONからの塩基配列データはFAST5という形式で出てくるらしい。そのFAST5から、一般によく使われる配列フォーマットFASTQへの変換をするには、poretoolsというのを使えばよい、となかのひとに教えてもらった。そのporetoolsはやはりHomebrewにあって、
[shell]
brew install -v poretools
[/shell]
でインストールできる。
Usage examplesにあるように
[shell]
poretools fastq test_data/*.fast5
[/shell]
で変換可能。
Continue reading »
Written by bonohu in misc on 火 25 4月 2017.
ぼうのブログによると、2009年の7月29日に初めてDBCLSに来てくれた当時M1の学生さんだった。その後、Research Assistantとして統合牧場でUNIXとしてのMacの使い方を身につけ、そして学んだことをブログや統合TVとしてまとめてくれた。その後、そのままDBCLSに残ってくれて、後輩RAの指導、ときには自ら統合TVのコンテンツとなり、統合DBプロジェクトを盛り上げる一方、Sequence Read ArchiveのデータをFastQCで計算してその結果を可視化する、というプロジェクトに挑み、それを論文という形にすべく頑張ってくれた。
そして本日、ついに筆頭著者論文をpublish。おめでとう!今後、ますますの活躍を期待しています。
Continue reading »
Written by bonohu in misc on 月 24 4月 2017.
ついにこの日がやってきた。各所に働きかけた末に。
That's one small step for (a) man, one giant leap for mankind.
ただの盛り上がったミーティングに過ぎないと思われるが、実はそうでない。
これから快進撃が続くことは、すべて「ゼーレのシナリオ通り」なのである。
Continue reading »
Written by bonohu in misc on 金 21 4月 2017.
配列データ解析、頑張らないと。それを後押しする事象発生。折しもそれ関連の計算を始めていたのは、虫が知らせたか?自分自身の解析からそういうのが見つけられるかどうか、正念場。しまっていこう。
Continue reading »
Written by bonohu in misc on 木 20 4月 2017.
昨日のエントリの続きで、複数のBAMファイルをCRAMに変換し、そのサイズを比べてみた。
それぞれのファイルサイズを同一行に来るように1行ごとにデータを作って(cram.txtとbam.txt)、その圧縮率をちゃらっとawkで計算。
[shell]
paste cram.txt bam.txt | awk '{ print $1,$2, $1/$2 }'
[/shell]
結果が以下の通り。見ての通り、3カラム目がその圧縮率となる。
5660626395 9438014937 0.599769
4654129817 7896095631 0.589422
5087289649 8493424101 0.598968
5002310872 8382697420 0.596742
4117487398 7031097146 0.585611
4507556734 7565563208 0.595799
502823996 …
Continue reading »
Written by bonohu in misc on 火 18 4月 2017.
BAMからCRAMへの変換を試してみた。samtoolsを使えば良いだけの模様。
[shell]
samtools view -@ 4 -T hogenome.fa -C -o hoge.cram hoge.bam
[/shell]
3,091,833,154byteあったファイルサイズが、2,325,565,061byteに。約75%になったとは、すごい!CPU時間的には
263.48s user 45.58s system 205% cpu 2:30.72 total
複数CPUを指定効果もあったようだ。これぐらいの時間なら普段はCRAMにしておいて使うときだけBAMというのが実現可能か?
また、逆にCRAMからBAMへの変換もsamtoolsで。
[shell]
samtools view …
Continue reading »
Written by bonohu in misc on 火 18 4月 2017.
dottupの説明を書くため、EMBOSSのチュートリアルを見ていたら、XL23808というアクセッション番号の配列が使われていた。チュートリアルに出ている配列だし、自分もこれで例を作ろうと、これが何かをDDBJ/ENA/GenBankで検索しても、ググっても出てこない。
別ページにあったFASTA形式のファイルのヘッダにある機能アノテーション情報
Xenopus laevis rhodopsin gene, complete cds.
からググって解決。なんと、IDにもinsertionが入ってXLU23808になっていたというオチ。
そういえば、RefSeqも始まった当初から数字の部分の桁数が増えて、システムによっては異なる遺伝子がリンクされたりということがあったな。
諸行無常である。
Continue reading »