まずとにかくyumでwgetとnanoだけは入れておく。そのうえでlinuxbrewを入れてしまえば、ほぼ問題なく殆どの環境が構築できてしまった。 http://linuxbrew.sh Linuxbrewのインストール方法は上のサイトのままでOK幾つかのソフトは $ brew tap brewsci/science…

$ samtools faidx genome.fasta chr1:15,000-50,000 こんな感じでfastaファイルから取り出したい配列を指定して取り出すことが可能。 ある遺伝子を特定したとして、その上流プロモーターの配列をfastaから取り出したいなどの用途で使える。 また $ samtools …

まず、SQLで抽出された遺伝子のIDを元にtranscriptのFASTAファイルから配列を抽出。 抽出した配列をfasta形式で一時保存し、Arabidopsisの遺伝子データベースに対してtblastxでサーチ。 結果をwebに表示。こんな流れを構築した。 https://pypi.python.org/py…

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

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

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

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

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

早速、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でやるのが効率的なのだけど、これを視覚的にわかりやすく提示するにはちょっとごちゃごちゃしすぎているので、これをtext抽出してSQLに入れてやることにする。INFOやsampleカラムには複数の情報が詰まっ…

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

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

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

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

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

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

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

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…

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

$ 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:…

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

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

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テンプレートに任せたおかげで、ほぼ１日半で研究プロジェクトの公式HPの形が出来上がった。まあ、もうちょっと色合…

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

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…

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のアダプター配列…

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