アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
6月26日(金)全社研修による休業のお知らせ
平素よりアメリエフ株式会社のサービスをご利用頂きまして、誠にありがとうございます。
全社研修による休業のお知らせです。

【休業期間】
2015年6月26日(金)12:00〜18:00
【営業開始】
2015年6月29日(月)より通常営業


お問い合わせにつきましては、6月29日(月)より順次対応させていただきます。
ご不便をおかけ致しますが、何卒ご了承くださいますよう宜しくお願いいたします。
| ymm | 会社のこと | 16:04 | comments(0) | - |
Pythonとくんずほぐれつ
hatです。
平成27年度NGSハンズオン講習会で使用するPythonの資料を作っています。

人に教えるとなると改めてPythonを基本から学ぶ必要があり、土日もPythonの参考書やCodecademyで勉強している今日この頃です。

根を詰めすぎて逃避したくなったのか、あるいはPythonPythonと思いつめすぎてしまったのか、先週末、気が付いたら特急電車に乗って群馬のジャパンスネークセンターに向かっていました。

ジャパンスネークセンターには、大きいのから小さいの、黒いのから白いの、毒があるのから毒がないの、ガラガラ言うの、ありとあらゆる蛇がいました。Python(ニシキヘビ)もたくさんいました。

そこらじゅうで蛇同士が絡み合ってどれがだれやらわからない状態になっており、このようなスパゲッティプログラムは書きたくないものだと、胸を衝かれる想いがいたしました。

蛇を見まくり、毒蛇の採毒実演で息を飲み、大蛇に巻かれて記念写真を撮り、売店で買ったまむしドリンクを飲んで、非常にリフレッシュしました。
お疲れのみなさんにおすすめしたいスポットです。

蛇でチャージしたエネルギーを使って、引き続きPythonの資料作成を頑張りたいと思います。
| hat | よもやま話 | 16:11 | comments(0) | - |
paste関数を使う
短い間ですが、弊社のトレーニングのTAを担当させていただいていた時期があります。
RのTAをやるなかで、とても便利なのに、初めてRに触れる方が戸惑われることが多いと感じたポイントのひとつに paste 関数があります。
今回は paste 関数をあれこれいじくりまわしてみます。

1.paste関数とは

複数の文字列同士を結合する関数です。

> paste("hello", "world")
[1] "hello world"

なぜpaste関数が必要かというと、Rではこのように文字列を単にprint関数に渡すことができないためです。

> print("hello","world")
以下にエラー print.default("hello", "world") : 'digits' 引数が不正です
追加情報: 警告メッセージ:
In print.default("hello", "world") : 強制変換により NA が生成されました


対話モードだとprint関数を省略することも多いので、こんな書き方をされるかもしれません。いずれにせよ、複数の文字列を並べただけでは結合できません。
> "hello","world"
エラー: 予想外の ',' です in ""hello","

paste関数は大事です。

2.paste関数の基本

使い方の簡単な例をお見せします。
数字も、文字列に変換されて結合されます。

> paste(108, "yen")
[1] "108 yen"

結合した文字列の間に挟まる文字は、sep オプションで指定可能です。指定がない場合は半角スペースが入ります。

> paste(108, "yen", sep="-")
[1] "108-yen"

文字を挟みたくないときは、sep=""と指定します。

> paste(108, "yen", sep="")
[1] "108yen"


引数には、複数の要素があるベクトルを渡すこともできます。この場合、それぞれの要素について結合されます。

■ 1対複数
> pre
[1] "gene"
> str
[1] "a" "b" "c" "d" "e"
> paste(pre, str)
[1] "gene a" "gene b" "gene c" "gene d" "gene e"

■ 複数対複数(要素の数が同じ場合)
> num
[1] 1 2 3 4 5
> suf
[1] "st" "nd" "rd" "th" "th"
> paste(num, suf, sep="")
[1] "1st" "2nd" "3rd" "4th" "5th"

■ 複数対複数(要素の数が一致しない場合)…要素が足りないほうのベクトルが使いまわされます。
> paste(num, suf[1:4], sep="")
[1] "1st" "2nd" "3rd" "4th" "5st"

ベクトル同士ではなく、ひとつのベクトルの中の各要素を結合する場合は、collapseオプションを使用します。
collapseオプションに、結合する要素同士の間に挟まれる文字を指定します。

> day
[1] "2009" "July" "1"
> paste(day, collapse="-")
[1] "2009-July-1"

sepオプションとcollapseオプションを両方指定したらどうなるでしょう?
> paste(num, suf, sep="", collapse="/")
[1] "1st/2nd/3rd/4th/5th"

ベクトル同士を結合した後、各要素を結合していますね。

もう少し実践的な例として、for文とグラフでの使い方をみていきます。
続きを読む >>
| kubo | 統計解析ソフトR | 15:16 | comments(0) | - |
prefetch すらっと落とす SRA
こんにちは、久保(kubor)です。
先日転んだ時に頭を打ったので、それがきっかけで、僕に異常な言動がないかhatさんに確認してもらったのですが、
「打つ前から変わらずに変なことを言ってるよ。」
とのことでした。ひどい。

そんなhatさんが「すら(SRA)っとクイック(Q)に変換」を公開していたので、僕はこれに対抗して、
「NCBIのSRAから複数ファイルをまとめてダウンロードする方法」を紹介いたします。
組み合わせると”すらすら”とfastqを手に入れることが叶います。
今回は柑橘のシーケンスデータ「DRR008634 - DRR008641」8ファイルをダウンロードします。

まずは、NCBIのSRAからファイルリストを手に入れましょう。
どれでもいいのでアクセッション番号で検索して、「Run Selector」を確認します。


目的の8ファイルに辿りつけましたので、これらのアクセッション番号を手に入れます。
「Accession List」からテキストファイルをダウンロードします。



このファイルには、改行区切りでアクセッション番号が記述されています。
DRR008634
DRR008635
DRR008636
DRR008637
DRR008638
DRR008639
DRR008640
DRR008641

これを使用してSRA Toolkitでまとめてダウンロードしましょう。
使うコマンドは、「prefetch」です。SRA Toolkitに含まれていますので、ご確認ください。
prefetch --option-file SRR_Acc_List.txt

「--option-file」の引数で指定したファイルの中身のアクセッション番号を、順にSRAデータベースからダウンロードしてくれます。

ダウンロードされたファイルは、以下に保存されます。
$HOME/ncbi/public/sra

fastq形式への変換は「すら(SRA)っとクイック(Q)に変換」を御覧ください。
| kubor | バイオインフォマティクス | 17:17 | comments(0) | - |
heatmapの横に色をつけるオプション
heatmap関数のRowSideColors/ColSideColorsオプションの紹介です。
gplotsパッケージのheatmap.2やheatmap3パッケージでも同様の機能があります。
ヒートマップのデンドログラムの横に、カラフルなサイドバーを付ける機能です。

テスト用のマトリクスを作ります。
(結果は実行するたびにランダムに変わります)
mat<-matrix(sample(100,70,replace=T),ncol=7)
rownames(mat)<-c("geneA","geneB","geneC","geneD","geneE","geneF","geneG","geneH","geneI","geneJ")
colnames(mat)<-c("sample1","sample2","sample3","sample4","sample5","sample6","sample7")

このマトリクスでヒートマップを作成するとこうなります。
heatmap(mat)


以下は左にあるデンドログラムの横に、サンプルごとに色を付ける場合の例です。

RowSideColorsオプションに、色名(16進数での指定も可)のベクトルを指定します。
この色名のベクトルは、マトリクスの行数と同じ長さである必要があります。
nrow(mat)
[1] 10

また、色はマトリクスの各行に対して指定していますので、ヒートマップではクラスタリングされた順に並びかえられて描画されます。
rownames(mat)
[1] "geneA" "geneB" "geneC" "geneD" "geneE" "geneF" "geneG" "geneH" "geneI" "geneJ"   ←この順に色が指定されます
rcol<-c("red","blue","green","goldenrod","violetred","salmon","orange","springgreen","pink","black")
heatmap(mat,RowSideColors=rcol)


クラスタリングされたグループごとに色をつけたい場合は、マトリクスをあらかじめクラスタリングしたうえでcutree関数を使います。
mat.d<-dist(mat)
mat.h<-hclust(mat.d)
mat.c<-cutree(mat.h,k=3)
cutree関数のkで指定した数値と同じ長さの色名のベクトルをRowSideColorsに指定します。
clust.col<-c("red","blue","orange")
heatmap(mat,RowSideColors=clust.col[mat.c])


ColSideColorsの使い方も同様です。
なお、heatmap.2のlmatオプションを使用した場合のColSideColors / RowSideColorsの作図順については、以前記事を書きましたのでご参照ください。
| kubo | 統計解析ソフトR | 15:39 | comments(0) | - |
役員変更のお知らせ
この度、弊社は2015年6月8日をもちまして、下記の通り役員の変更がございましたのでお知らせいたします。
今後も、生命科学の課題解決に最適なITソリューションを提供する「プロフェッショナルチーム」としてお客様の研究のお役に立てますよう、役員並びに従業員一同、鋭意努力して参ります。
今後とも、より一層のご愛顧を賜りますようよろしくお願い申し上げます。

                 記

代表取締役社長CEO 山口 昌雄
取締役COO 金 景順
取締役CSO 須田 仁之
経営顧問 丸 幸弘 (株式会社 リバネス 代表取締役CEO)
アドバイザー 本多 喜久雄 (株式会社 グローバルコーチング 代表取締役)

                               以上
| ymm | 会社のこと | 14:44 | comments(0) | - |
ChemmineRを使ってみよう【4】
ChemmineRの紹介連載4回目です。
前回は、ChemmineOBというパッケージを使って化合物のPubChem fingerprintを取得しました。
今回では、そのfingerprintを使って、類似比較・クラスタリング解析を行います。

fingerprintによる類似検索

fpset[1]と類似したものをfpset全体から検索します。
fpSim(fpset[1], fpset, method="Tanimoto", cutoff=0.4, top=4)
650001
1
化合物のID類似度が出力されます。今回解析した10サンプル中に、cutoff以内で類似している化合物がなかったようです。

fingerprintを用いたクラスタリング解析
simMA <- sapply(cid(fpset), function(x) fpSim(fpset[x], fpset, sorted=FALSE))
hc <- hclust(as.dist(1-simMA), method="single")
plot(hc,hang=-1)


ヒートマップの作成
クラスタリングまでしたら、ヒートマップだって描きたいですよね。
library(gplots)
heatmap.2((1-simMA), Rowv=as.dendrogram(hc), Colv=as.dendrogram(hc), col=bluered(256), density.info="none", trace="none")



これで化合物間の類似関係が視覚的にわかりますね!
(10サンプルなので少々つまらない結果になりましたが……)


非常に簡単な紹介となりましたが、以上でChemmimeRに関する連載は終了です。
私にとって、あまり触れたことのないケモインフォマティクス解析の良い勉強となりましたが、皆様の研究の一助にもなれば幸いです。
| kubo | バイオインフォマティクス | 16:11 | comments(0) | - |
ChemmineRを使ってみよう【3】
前回までで、SDFデータの読み込みと、データのvalidationをご紹介しました。
今回は、読み込んだ化合物のfingerprint/atom pair descriptorを取得する方法です。

SDFsetからfingerprint/atom pair descriptorへの変換には、ChemmineOBという新しいパッケージをインストールします。
Open Babelというツールを扱うためのRインターフェースです。
biocLite("ChemmineOB")
library(ChemmineOB)
※ Open Babelのデータベースの関係で、OSによって挙動が違うようです。詳細はChemminOBのREADMEをご覧ください。

化合物のFingerprintの取得
fingerprintを取得します。
fpset<-fingerprintOB(sdfset,"FP2")
fpset
An instance of a 1024 bit "FPset" of type "FP2" with 10 molecules
fingerprint以外にもatom pair descriptorでも解析できます。
apset <- sdf2ap(sdfset)
apset

An instance of "APset" with 10 molecules

次回では、fingerprintを使った類似検索・比較を行います。
| kubo | バイオインフォマティクス | 16:25 | comments(0) | - |
My Favorite Things
みなさん、好きなPerlの演算子はなんですか?

私は、qw演算子が好きです。でも、「1..100と書くと、1から100までの整数を指すアレ」はもっと好きです。

これ(..)の名前を今日初めて知りましたが、「範囲演算子」と言うそうです。

これ(..)に名前があるなんて今まで考えたことがなかったのですが、人間が考えたものにはやはり何にでも名前がついているのだなあと思いました。

そして、(...)というのもあることを初めて知りました。世界は広いです。
| hat | よもやま話 | 15:19 | comments(0) | - |
 123456
78910111213
14151617181920
21222324252627
282930    
<< June 2015 >>

このページの先頭へ