kuroの覚え書き

96の個人的覚え書き

RNAseq

bowtie2-RSEM

bowtie2によるローカルアライメントを使い、Poly-Tプライマーで作成したcDNAライブラリのRNA-seq解析を行う。まずはbowtie2をインストールする。 $ git clone https://github.com/BenLangmead/bowtie2.git $ cd bowtie2/ $ make $ sudo make install次にリフ…

bwa index

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

RNA-seq のde novo assembly

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

bedtoolsでread coverage

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

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はバイナリなので中身が見れないが…

fastaから配列を取り出す

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

DDBJスパコンのbcftoolsが古い!

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

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

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

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…

データベースアプリ完成

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

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