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で開くことができる。

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

MACSでChIP-seqデータを解析

MACSはピークを検出するだけでなく、wiggle形式のデータも出力してくれる。

必要なデータ

　<treat.bam>　サンプル１

　<control.bam>　サンプル２

BWAとsamtools samseで、FASTQからBAMを得ておく。

他のパラメータ

　<genome.size>　ゲノムサイズを「1.5e+9」などの形で与える必要がある

　<name>　この名前のフォルダが作られる

コマンド

　macs14 -t <treat.bam> -c <control.bam> -f BAM -g <genome.size> -n <name> -w

作られるファイルとフォルダ

　<name>（フォルダ）

　　<name>_peaks.bed

　　<name>_peaks.xls

　　<name>_negative_peaks.bed

　　<name>_summits.bed

　　<name>_MACS_wiggle（フォルダ）

　　　control（フォルダ）

　　　treat（フォルダ）

＊ controlとtreatのフォルダに染色体毎のwiggle（.wig.gz）が作られる

＊ -wを与えなければwiggleは作られない

＊ コントロール（-c <control.bam>）は必ずしも必要ではない

.wigと.bedはIGVで開くことができる

SRAをFASTQに変換

NCBI SRA からダウンロードした.sraを.fastqにする。

SRA Toolkitを使う（たとえば、ver 2.1.9）

必要なデータ

　<in.sra>

コマンド

　fastq-dump.2.1.9 <in.sra>

作られるデータ

　<in.sra>と同じprefixの.fastq（in.fastq）

paired-endのfastqがセットになった.sraの場合、

fastq-dump.2.1.9 に 「--split-3」を与える

コマンド

　fastq-dump.2.1.9 --split-3 <in.sra>

作られるデータ

　in_1.fastqとin_2.fastq

SRAをFASTQに変換できない？
そう困ったときは、paired-endかもしれないと疑ってみるのがいいみたい。

BAMからdepthを得る

.bamからdepthを得る

必要なデータ

　ref.fasta　リファレンスゲノム

　in1.bam　サンプル＃１のデータ

　in2.bam　サンプル＃２のデータ
　.bamはsamtools view（必要ならrmdup）で得る

コマンド

　samtools mpileup -D -f <ref.fasta> <in1.bam> <in2.bam> > <mpileup.result1.txt>

生成されるデータ

　ref.fasta.fai

　mpileup.result1.txt

mpileup.result1.txtの内容

　１行目：染色体名

　２行目：塩基番号（1~）

　３行目：塩基（A, T, G or C）

　４行目：depth <in1.bam>

　５行目：リードの方向（"." or ","） <in1.bam>

　６行目：対応するクオリティー <in1.bam>

　７行目：depth <in2.bam>＊

　８行目：リードの方向（"." or ","） <in2.bam>＊

　９行目：対応するクオリティー <in2.bam>＊

　＊インプットが１つの場合、６行目まで

生成されるデータを軽量化する（リードの方向やクオリティーを省く）

コマンド

　samtools mpileup -D -f <ref.fasta> <in1.bam> <in2.bam> | awk 　'{print($1,$2,$3,$4,$7)}' > <mpileup.result2.txt>

mpileup.result2.txtの内容

　１行目：染色体名

　２行目：塩基番号（1~）

　３行目：塩基（A, T, G or C）

　４行目：depth <in1.bam>

　５行目：depth <in2.bam>

高速シーケンスデータからBAMを得てバリアントをコールする

BWAによるアライメントの続き

samtools & bcftools

アライメントデータ（.sam）をバイナリデータ（.bam）に変換し、バリアント（SNP, indel）をコールする

bcftoolsはsamtoolsについてくる

必要なファイル

　ref.fasta

　fastq1.txt

　samse.result.sam

　sampe.result.sam

samtools view

　samtools view -uS <input.sam> | samtools sort - <name>

　→　.bamファイルが作られる（name.bam：sort済みのuncompressedな.bam）

samtools rmdup

　samtools rmdup -s <imput.bam> <output.bam>

　→　PCR duplicatesが取り除かれる

samtools index

　samtools index <input.bam>

　→　.baiファイルが作られる

samtools flagstat

　samtools flagstat <input.bam> > <output.txt>

　→　リードのマップ状況などに関する統計が表示される

samtools idxstats

　samtools idxstats <input.bam> > <output.txt>

　→　各染色体にマップされたリードの数が表示される

samtools pileup

　廃止された（現在はmpileupを使用する）

samtools mpileup（variantのコール）

　samtools mpileup -uf <ref.fasta> <input.bam> | bcftools view -bvcg - > <output.raw.bcf>

　→　.bcfファイルが作られる（uncompressed BCF: binary call format)

bcftools view

　bcftools view <input.raw.bcf> | vcfutils.pl varFilter -D10000000 > <output.var.flt.vcf>

　→　フィルター済みの.vcfファイルが作られる

高速シーケンスのデータをリファレンスにアライメント

BWA（Burrows-Wheeler Aligner）を使って

生データ（.fastq）をリファレンス（.fasta）にアライメントする（.sam）

必要なファイル

　ref.fasta リファレンスゲノム（multiple fasta）

　fastq1.txt FASTQデータ

　fastq2.txt FASTQデータ（paired-endの場合はFASTQを２つ）

bwa index

　bwa index -a is <ref.fasta>

　→　８種類のファイルが作られる

　　　（.amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa）

bwa aln

　bwa aln <ref.fasta> <fastq1.txt> > <aln.result1.sai> | 

　bwa aln <ref.fasta> <fastq2.txt> > <aln.result2.sai>

　→　single-endの場合は上の１行のみ

　　　.saiファイルが作られる

bwa sampe（paired-endの場合）

　bwa sampe <ref.fasta> <aln.result1.sai> <aln.result2.sai> <fastq1.txt> <fastq2.txt> > <sampe.result.sam>

　→　.samファイルが作られる

bwa samse（single-endの場合）

　bwa samse <ref.fasta> <aln.result1.sai> <fastq1.txt> > <samse.result.sam>

　→　.samファイルが作られる

MEMEでモチーフ探し

ある一群のDNA配列が特定のエレメントをもつかどうかを知りたい。例えば、転写因子が結合するエレメントがあると期待している。そんなときは、MEMEが便利。


MEME (Multiple Em for Motif Elicitation)


MEMEに一群の配列を渡すと、その中に存在するエレメントをLOGOとして返してくれる。


 MEMEは、ダウンロードして自分のmacで走らせることができる。サーバの混雑に影響されずに済むし、CUIだとパラメータを少しずつ変更しながらの作業が楽。 


自分のmacで走らせると、指定したディレクトリにHTML書類（meme.html）を出力してくれる。同じディレクトリに発見されたモチーフの順鎖と逆鎖のPNG画像（logo1.png, logo_rc1.png等）があり、HTMLをブラウザで開くとそれらが参照される。


meme.txtと名付けられたファイルに重要な情報がまとめられている。こちらからコピペすれば、HTML越しに作業をするより楽。 meme.txtにはモチーフの定義が記述されており、これを切り抜いて個別のファイルにしておくとFIMO (Find Individual Motif Occurences)に渡すことができる。




アウトプットについての情報

MEME Suite Motif File Formats

モチーフのフォーマットはこちら
MEME Minimal Motif Format (sample DNA motif)

ファイルの行数を調査する

「テキトーな表」を使う。
<dir_name>と名付けたディレクトリに「テキトーな表」に由来するファイルが入っているとする。


ディレクトリ（フォルダ）に移動する。
cd <dir_name>


カレントディレクトリ内の、ファイル名が「.txt」で終わるファイルの行数を調べる
wc -l *.txt


$ wc -l *.txt
       8 dataA.txt
       8 dataA_col23.txt
       5 dataA_2over0.txt
       7 dataA_2plus3.txt
       8 dataB.txt
      36 total


「.csv」で終わるファイルなら、
wc -l *.csv


$ wc -l *.csv
       8 dataA.csv
       8 dataB.csv
      16 total


「dataA」で始まるファイルなら、
wc -l dataA*


$ wc -l dataA*
       8 dataA.csv
       8 dataA.txt
       5 dataA_2over0.txt
       7 dataA_2plus3.txt
       8 dataA_col23.txt
      36 total


特定のファイルの行数を調べる（例えば「dataA.txt」）
wc -l dataA.txt


$ wc -l dataA.txt
       8 dataA.txt

awkで特定の列や行を抽出する

「テキトーな表」を使う


awkでタブ区切りファイルを開き、２列目と３列目だけを抽出して新規ファイルとして保存する。
awk 'BEGIN{FS="\t"; OFS="\t"} {print $2,$3}' dataA.txt > dataA_col23.txt


上を実行する際に、１行目（ヘッダー行）だけを省く
awk 'BEGIN{FS="\t"; OFS="\t"} {if(NR >= 2) {print $2,$3}}' dataA.txt > dataA_col23_noHead.txt


あるいは
awk 'BEGIN{FS="\t"; OFS="\t"} {if(NR != 1) {print $2,$3}}' dataA.txt > dataA_col23_noHead.txt


あるいは
awk 'BEGIN{FS="\t"; OFS="\t"} {if(NR > 1) {print $2,$3}}' dataA.txt > dataA_col23_noHead.txt




２列目がゼロより大きい行を抽出して保存する。
awk 'BEGIN{FS="\t"; OFS="\t"} $2 > 0' dataA.txt > dataA_2over0.txt


２列目と３列目を足した値を４列目として保存する。
awk 'BEGIN{FS="\t"; OFS="\t"} $4=$2+$3 {print $0}' dataA.txt > dataA_2plus3.txt

awkによる、タブ区切りからcsvへの変換、およびその逆

「テキトーな表」を使う

awkでタブ区切りファイルを開き、CSVファイルとして保存する
awk 'BEGIN{FS="\t"; OFS=","} {print $0}' dataA.txt > dataB.csv

あるいは
awk 'BEGIN{FS="\t"; OFS=","} $0' dataA.txt > dataB.csv

awkでCSVファイルを開き、タブ区切りファイルとして保存する
awk 'BEGIN{FS=","; OFS="\t"} {print $0}' dataA.csv > dataB.txt
awk 'BEGIN{FS=","; OFS="\t"} $0' dataA.csv > dataB.txt

Rでタブ区切り（.txt）やcsv（.csv）を開く

「テキトーな表」を使う


タブ区切りファイルを読み込む
data.B <- read.table(file="dataA.txt",header=T,sep="\t")


あるいは
data.B <- read.table(file="dataA.txt",header=T)
data.B <- read.delim(file="dataA.txt",header=T)


read.delim() に header=T を与えなくても、ヘッダーがカラム名になる
data.B <- read.delim(file="dataA.txt")


> data.B
   name1  name2  name3
1 -1.220 -0.321 -0.585
2  1.584 -1.410 -0.378
3 -0.305  1.061 -0.392
4 -1.344  0.814  0.974
5 -0.798  0.139  0.742
6  0.688 -1.462 -0.149
7  0.023  0.342  2.005


read.table() に header=T を与えなければ、１行目にヘッダーが入ってしまう
data.B <- read.table(file="dataA.txt")


> data.B
      V1     V2     V3
1  name1  name2  name3
2  -1.22 -0.321 -0.585
3  1.584  -1.41 -0.378
4 -0.305  1.061 -0.392
5 -1.344  0.814  0.974
6 -0.798  0.139  0.742
7  0.688 -1.462 -0.149
8  0.023  0.342  2.005

Rで表をつくり、カラム名をつけ、タブ区切り（.txt）やcsv（.csv）として保存する

Rでテキトーな表を作る
data.A <- matrix(round(rnorm(21),digits=3),ncol=3)


rnorm: 乱数発生
round: 丸め処理（digits: 小数点以下の桁数）
matrix: リストから行列をつくる（ncol: カラムの数）


> data.A
       [,1]   [,2]   [,3]
[1,] -1.220 -0.321 -0.585
[2,]  1.584 -1.410 -0.378
[3,] -0.305  1.061 -0.392
[4,] -1.344  0.814  0.974
[5,] -0.798  0.139  0.742
[6,]  0.688 -1.462 -0.149
[7,]  0.023  0.342  2.005


カラム名をつける
colnames(data.A) <- c("name1","name2","name3")


> data.A
      name1  name2  name3
[1,] -1.220 -0.321 -0.585
[2,]  1.584 -1.410 -0.378
[3,] -0.305  1.061 -0.392
[4,] -1.344  0.814  0.974
[5,] -0.798  0.139  0.742
[6,]  0.688 -1.462 -0.149
[7,]  0.023  0.342  2.005


タブ区切りファイルとして保存するなら、
write.table(data.A,file="dataA.txt",sep="\t",col.names=T,row.names=F,quote=F)


csvとして保存するなら、
write.table(data.A,file="dataA.csv",sep=",",col.names=T,row.names=F,quote=F)


dataA.txtをエディタで開くと、


name1 name2 name3
-1.22 -0.321 -0.585
1.584 -1.41 -0.378
-0.305 1.061 -0.392
-1.344 0.814 0.974
-0.798 0.139 0.742
0.688 -1.462 -0.149
0.023 0.342 2.005


dataA.csvをエディタで開くと、


name1,name2,name3
-1.22,-0.321,-0.585
1.584,-1.41,-0.378
-0.305,1.061,-0.392
-1.344,0.814,0.974
-0.798,0.139,0.742
0.688,-1.462,-0.149
0.023,0.342,2.005