毎度忘れるので覚え書き $ bwa index -p dir/index_name reference.fasta

通常RNA-seqしたらreferenceのfastaファイルを使ってmappingして発現解析なりするわけだが、referenceが完備されていない種のseqはどうするのか？ 近縁種のreferenceを使う、というのが簡単な手段なわけだが、今回mappingしてみるとmap rateが30％くらいしか…

bedtoolsでbamファイルのread coverageを求める方法は以前にも書いた。 現在のbedtoolsのバージョン(2.27.1)ではbedファイルではなくbamファイルから直接read coverageを出すことが可能になっている$ bedtools coverage -a referense.bed -b sample.bam > sa…

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

以前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という感…

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

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

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

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

さて、データにアノテーションを付けていく作業は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…

いやあ、一つテンプレートとなるアプリを作っておけば、内容が違っていてもそれなりのデータベースアプリがすぐに作成できてしまうな。 前のやつには結局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…

$ 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を解析することになったのだが、やはりと言うか、色々と環境が整備されていなくて、いちいち対応していく必要が…