kuroの覚え書き

96の個人的覚え書き

science

BLASTでトップヒットだけを抽出する

アノテーションデータベースが作りたいので、A種のgene aについてB種の遺伝子をBLASTで調べて、一番上に来るものを1:1対応させたい。http://bioinfo-dojo.net/2016/03/25/blast_besthit_outfmt7/ ここにあるように、blastの結果から1行目だけをawkで取り出す…

BLASTを自前で

バイオインフォマティックスをやってるならBLASTくらい自分で構築しないとな、 というわけでもないけど、やはりアノテーションデータベースを構築したいので、blastをMacに入れることにする。 さらっとぐぐってみたらHomebrewで行けそうなので $ brew instal…

BLASTとの連携

生物種間で発現を比較しようとすると、アノテーションデータベースで橋渡ししてデータを見るということになるわけだが、アノテーションにはどうしても限界がある。 例えば生物Aで発現変動が見られた遺伝子について、生物Bでの発現を見る分にはアノテーション…

DDBJスパコンのbcftoolsが古い!

パイプラインが走らないとウンウンと苦悩したのだが、なんとbcftoolsのバージョンが0.1.17だった。ローカルにインストールしているのが1.7なので、古すぎる。どうりでbcftools callがunknown commandと言われるわけだ。2時間返せ〜

データの抽出作業でuniqをつかう

uniqは便利なコマンドなんだけど、OSXにプリインストールされているuniqはBSD版で機能がちょっと物足りない。他にもBSD版コマンドはちょくちょくGNU/Linux版と違うことがあるので、この際ということでGNU版をインストールしてしまうことにした。例によってHo…

DDBJのスパコンアカウントが使えるようになった。

早速、ssh→qloginでログインノードまで進む。 何が使えるかな? Python3 OKだな。 $ pip3 list biopython (1.66) cutadapt (1.9.1) cycler (0.10.0) drmaa (0.7.6) Genshi (0.7) matplotlib (1.5.1) nose (1.3.7) numpy (1.10.4) Pillow (3.1.1) pip (8.1.0)…

vcfからSQLに流し込む

さて、データにアノテーションを付けていく作業はvcfでやるのが効率的なのだけど、これを視覚的にわかりやすく提示するにはちょっとごちゃごちゃしすぎているので、これをtext抽出してSQLに入れてやることにする。INFOやsampleカラムには複数の情報が詰まっ…

vcfにアノテーションを付ける2

昨日の作業でsnpにアノテーションを付けたわけだが、snpデータを複数のcontrolサンプルからのデータにしたほうがいいかもしれない、と思い、vcfをmergeしてアノテーションをつけなおすことにした。複数サンプルのvcfをマージする方法はvcf-mergeを使うことに…

vcfにアノテーションを付ける

さて、無事samtools mpileup | bcftools callにてsnv callingが実施できたところで、これを解析するためのアノテーションをつけようと思う。 今回は自前で用意したcontrolのサンプルで検出されたsnp/indelをバックグラウンドのSNVとして対象のvcfにマーキン…

RNAseqからSNV解析

さて、ちょっと変則的にRNAseqデータからsnpを検出すべく、あれこれやっているわけだが、 普通にゲノムデータでSNV解析するならBWAでマッピングしてsamtools mpileupと行くのだろうが、ここをHisat2でマッピングしてmpileupと持っていこうと思う。 で、samba…

sambambaでmpileupを並列化

snp callを行うのにBAMファイルをsamtools mpileupで処理するのは以前にやった。(4年前) しかしこのmpileupは4コアあっても有効に使ってくれないのですごく時間がかかった記憶がある。これを高速に処理する方法はないものかと検索したところ、sambambaなる…

データベースアプリ完成

いやあ、一つテンプレートとなるアプリを作っておけば、内容が違っていてもそれなりのデータベースアプリがすぐに作成できてしまうな。 前のやつには結局5ヶ月ほどかけたが、今回のはほぼ2週間で出来上がった。細部少々荒いし、マルチユーザーを意識せず、ロ…

cummeRbundによるクラスター解析

以前に一通り試したクラスタリングを実際にやってみた。 忘れていたことも含めて覚書 > library("cummeRbund") #ライブラリ読み込み > cuff <- readCufflinks("cuffdiff_resultsのディレクトリ") #cuffdiffの出力ディレクトリからデータを読み込み) > setwd…

StringTieを入れる

TopHat→cufflinks→cuffdiff→cummeRbundというのがRNAseq解析の王道であるが、最近はHISAT2→StringTie→Ballgownというのが流行りはじめているらしい。これはやっとかないと。ということでまずは環境構築から。 StringTieのインストールはbrewでできるらしいの…

cuffmergeがエラー

Error (GFaSeqGet): end coordinate (515) cannot be larger than sequence length 504 Error (GFaSeqGet): end coordinate (515) cannot be larger than sequence length 504 Error (GFaSeqGet): subsequence cannot be larger than 506 Error getting subs…

cufflinksのFPKMを遺伝子ごとに集計して、サンプル間で比較する。

常法的にはcuffdiffでFPKMを比較し、cummeRbandでグラフ化するのだけど、サンプル数が多いと処理もけっこう大変になってくる。とりあえずFPKMを遺伝子ごとに集計して並べて見れるようにしてみることにする。cufflinksで出力される値に幾つか項目を加えてテー…

Nbのリファレンス

$ head Niben.genome.v1.0.1.contigs.fasta >Niben101Scf00001Ctg001 AAAAAAAGGATTAAGTGTCATAAATGTGGTAAATTTGGTCATTATGCAAGTGAGTGTAAA ACTCAGGAAAATATTAAGAGTCTAGATTTAGATGATAAACTTAAGGATTCTTTGTGTAAG ATTCTACTAAATTCTGATTATAGTTCTGATGTATCTGATTCCTCTTCTACTG…

リファレンスとアノテーションファイルの中身を理解する。

とりあえずちゃんとできているArabidopsisで使っているファイルを見てみると マッピングに使うリファレンス:Arabidopsis_thaliana.TAIR10.dna.toplevel.fa (IGVのreferenceにもそのまま使う) $ head Arabidopsis_thaliana.TAIR10.dna.toplevel.fa >1 dna:…

Nicotiana benthamiana

イネやアラビはゲノムデータが充実していて定法通りでNGSデータをマッピングできるのだけど、それ以外の植物だとどうだろう。 タバコ属植物Nicotiana benthamianaの場合、ゲノムがきちんと解析完了しているわけではなく、いまだcontig情報の塊でしかない。 t…

arabidopsisでhisat2

さて、これまでヒトゲノム研究でRNAseqパイプラインをせっせと作っていたのだが、縁あって植物科学の世界に舞い戻ってきた。そこで植物のRNAseqを解析することになったのだが、やはりと言うか、色々と環境が整備されていなくて、いちいち対応していく必要が…

BAMファイルを分割、結合

BAMファイルの一部分を切り取ってサイズの小さいBAMファイルを作ってみる。 (デモンストレーション用にコンパクトにしておく) $ samtools view -h -b 0123.bam 4:1-10000000 > 0123_4.bamこれで0123.bamのchr4の1から10000000の領域を切り出してあらたなba…

遺伝病

http://medi.atsuhiro-me.net/entry/2015/09/09/204159

ホームページ制作つづき

昨日着手したホームページだが、とりあえず入力完了。結局テキストエディタでひたすらタグを打って完成。 cssとjavascriptの大半をBootstrapテンプレートに任せたおかげで、ほぼ1日半で研究プロジェクトの公式HPの形が出来上がった。まあ、もうちょっと色合…

vcfファイルをいろいろなところからダウンロード

病気サンプルの原因遺伝子を特定するためのふるいとして、既知のSNPを除く、もしくは調べるという作業はまず必要となってくる。 そのために既知の病気の原因になっていると思われるSNPや、逆に普通の人が何%かの割合で持っていて、特に病気の原因とは考えら…

StringTieでread coverageを計算

Cufflinksに近いプログラムStringTieをHisat2でマッピングしたRNA-seqデータに使う。 Cufflinksにくらべて It uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts repr…

Trimmomaticのilluminaclipについて

Suggested adapter sequences are provided for TruSeq2 (as used in GAII machines) and TruSeq3 (as used by HiSeq and MiSeq machines), for both single-end and paired-end mode. とあるので、HiseqでシークエンスしたデータはTruSeq3のアダプター配列…

exome analysisやってみた4

一通り最後までできるようになったので、時間を測りながら一気に流してみたい。 解析に用いるファイルはfastq.gzの状態で4.5Gと4.6Gのサイズである。 qsubのオプションで -pe def_slot 16として16スロット確保。 FastQCはそれぞれ8スロットでペアを同時に…

exome analysis やってみる3

リアライメントほぼ「教本」通りで完走。 モディファイした点は -Xmx4gのところを-Xms2g -Xmx2gと減らした samtools indexはやらなくてもbaiを作ってくれているような?RealignerTargetCreatorは-nt 12とスレッドの並列化オプションが使える模様IndelRealign…

並列化

現在のところ処理したいファイルが多数あり、それを同時進行で一気に処理を進める方向で並列化を設定している。 つまり、QC→トリミング→マッピングという流れを実験1〜nまで同時並行で実施するスクリプトを書いている。 実験をcsvファイルの表にしておき、…

exome analysis やってみる2

このスクリプトだと多分数時間以上かかると思われるのでもっと並列化を目指す。 ちょうどHGCのサイトに参考になるスクリプトがあった https://supcom.hgc.jp/internal/mediawiki/アレイジョブでbwaの重い処理を分割手順としてはfastqファイルを分割 fastqは4…