science
リアライメントほぼ「教本」通りで完走。 モディファイした点は -Xmx4gのところを-Xms2g -Xmx2gと減らした samtools indexはやらなくてもbaiを作ってくれているような？RealignerTargetCreatorは-nt 12とスレッドの並列化オプションが使える模様IndelRealign…
現在のところ処理したいファイルが多数あり、それを同時進行で一気に処理を進める方向で並列化を設定している。 つまり、QC→トリミング→マッピングという流れを実験１〜ｎまで同時並行で実施するスクリプトを書いている。 実験をcsvファイルの表にしておき、…
このスクリプトだと多分数時間以上かかると思われるのでもっと並列化を目指す。 ちょうどHGCのサイトに参考になるスクリプトがあった https://supcom.hgc.jp/internal/mediawiki/アレイジョブでbwaの重い処理を分割手順としてはfastqファイルを分割 fastqは4…
「教本」に従ってサンプルファイルをダウンロードしてくる。諸般の都合でDDBJのサーバではなくhttp://www.ebi.ac.uk/ena/data/view/DRR006760 こっちからsraじゃなく、fastq.gzファイルを落とす。（若干手間省略）FastQC〜Trimmomaticの流れはこれまでのRNA-…
HGCのスパコンにてJava PATHの通っているバージョンはjava version "1.7.0_131" 1.8が必要なソフトのためにaliasを設定しておく alias java1.8="/usr/local/package/java/current8_32/bin/java"java version "1.8.0_141"csvkit (python) $ pip search csvkit…
メイン環境のiMacでcummeRbundが言うことを聞かなくなって暫く経つ。 現在、諸般の事情でサブ環境のMacbookAirで作業をしているのだが （主要な解析はスパコンやクラスタマシンを使用、グラフィカルな出力のあるRの部分はできるだけ手元のMacでやりたい） こ…
さて、ちっともまともに動いてくれないツール群の中で唯一走りきったっぽいCASHのデータを紐解いてみよう。こんな感じにexcelでも表示できる出力がtxtでできるので見てみる。左からID、Chr:position、exon numberと並んでいる。続く4カラムがちょっとよくわ…
今回は完全にDEXSeqだけのレシピでまずやってみる。 http://bioconductor.org/packages/release/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.pdf 環境はHGCのスパコン。pythonの環境は構築したのでRについてちょっと整備。 これまでRはローカルのMacでやっていた…
rMATS http://rnaseq-mats.sourceforge.net/user_guide.htm#output これに従ってやってみる。 サンプルデータはSRSF10のknockdownとcontrol STARのindexが必要とのことなのでダウンロード http://rmaps.cecsresearch.org/STAR/STARindex.tgz 70GB以上ある巨…
どれもこれもうまく動かないのでふてくされてても仕方がないので、一番最新のCASHを試してみる。 javaなようなので java -jar -Xmx8g ~/local/bin/cash_v2.1.0-alpha1/cash.jar --Case:prefix1 case.bam --Control:prefix2 Control.bam --GTF genes.gtf --Ou…
なにやらcummeRbundが使えなくなった。急に。 cuffdiffのデータ処理をミスったかと思ったが、以前に解析したデータでも同じエラーが出て止まってしまう。 > cuff <- readCufflinks("~/expression/SRSF10/hisat_results/cuffdiff_result") Creating database …
散漫的にやっていても埒が明かないので 実績のあるサンプルデータを使って、論文に書かれた方法を一からトレースしてみることにする。https://academic.oup.com/bib/article/doi/10.1093/bib/bbx034/3108818/CASH-a-constructing-comprehensive-splice-site…
featureCountの方もまだ不完全だがDEXSeqも手付かずなのでまずは環境構築 $ dexseq_prepare_annotation2.py -f DEXSeq.gtf expression/Homo_sapiens/NCBI/build37.2/Annotation/Archives/archive-2014-06-02-13-47-29/Genes/genes.gtf genes.DEXSeq.gtf Coul…
http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf を参考に。 $ ~/local/bin/subread-1.5.2-MaxOSX-x86_64/bin/featureCounts -T 4 -a ~/expression/Homo_sapiens/NCBI/build37.2/Annotation/Archives/archive-2014-06-02-13-47-29/Genes/g…
python php+SQL R これらを実用域で使えるようにする。pythonはまあshellよりはperlやRに近い文法っぽいな。 php+SQLは昔、MAMP(もしくはXAMPP)に手を出していたから大まかな概念は知っているつもりだったが、まあ既成のツールをいじっていただけだったので…
実際の解析を行うにあたりまずはアノテーションファイル（alternative splicingのカタログのようなファイル）を用意する。 いちおうMISOのサイトにもアップされているが、GRCh38バージョンのヒトゲノム情報とかはないので、自前で作製する方法を覚え書き。 r…
というわけで改めてやり直しまずはpython2.7.13をインストーラをダウンロードしてインストール（3.6もHomebrew使ってないので） 次に $ pip2 install misopy Collecting misopy Using cached misopy-0.5.3.tar.gz Collecting matplotlib (from misopy) Downl…
Alternative splicing, splice variant, isoformを検出、定量、比較する。 段取りとしてはFASTQのQC→トリミング→マッピング→AS analysis (& visualization)QC,トリミング,マッピングまではこれまでにやってきた組み合わせでいいはず。cufflinksで発現量比較…
$ cat A.bed chr1 0 20 chr2 0 20$ cat B.bed chr1 0 2 chr1 0 4 chr1 0 10 chr1 5 15 chr1 12 20 chr2 1 4 chr2 3 10 chr2 2 15というbedファイルが有ったとする。 $ coverageBed -d -a A.bed -b B.bed > cov.bed もしくは $ bedtools coverage -d -a A.bed…
Sun N1 Grid Engineで並列演算させる。 job投入用スクリプトではオプションで色々設定を与えることができる。 不用意に負荷の高いアレイジョブなどを投げるとリソースを使い切ってサーバーをダウンさせるので適切に制限をかけたほうがいい。 #! /bin/sh #お…
#!/bin/sh #csvファイルから実験名（alias, 2列目）を取り出し、コンマ区切りで並べて変数に格納 cut -d ',' -f 2 csvtest.csv | sed -e '1d' >tmp.txt list=$(paste -s -d ',' tmp.txt) echo $listcsvから変数を取り出す アレイジョブ編 #!/bin/sh #$ -cwd …
ということなのでgunzipとパイプで連携を考える mkfifo temp1.fastq mkfifo temp2.fastq gunzip -c hogehoge_1.fastq.gz > temp1.fastq & gunzip -c hogehoge_2.fastq.gz > temp2.fastq & bowtie -1 temp1.fastq -2 temp2.fastqこれでうまく走ることを確認し…
どっちもMacならhomebrewからインストールできる。 Trimmomaticは素でインストールするとjarなので $ java -jar /hogehoge/trimmomatic-0.36.jarとしなければならないがそれらもまとめてスクリプト化されているようで、普通に trimmomatic PE hogehogeで動い…
FastQCのコマンドを入れる際 #!/bin/sh list="exp1 exp2 exp3" mkdir -p ./FastQC for file in $list do fastqc -t 8 -q --nogroup -o ./FastQC ./"$file"_1.fastq.gz fastqc -t 8 -q --nogroup -o ./FastQC ./"$file"_2.fastq.gz doneこんな感じだけど "$fi…
exp1という実験のfastqファイルをQCすると exp1_1_fastqc.html exp1_1_fastqc.zip exp1_2_fastqc.html exp1_2_fastqc.zip の４つのファイルができる。 htmlファイルをブラウザで開けばそれぞれ個別にチェックできるが、大量のデータをいちいち開くのは勘弁し…
SRR000001~SRR000003の3つのファイルに処理を施したいとき、それぞれの処理が独立した処理なのでクラスタで並列に実施できる。 そのためにアレイジョブとしてスクリプトを記載。 #$ -t 1-3:1 でジョブ１〜３までプロセッサに空きがあったら順次実行する。 そ…
qsubでスクリプトを投げると hogehoge/local/bin/bowtie2-2.3.2/bowtie2-align-s: error while loading shared libraries: libtbb.so.2: cannot open shared object file: No such file or directory (ERR): Description of arguments failed!という感じでエ…
SRAファイルを展開する際 fastq-dump --split-files hoge.sra とする。
tophat -fusiontophat-fusion-post融合遺伝子の探索
サンプル間で有意な遺伝子発現の差（p > mySigMatrix <- sigMatrix(cuff, level = 'genes', alpha = 0.05) > mySigMatrix > library("psych") > setwd("/Users/hogehoge/expression/cuffdiff_results") > cuff_data <- read.table('/Users/hogehoge/expressi…
