1. Mark duplicate란?

library prep과정에서 DNA조각이 여러번 복제되는 PCR 중복이 일어남. 이렇게 중복된 read는 특정 클러스터에 많이 나타나는 것으로 인식하게 되어 해당 부위의 유전자의 열림 정도를 과대해석하게 함. 따라서, 기술적으로 이를 보정하는 과정임.

2. Peak calling란?

ATAC-seq 데이터에서 열린 유전자 영역을 식별하는 과정. 이를 위해 매핑된 리드들의 밀도 분포를 분석하여 각 유전자의 열림 정도를 추정함.

3. Mark duplicate processing

① 자동화 script 확인하기

$ more bowtie_ATAC_mdup.pl

 more bowtie_ATAC_mdup.pl
@files = <*.raw.bam>;

foreach $file (@files)
{
$file=~s/.raw.bam//;

print "`samtools sort -@ 4 -o ".$file.".srt.bam ".$file.".raw.bam`;\n";
print "`samtools index ".$file.".srt.bam`;\n";
print "`bammarkduplicates markthreads=4 I=".$file.".srt.bam O=".$file.".rmdup.bam M=".$file.".rmdup.log index=1 rmdup=0`
;\n";
print "`rm ".$file.".srt.bam ".$file.".srt.bam.bai`;\n";
print "`samtools view -@ 4 -F 1804 -f 2 -q 30 -b ".$file.".rmdup.bam > ".$file.".final.bam`;\n";
print "`samtools index ".$file.".final.bam`;\n";
print "`rm ".$file.".rmdup.bam ".$file.".rmdup.bam.bai`;\n\n";

② 자동화 script를 이용해 모든 sample에 대한 명령어 생성

$ perl bowtie_ATAC_mdup.pl > bowtie_ATAC_mdup_auto.pl

$ more bowtie_ATAC_mdup_auto.pl
`samtools sort -@ 4 -o mESC_ATACseq_day0_1_atac2.srt.bam mESC_ATACseq_day0_1_atac2.raw.bam`;
`samtools index mESC_ATACseq_day0_1_atac2.srt.bam`;
`bammarkduplicates markthreads=4 I=mESC_ATACseq_day0_1_atac2.srt.bam O=mESC_ATACseq_day0_1_atac2.rmdup.bam M=mESC_ATACse
q_day0_1_atac2.rmdup.log index=1 rmdup=0`;
`rm mESC_ATACseq_day0_1_atac2.srt.bam mESC_ATACseq_day0_1_atac2.srt.bam.bai`;
`samtools view -@ 4 -F 1804 -f 2 -q 30 -b mESC_ATACseq_day0_1_atac2.rmdup.bam > mESC_ATACseq_day0_1_atac2.final.bam`;
`samtools index mESC_ATACseq_day0_1_atac2.final.bam`;
`rm mESC_ATACseq_day0_1_atac2.rmdup.bam mESC_ATACseq_day0_1_atac2.rmdup.bam.bai`;

`samtools sort -@ 4 -o mESC_ATACseq_day0_2_atac2.srt.bam mESC_ATACseq_day0_2_atac2.raw.bam`;
`samtools index mESC_ATACseq_day0_2_atac2.srt.bam`;
`bammarkduplicates markthreads=4 I=mESC_ATACseq_day0_2_atac2.srt.bam O=mESC_ATACseq_day0_2_atac2.rmdup.bam M=mESC_ATACse
q_day0_2_atac2.rmdup.log index=1 rmdup=0`;
`rm mESC_ATACseq_day0_2_atac2.srt.bam mESC_ATACseq_day0_2_atac2.srt.bam.bai`;
`samtools view -@ 4 -F 1804 -f 2 -q 30 -b mESC_ATACseq_day0_2_atac2.rmdup.bam > mESC_ATACseq_day0_2_atac2.final.bam`;
`samtools index mESC_ATACseq_day0_2_atac2.final.bam`;
`rm mESC_ATACseq_day0_2_atac2.rmdup.bam mESC_ATACseq_day0_2_atac2.rmdup.bam.bai`;

③ script 실행

$ ./bowtie_ATAC_mdup_auto.sh

$ more bowtie_ATAC_mdup_auto.sh
#!/bin/bash
module load SAMtools/1.3.1-foss-2016b
module load biobambam2/2.0.76-foss-2016b
perl bowtie_ATAC_mdup_auto.pl

4. Peak calling processing

① 자동화 script 확인하기

$ more macs_automate_ATAC_2018_narrow.pl

$ more macs_automate_ATAC_2018_narrow.pl
@files = <*.final.bam>;

foreach $file (@files) {
`macs2 callpeak -t $file -n $file -g mm -q 0.05 --nomodel --shift -100 --extsize 200 -B -f BAMPE`;


}

callpeak 옵션 (MACS2)
-t : 입력 BAM file
-n : 출력 파일명
-g : 맵핑된 유전체의 크기설정 (mouse는 mm, human은 hg)
-q : false discovery rate
--nomodel : 모델링 x
--shift : DNA 단편화 시점과 바인딩 사이트가 일치하지 않는 경우를 보정하기 위해 peak 위치조정
--extsize : 단일 염색체 read의 크기
-B : peak에 해당하는 BAM파일 출력
-f : BAM파일의 유형 (BAMPE : paired-end)

② script 실행

./macs_automate_ATAC_2018_narrow.sh

$ more macs_automate_ATAC_2018_narrow.sh
#!/bin/bash
module load MACS2/2.1.0.20150731-foss-2016b-Python-2.7.12
perl macs_automate_ATAC_2018_narrow.pl

③ 출력파일

filename_peaks.narrowPeak: peak 위치 정보를 가진 좁은 범위 peak 파일
filename_peaks.xls: peak 위치 및 통계 정보를 가진 엑셀 파일
filename_summits.bed: peak summit 위치 정보를 가진 BED 파일
filename_model.r: peak 모델을 나타내는 R 스크립트 파일
filename_control_lambda.bdg: control lambda 값을 포함하는 BDG 파일