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