MEMEが見つけるモチーフの中身

たとえリピートをマスキングしてから実行したとしても、MEMEは沢山のモチーフを提示してくれる。

ランダムにサンプリングしたプロモーター配列をMEMEに渡すと、MEMEが発見するモチーフの傾向がつかめる。

転写開始点付近のGCに富む領域、Alu配列の下流にあるA-richな領域などが目立つ。A-richな領域はリピートマスキングでもマスクされないのでどうしてもMEMEが拾ってしまう。

ハズレらしきモチーフは極端に長い傾向にある。特にGCに富む領域の場合、ひとつの配列に類似のモチーフが複数見つかる。もちろん、ハズレと断定できるわけではないのだが、このような傾向があることを知っておくのは重要だと思う。

PowerPointの低い解像度について

macでKeynoteを使っていると、PowerPointが使い難くてたまらない。数ある問題のひとつ、「貼付けたPDFの画質が落ちる」という問題への対処法。せっかくPDFで貼付けているのに、どうしてガビガビにしてしまうのか。液晶プロジェクタで出力するなら高解像度は必要ないが、配布資料の画質が荒いのは困る。


使っているプレゼンソフト
Keynote'09
PowerPoint 2008


解決法
PowerPointの「環境設定」の「保存オプション」の「解像度の設定」で解像度を高める。「グラフィックファイルの圧縮」を解除する。


PDFをソースからコピーし、プレビューで「クリップボードから新規作成」し、PNGで保存する。このときのPNGの解像度が重要。解像度 300 ピクセル／インチ以上が望ましい。


.pngファイルをPowerPointに貼付ける。


Keynoteにはめ込んだPDFほどではないが、それなりにクリアになる。

for()ループを入れ子にして検索するより、Rならではの方法で

IDと数値をもつ表があり、それとは別にIDのリストがあるとする。そのとき、リスト中のIDに対応する数値を表から抽出するにはどうしたらいいだろうか。

Rを使う。

for()ループでIDリストを１つずつ取り出し、そのループの中で更にfor()ループを使って表を１行ずつ取り出し、対応するIDをif()で探すという方法が思い浮かぶ。for()ループを入れ子構造にするやり方だ。方法を下に示す。

この入れ子構造でやってみるとものすごく遅い。２つ目のfor()ループは表の行数マイナス１行分も無駄に動いているわけだから、遅いのも頷ける。でも、perlだったらもっと速くに結果がでるような気がする。どうやらRはループが苦手らしい。

作業効率がとても悪いので、for()ループをひとつ減らして同じ結果を得る方法を考えた。if()で条件検討をせず、セルが一致する行をそのままRっぽく拾う。方法を下に示す。

ここに示す小さなデータではスピードの違いを実感できないが、数千行の表の中から数百のリストに合致する行を得ようとすれば、その違いは雲泥の差。前者ならかけっぱなしで散歩にでも行きたくなるが、後者ならあくびをしている間に終了する。


とりあえず、説明用にデータを用意する。

データフレーム（df.1）を用意する。

id1 <- c("A","B","C","D","E","F")
value1 <- rnorm(6,mean=3,sd=1)
df1 <- data.frame(id1,value1)

> df1
  id1   value1
1   A 2.484163
2   B 1.539250
3   C 2.773143
4   D 1.041669
5   E 2.409677
6   F 2.532448


IDリスト（id2）を用意する。

id2 <- c("B","D","E")

> id2
[1] "B" "D" "E"


入れ子構造のやり方。
forループを２つ重ねて検索する。

found1 <- df1[0,]
for (i in 1:length(id2)){
for (j in 1:nrow(df1)){
if (id2[i]==df1[j,1]){
hit1 <- df1[j,]
found1 <- rbind(found1,hit1)
}
}
}

> found1
  id1   value1
2   B 1.539250
4   D 1.041669
5   E 2.409677



Rの強みを活かすやり方。
forループをひとつだけにする。

found2 <- df1[0,]
for(x in 1:length(id2)){
hit2 <- df1[df1[,1]==id2[x],]
found2 <- rbind(found2,hit2)
}


> found2
  id1   value1
2   B 1.539250
4   D 1.041669
5   E 2.409677

MEMEの前にリピートをマスク

MEMEを使ってプロモーター中のモチーフ探しをすると、意味がありそうな、なさそうな、保存性の高いエレメントが多数ヒットすることがある。これらは通常とても長い。ゲノムにはAluなどのリピート配列が多く、それらが解析対象のプロモーター配列セットの中に沢山含まれてしまっているためだ。

MEMEでDNAモチーフを探すときは、あらかじめリピートをマスクしたFASTAを用いるべきだろう。


リピートのマスキングには、giri（Genetic Information Research Institute）のRepeat Maskingツール「CENSOR」が使える。生物種（Sequence source）を指定し、Report simple repeatsをチェックし、FASTAを与えて実行すると、リピート部分の塩基を「X」に置き換えたFASTAを返してくれる。もちろん、どのようなリピートがどこにヒットするかも教えてくれる。

都合のいいことに、CENSORが返してくれたFASTA（「X」でリピートがマスクされたもの）は、そのままMEMEに与えることができる。

.wigのtag countをmedianでノーマライズする

MACSが出力したwiggleのtag countや、samtools mpileup -Dあるいはsamtools depthで得られたdepthを比較解析したいとき、ノーマライズのひとつのやり方として、medianで割るという方法が考えられる。

Rを使えば簡単。

wiggleを読み込む。ヘッダー２行は省く（skip=2）。

wig <- read.table("<in.wig>",skip=2)

> head(wig)

  V1  V2

1  1   3

2 11  34

3 21  49

4 31  65

5 41  81

6 51 105

２列目のデータを、２列目のmedianで割った値で置き換える。

wig$V2 <- wig$V2/median(wig$V2)

> head(wig)

  V1          V2

1  1 0.007575758

2 11 0.085858586

3 21 0.123737374

4 31 0.164141414

5 41 0.204545455

6 51 0.265151515

小数点以下２桁にする

wig$V2 <- round(wig$V2,digits=2)

> head(wig)

  V1   V2

1  1 0.01

2 11 0.09

3 21 0.12

4 31 0.16

5 41 0.20

6 51 0.27

まとめてやるなら、

wig <- read.table("<in.wig>",skip=2)

wig$V2 <- round(wig$V2/median(wig$V2),digits=2)

必要に応じてread.tableで書き出したり、プロットしたりする。

IGVの調整によりwiggle等を開けるようにする

IGVの調整によりwiggleを開けるようにする

MACSが作成したwiggle（.wig）のヘッダー情報にある染色体名（chrom=）がIGVで用いるアセンブリと一致しない場合、wiggleのヘッダーを書き換える以外に、IGVの設定を調整するという解決策もある。

macなら、ホームフォルダに「igv」というフォルダがあり（IGV 2.0のときは不可視化されていた？）、その中の「genomes」に、エイリアスファイルなるものを作って入れてやるだけで良い。

詳しい説明はこちら。

Creating a Chromosome Name Alias File

.wigと染色体名が一致しない問題は、BWA等によるアライメントの時に使用するリファレンスゲノム（.fasta）の染色体名をあらかじめIGVに合わせておくことで回避できるが、様々なソースから得たGFFやBEDなどもまた、染色体名の記述方法が異なる場合がある。その都度書類を書き換えるより、上記のごとくIGVの設定を変更しておくと楽。

IGVでwiggleを開く

MACSが解析したChIP-seqのデータ（.wig）をIGVで眺める。

必要なデータ

　<in.wig.gz>　MACSで作成したwiggle（.gzを解凍してもしなくても）

IGV（Integrative Genomics Viewer）はJavaで動くゲノムブラウザ。IGVのダウンロードページから、Binary distributionをダウンロードする。

IGV_2.1.17.zip（14.6 MB）を解凍し、igv.jarを起動。

生物（アセンブリ）を指定する。染色体名を確認する。「chr1」や「I」など、アセンブリにより染色体名の記述の仕方が異なる。

MACSが作成した.wig.gzを解凍して.wigにし、エディタで開いてヘッダーの２行目の「chrom=」の後ろを確認する。ここがIGVの記述方法と一致している必要がある。一致していなければ書き換える。

染色体名が一致していれば、.wig.gzのままでもIGVで開くことができる。

もうひとつのやり方はこちらを参照。