Mapped reads (.bam) 파일을 갖고 있다고 가정한다. (36bp짜리 짧은 read들 mapping에는 bowtie1이 좋고, 긴 read들 mappin에는 bowtie2나 bwa-mem이 좋다고 한다. 특히 bwa-mem은 error rate이 높을 때도 꽤 잘 작동한다고 한다.)
보통 ChIP-seq 실험에는 IP하지 않은, 즉 antibody 처리 안 한 control (input) sample과, 관심 있는 protein에 대한 antibody 처리를 한 IP sample 두 개를 사용한다. → 현재 control.bam과 IP.bam 두 파일이 있다고 하자.
MACS2 개요
Removing redundancy
WXS에서와 같이 exact same location에 있는 read들을 하나로 친다.
The bad kind of duplicates
The good kind of duplicates
→ 결국 그냥 peak calling 전에 duplicate을 없애는 게 가장 좋다. Differential binding analysis 할 때는 duplicate을 retain해도 좋고, repetitive region에 binding하는 것이 예상된다면 duplicate이라 multimapped read를 retain하는 게 좋다.
Modeling the shift size
True binding site 주변의 read density를 보면 bimodal한 enrichment pattern이 나타난다.
MACS는 이러한 bimodal한 패턴을 보고 각 봉우리가 얼마나 shift되어야 '잘 합쳐진' 하나의 peak을 형성할지 모델링을 통하여 알아낸다.
가장 먼저, 어떤 지역을 분석할지를 정하는 게 필요하다. MACS는 이 과정에서 오직 ChIP 샘플만을 사용해서 (input/control 말고) significantly enriched region을 찾아낸다. bandwidth 길이의 슬라이딩 윈도우를 가지고 쭉 스캔하면서 genome에 random 하게 read가 mapping되었을 때의 density보다 mfold 배 이상 enrich되어 있는 지역을 골라내는 것이다.
그 다음, 이러한 high-quality peak 중 1000개를 random sample한다. 그 다음, positive (+) read와 negative (-) read를 분리하고, 각각을 따로 생각하여 두개의 봉우리 (mode)를 찾아낸다. 두개의 봉우리 사이의 거리를 d라고 정의하면 이것이 estimated fragment length가 된다. 각 read의 3' 방향으로 d/2 만큼 이동하면 그곳이 가장 그럴듯한 protein-DNA interaction site가 된다. → 1000개의 샘플을 가지고 d를 찾아내는 것이 핵심!
Scaling libraries
Input과 IP 샘플에서 sequencing depth가 차이가 나는 경우 MACS는 두 샘플의 total read count가 같아지도록 linear scale한다. 기본적으로 read count가 많은 샘플이 scaled down된다.
Peak detection
각 read를 d/2만큼 움직인 다음 MACS는 다시 2d 길이의 윈도우를 가지고 스캔하여 candidate peak를 찾아낸다. Genome 상의 read 분포는 푸아송분포로 모델링될 수 있다. 이때 lambda 는 그 window에서의 read 개수의 기댓값이다.
Genome 전체에서 균일한 lambda를 사용하는 것이 아니라, MACS는 lambda_local이라는 dynamic parameter를 사용한다. 이것은 각각의 candidate peak마다 정의된다. Lambda parameter는 control sample으로부터 계산되며, 다양한 window size를 시험해보고 가장 큰 값을 취한다. 즉, lambda_local = max(lambda_BG, lambda_1k, lambda_5k, lambda_10k)
결과적으로 lambda를 이용하면 local biase를 capture할 수 있고, 작은 region에 간헐적으로 발생하는 low tag count에 대해서 robust해진다. (low tag count는 local chromatin structure, DNA amplification, sequencing bias, genome CNV 등에 의해 발생한다.)
각 region은 p-value < 10e-5일 때 significant하게 read가 enrich되어 있다고 판단되며, p-value는 lambda를 parameter로 갖는 푸아송분포를 이용하여 계산된다.
끝으로 겹치는 peak들은 merge되며, 각 read position은 그 중앙으로부터 d만큼 extend된다.
가장 많은 fragment들이 mapping된 base 위치를 'summit'이라고 부른다. → Precise binding location이라고 볼 수 있다.
Fold enrichment는 ChIP-seq read count와 lambda_local의 비로 계산된다.
MACS2 사용법
MACS2 에는 다양한 subcommand가 있다.
•
callpeak
•
bdgpeakcall
•
bdgbroadcall
•
bdgcmp
•
bdgopt
•
cmbreps
•
bdgdiff
•
filterdup
•
predictd
•
pileup
•
randsample
•
refinepeak
그 중 메인은 callpeak이다. ChIP-seq 실험에 알맞게 다른 subcommand를 잘 조합해서 실행해준다.
macs2 callpeak -t H3K36me1_EE_rep1.bam -c Input_EE_rep1.bam \
-B --nomodel --extsize 147 --SPMR -g ce -n H3K36me1_EE_rep1
C++
복사
--nomodel and --extsize 147 tell MACS2 use 147bp as fragment size to pileup sequencing reads. -g ce lets MACS2 consider C elegans genome as background. Set it as 'dm' for fly or 'hs' for human. -B --SPMR ask MACS2 to generate pileup signal file of 'fragment pileup per million reads' in bedGraph format.
Callpeak options
•
Input file options
-t : IP 데이터 파일 (이것만 있으면 돌아감. The only required parameter)
-c : Control 데이터 파일
-f : Format of input file (default = AUTO라서 자동으로 찾아줌)
-g : Mappable genome size → Sequencing 가능한 genome 크기.
•
Output arguments
--outdir : 아웃풋 저장할 디렉토리
-n : 아웃풋 파일의 prefix
-B/--bdg : bedGraph 파일에 fragment pileup 정보와, control labmda, -log10pvalue, -log10qvalue를 저장한다.
•
Shifting model arguments
-s : Sequencing read의 size. 안 주어지면 처음 10개 read가지고 결정함.
--bw : Expected sonication fragment size. Model building할 때 윈도우 크기로 사용됨.
--mfold : High-quality peak 검출할 때 upper and lower limit으로 사용.
•
Peak calling arguments
-q : q-value (minimum FDR) cutoff
-p : p-value (instead of q-value cutoff)
--nolambda : lambda_local을 계산 안함
--broad : broad peak calling
Example commands
macs2 callpeak -t bowtie2/H1hesc_Nanog_Rep1_aln.bam \
-c bowtie2/H1hesc_Input_Rep1_aln.bam \
-f BAM -g 1.3e+8 \
-n Nanog-rep1 \
--outdir macs2
Shell
복사
이러면 시끄러우니까 stderr 파일로 저장하자
$ macs2 callpeak -t bowtie2/H1hesc_Nanog_Rep1_aln.bam \
-c bowtie2/H1hesc_Input_Rep1_aln.bam \
-f BAM -g 1.3e+8 \
-n Nanog-rep1 \
--outdir macs2 2> macs2/Nanog-rep1-macs2.log
Plain Text
복사
Output format
•
narrowPeak 포맷 (_peaks.narrowPeak)
BED 6+4 포맷이다.
•
나머지들
_peaks.xls : called peaks, pileup, fold enrichment 정보를 갖는 엑셀파일.
_summits.bed : 모든 peak의 summit 정보. Binding site motif 찾으려면 이 정보 활용하면 된다.
_model.R : model에 관한 PDF 이미지 뽑아내는 R 스크립트.
_control_lambda.bdg : bedGraph format for input sample.
_treat_pileup.bdg : bedGraph format for treatment sample.