kuroの覚え書き

96の個人的覚え書き

science

DNAのmultiple-FASTAファイルを翻訳してアミノ酸の配列を得る

個々のDNA配列からアミノ酸配列に変換するツールは数あれど、multi-FASTAにまとまっているDNA配列、それもATGーstopじゃなく5UTR-3UTRなどのcDNA配列から、読み枠をすべてサーチして最長のORFだけをリストにしてくれるツールはなかなかない。http://shigen.n…

bwaでtranscript.faにマッピングしたあとcufflinksで解析するには

bwaを使ってRNAseqのリードをtranscript.faをリファレンスとしてマッピング、出来上がったBAMファイルからcufflinksでFPKMを得ようとしたが、gtfファイルが無いので、FPKMを計算すると、マッピング領域のゆらぎを反映して1つのtranscriptの中にいくつかのloc…

塩基単位でのread depthを求める

以前BEDToolsのcoverageBedでdepthを各塩基ごとに求める方法を書いた。 $ coverageBed -a genome.bed -b sample.bed -dというのが基本書式だ。 genome.bedはreference.fasta.faiから $ awk '{print $1"\t0\t"$2}' reference.fasta.fai > genome.bedという感…

bamファイルから抽出

mappingを行ったbamファイルから更に絞り込みをかけてやる。 具体的にはSNV、INDELを持つreadを取り除いて100%matchしたreadだけにしたい。基本はsamtools viewでbamを開いて、grepやawkで整形をすることになるだろう。bamはバイナリなので中身が見れないが…

テスト運転

さて、早速 見せてもらおうか、富士通の中古サーバーの性能とやらを (注:ガ◯ダ◯って実は殆どみたことないんだが)大体270MBくらいのfastq.gzファイルを4つ trimmomatic SE \ -threads $cpu -phred33 -trimlog log.trimlog \ ./fastq_files/"$samp"_"$exp"…

CentOS7にNGS関連のソフトウェアをインストール

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

fastaから配列を取り出す

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

SQLとBLASTの連携

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

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