히스톤은 진핵세포의 핵 내에서 DNA와 결합하여 크로마틴(chromatin)을 구성하는 구조 단백질이다(Wolffe and Guschin, 2000; Cutter and Hayes, 2015). 4종류의 히스톤 단백질(H2A, H2B, H3, H4)은 2개씩 결합하여 히스톤 팔량체(histone octamer)를 이루며, 이 팔량체는 약 146 bp의 DNA에 감겨 크로마틴의 구조적 기본 단위체인 뉴클레오좀(nucleosome)을 형성한다. 이 때, 히스톤과 DNA 사이의 정전기적인 인력이 뉴클레오좀을 형성하는 기본 원리이며, 이를 바탕으로 긴 DNA 가닥을 체계적으로 응축하여 크로마틴 구조를 만든다. 하지만 전사인자(transcription factor)나 RNA 중합효소(polymerase) 입장에서 보면 뉴클레오좀은 이 단백질들이 DNA에 결합하는 것을 막는 장애물이다. 이를 해결하기 위해 DNA와 히스톤 사이의 결합을 약화시키는 방법이 존재하는데, 히스톤 단백질의 화학적 변형이 중요한 부분을 담당한다.

히스톤 단백질의 N-말단 아미노산에는 아세틸기나 메틸기 같은 화학 그룹이 가역적으로 결합할 수 있는데, 이 화학적 변형들을 히스톤 변형(histone modifications)이라 일컫는다(Bannister and Kouzarides, 2011). 히스톤 변형은 아미노산의 위치와 결합하는 화학 그룹의 종류에 따라 다양하며, 직접적으로 크로마틴 구조를 변화시키거나 구조를 변화시키는 단백질을 모집하는데(recruit) 관여한다. 아세틸화(acetylation)와 메틸화(methylation)는 대표적인 히스톤 변형으로 주로 리신(lysine, K) 잔기에서 일어난다. 아세틸화의 경우, 히스톤 단백질의 양전하를 감소시키며, 그 결과 DNA와의 결합을 약화시켜 크로마틴 구조를 활성화시킨다. 이에 반해 메틸화는 변형이 일어나는 아미노산의 위치에 따라 그 역할이 다르다. 히스톤 H3의 4번 리신의 3-메틸화(H3K4me3)는 전사가 활발히 일어나는 프로모터(promoter)에서 전사인자의 모집을 도와주지만(Vermeulen et al., 2007; Lauberth et al., 2013), H3K27me3는 크로마틴 구조를 불활성화시키는 단백질을 불러들여 유전자 전사 억제에 기여한다(Cao et al., 2002; Min et al., 2003). 특히 히스톤 H3K27ac와 H3K27me3는 같은 아미노산에서 일어나지만 크로마틴 구조에 미치는 영향이 상반되기 때문에 이들 변형에 대한 많은 연구들이 수행되어 왔다(Kim and Kim, 2011; Bowman et al., 2014; Chrysanthou et al., 2022). 또한 히스톤 변형이 일어나는 유전체 상의 길이, 즉 분포 형태도 히스톤 변형의 종류에 따라 다르다(Landt et al., 2012). 예를 들어 H3K4me3와 H3K27ac는 유전자의 전사시작부위(transcription start site) 주변의 짧은(좁은) 지역에서 관찰되지만, H3K9me3나 H3K27me3는 상대적으로 긴(넓은) 지역에 연속적으로 분포하고 있다.

히스톤 변형은 암을 포함한 다양한 질병과 연관성을 보인다(Chi et al., 2010; Fernandes et al., 2020). 비정상적인 히스톤 변형은 유전자의 전사 상태를 변화시킬 수 있으며, 그 결과 질병을 초래할 수 있다. 제1형 당뇨 환자와 고혈압 유발 생쥐 모델에서 H3K27ac를 포함한 히스톤 H3의 아세틸화 수준이 비정상적으로 증가하였으며(Lee et al., 2012; Wang et al., 2019), 전립선암, 위암, 췌장암에서는 H3K27me3 증가에 의한 암 억제 유전자들(tumor suppressor genes)의 발현 감소가 보고되었다(Yu et al., 2007; Fujii and Ochiai, 2008; Ougolkov et al., 2008). 이러한 연구 결과를 바탕으로 히스톤 변형 효소 유전자를 표적으로 한 치료제가 개발되고 있으며, 실제로 히스톤 H3K27 변형을 담당하는 enhancer of zeste homolog 2 (EZH2)와 histone deacetylase (HDAC) 효소의 억제제가 항암제의 효과를 상승시킨다는 보고도 있다(Fiskus et al., 2009). 또한 유전체 수준의 히스톤 변형 연구를 통해 새로운 질병 관련 유전자를 동정해낼 수도 있다. 천식, 알츠하이머, 갑상선암 환자의 세포에서 수행한 유전체 전체의 H3K27ac 분포 연구는 새로운 질병 관련 유전자를 선별해냈으며, 치료제 개발의 표적 유전자 후보를 제시하였다(Marzi et al., 2018; McErlean et al., 2020; Zhang et al., 2021). 이처럼 히스톤 변형에 대한 연구 결과는 질병의 진단 및 예후 예측, 그리고 치료제 개발 등 다양한 분야에서 활용되고 있다.

히스톤 변형은 주로 면역침강법(chromatin immunoprecipitation, ChIP)을 이용하여 연구된다(Kuo and Allis, 1999). 이는 포름알데히드로 고정된 크로마틴을 뉴클레오좀 단위로 절단한 후, 원하는 히스톤 변형 항체와 반응시키고, 침전된 뉴클레오좀의 DNA를 정제하여 분석하는 기법이다. 과거에는 PCR 기법으로 ChIPed DNA의 특정 부분만을 분석하였으나, 최근에는 차세대염기서열분석(next-generation sequencing, NGS)이라고 일컫는 시퀀싱 기법이 개발되어 유전체 수준의 분석이 가능하게 되었다. 본 연구는 분포 형태가 다른 히스톤 변형들의 ChIP-seq 데이터를 적절히 분석하기 위하여, H3K27ac와 H3K27me3를 대상으로 NGS 데이터 peak 선별 프로그램들(MACS2와 SICER)을 비교하였으며, 이들 방법의 적절성을 peak로 선별되지 않은 부분과 비교하여 검증하였다. 또한 ChIP-seq 데이터를 다양한 용량(read 수)에서 분석하여 히스톤 변형 분포 형태에 따른 ChIP-seq 데이터의 질적 수준과 시퀀싱 용량(sequencing depth) 사이의 상관관계를 조사하였다.

Public K562 ChIP-seq 데이터 분석

사람 적혈구 유래 세포주 K562의 H3K27ac와 H3K27me3 ChIP-seq raw fastq 파일은 GEO database에서 확보하였으며(H3K27ac, GSE66023; H3K27me3, GSE167869), Galaxy 플랫폼을 이용하여 분석하였다(https://usegalaxy.org). 염기서열의 정확도가 낮은 reads는 Filter by quality와 FASTQ Quality Trimmer라는 분석 도구를 이용하여 제거하였고, 질적 수준(quality)이 보정된 reads를 Bowtie2를 이용하여 사람 hg19 canonical reference genome에 정렬(alignment, mapping)하였다(Langmead and Salzberg, 2012). Reads를 정렬한 bam 파일에서 mapping quality (MAPQ)가 20 미만이거나 PCR duplicate로 추정되는 reads를 제거하여 정확히 정렬된 reads만을 이후 분석에 사용하였다. ChIP-seq 데이터를 시각화하기 위해 BamCoverage를 이용하여 bigwig 파일을 생성하였으며(Ramírez et al., 2016), H3K27ac와 H3K27me3의 분포는 Integrative genomics viewer (IGV)에서 관찰하였다(Thorvaldsdóttir et al., 2013). UCSC database에서 제공하는 사람 유전자에 대한 bed 파일을 기준으로 전사시작부위 주변의 H3K27ac와 H3K27me3 분포 수준을 deeptools의 computeMatrix로 분석하고 plotHeatmap으로 나타내었다.

히스톤 H3K27 변형의 peak 선별

H3K27ac와 H3K27me3가 많이 관찰되는 지역(peak)을 선별하기 위해 총 세 가지 분석 방법을 사용하였다. 첫 번째로 MACS2에서 default parameter로 두 변형의 narrow peak를 선별하였고, 동일한 bam 파일에서 MACS2의 broad option을 사용하여 데이터를 분석하였다(Zhang et al., 2008). 마지막으로 SICER에서 default parameter를 이용하여 peak를 선별하고 비교하였으며(Zang et al., 2009), SICER로 분석하기 위해 bam 파일을 bed 파일로 변환하여 사용하였다. 모든 peak 선별 과정에서 ChIP 실험에 대한 input ChIP-seq 데이터를 control로 사용하였다.

Signal-to-background ratio 계산

Peak 선별 방법의 적절성을 판단하기 위해 signal-to-background ratio를 비교하였다. 세 가지 분석 방법에 의해 선별된 peaks를 기준으로 히스톤 변형이 일어나지 않는 지역(non-peak regions)을 ComplementBed로 지정하였다. Peaks와 non-peak regions에 정렬된 read 수를 측정하여 counts per million (CPM) 값을 계산하였으며(CPM=(특정 지역의 read 수)/((전체 read 수)/10⁶)), 이 때 전체 read 수는 flagstat을 이용하여 측정하였다. 그리고 Peaks의 평균 CPM 값을 non-peak regions의 평균 CPM 값으로 나눔으로써 signal-to-background ratio를 도출하였다.

세포 배양

본 실험에 사용한 K562 세포주는 ATCC에서 구입하였으며, 배지는 Gibco에서 구입한 시약들로 제작하였다. 사용한 배지는 10% FBS, 1% Penicillin/Streptomycin, 2 mM L-Glutamine, 20 mM HEPES가 포함된 RPMI 1640이었으며, 세포는 37℃, 5% CO₂ 조건에서 배양되었다.

크로마틴 면역침강법(ChIP)

ChIP 실험은 다음과 같이 수행하였다(Kuo and Allis, 1999). K562 (1×10⁷) 세포를 1% formaldehyde (Sigma) 조건에서 10분간 처리한(crosslinking) 후, 0.125 M Glycine (Invitrogen) 조건으로 반응을 중단하였다. Cell lysis buffer를 처리하고 10분간 얼음에서 반응한 후, 원심분리를 통해 핵을 분리하였으며, 분리한 핵에 MNase (Worthington Biochemical Corp)를 37℃, 15분간 처리하여 크로마틴을 단일 뉴클레오좀 수준으로 절단하였다. 핵을 분해한 후, 절단된 크로마틴을 Protein A agarose beads (Millipore)와 반응시켜(4℃, overnight) 면역침강 전에 beads에 비특이적으로 결합하는 크로마틴을 제거하였다(preclearing). 이후, 히스톤 H3K27ac 또는 H3K27me3에 특이적인 항체와 반응시키고(4℃, 3시간), Protein A agarose beads를 이용하여 항체에 결합한 크로마틴을 선별 분리하였다. 끝으로 phenol extraction 기법을 이용하여 DNA를 분리 정제하였다. ChIP 실험에 사용한 항체는 Abcam과 Millipore에서 구입하였다(H3K27ac, ab4729; H3K27me3, 07-449).

ChIP-seq 라이브러리 제작

ChIP-seq 라이브러리는 NEBNext Ultra II DNA Library Prep Kit (New England Biolabs)를 이용하여 제작하였으며, H3K27ac와 H3K27me3에 대한 ChIP-seq 라이브러리 제작은 본 연구실에서 기존에 수행했던 조건과 동일하게 진행하였다(Kang et al., 2021a; Kim et al., 2021). ChIPed DNA (H3K27ac, 10 ng; H3K27me3, 1 ng)의 양 말단을 비점착성 말단(blunt end)으로 만들고, 3' 말단에 아데닌을 추가한 후 AT base pairing으로 NEBNext Adapter를 연결하였다. Adapter가 연결된 DNA는 NEBNext Multiplex Oligos for Illumina (96 Unique Dual Index Primer Pairs)로 PCR을 수행하여 증폭시켰고 (H3K27ac, 7 cycles; H3K27me3, 11 cycles), NEBNext Sample Purification Beads로 약 175 bp 길이의 DNA를 선별한 후, 라이브러리를 제작하였다. 제작한 라이브러리의 농도는 Qubit dsDNA HS assay kit (Invitrogen)로 측정하였으며, NGS 시퀀싱은 NovaSeq 6000에서 100 bp paired-end 조건으로 수행하여 약 6천만개(60 M)의 reads를 얻었다.

ChIP-seq 데이터 분석

Paired-end 시퀀싱을 통해 얻은 H3K27ac와 H3K27me3 ChIP-seq 데이터 분석은 본 연구실에서 기존에 수행하던 분석 과정과 동일하게 진행하였으며(Kim et al., 2020; Kang et al., 2021b; Kang et al., 2021c), single-end 시퀀싱으로 얻은 public 데이터 분석 방법을 일부 수정하였다. Fastq 파일의 quality는 FastQC로 확인하였고, read 내 quality가 낮은 염기는 Trimmomatic으로 제거하였다(Bolger et al., 2014). Quality를 보정한 reads는 hg19 canonical reference genome에 Bowtie2로 정렬하였고, MAPQ가 20 미만이거나 PCR duplicate로 의심되는 reads는 제거하였다. ChIP-seq 데이터의 시각화를 위해 BamCoverage를 이용하여 bigwig 파일을 생성하였고, IGV를 통해 변형의 분포를 관찰하였다. ChIP-seq 데이터의 전체 read 수를 제한하기 위해 Downsample SAM/BAM을 이용하여 기존 bam 파일로부터 약 15 M, 30 M, 60 M, 120 M의 reads에 해당하는 bam 파일을 새로 생성하였으며, 생성된 bam 파일을 bigwig 파일로 변환하여 각 히스톤 변형의 분포를 관찰하였다. Read 수를 제한한 데이터에서 H3K27ac와 H3K27me3의 peak는 MACS2 broad와 SICER로 선별하였고, 선별한 peaks 사이에 중첩(overlap)되는 peak의 개수는 intervene으로 측정하여 venn diagram으로 나타내었다(Khan and Mathelier, 2017).

히스톤 H3K27 변형의 분포 형태 비교

히스톤 H3K27에서 일어나는 두 변형, H3K27ac와 H3K27me3는 크로마틴에서 그 분포 형태가 다른 것으로 보인다. 즉 H3K27ac는 주로 전사시작부위 같은 좁은 부분에서 관찰되지만, H3K27me3는 상대적으로 긴 부분에서 연속적으로 나타난다. 이러한 특징을 시각화하고 그 차이를 비교하기 위해 public K562 ChIP-seq 데이터를 분석하였다. Fig. 1A에서 보여지는 것처럼, H3K27ac는 짧은 길이에 걸쳐 분포하며, 내부에 밀집 부위가 별도로 존재하는 양상을 보였다. 하지만 H3K27me3는 상대적으로 긴 부위에 균일하게 분포하고 있었으며, 하나의 영역(domain, 도메인)을 형성하는 것처럼 보였다. 이러한 분포 형태가 유전체 전체에서 나타나는지 조사하기 위하여 유전자의 전사시작부위를 기준으로 H3K27ac와 H3K27me3의 분포 상태를 비교하였다(Fig. 1B). H3K27ac 수준을 기준으로 전사시작부위를 나열했을 때, H3K27me3는 H3K27ac 수준이 낮은 부위에서 관찰되었다. 두 변형의 배타적인 분포는 전사시작부위를 H3K27me3 수준으로 정렬했을 때도 관찰되었다. 그러나 이들 H3K27 변형의 분포 길이를 비교하면 H3K27ac는 전사시작부위의 ±2 Kb 이내에 주로 분포하지만, H3K27me3는 ±20 Kb까지 넓은 지역에 걸쳐 분포하고 있는 것이 관찰되었다. 이는 히스톤 H3K27 변형이 나타나는 길이, 즉 분포 형태가 H3K27ac와 H3K27me3 사이에서 크게 다름을 보여준다.

Fig. 1. Genome-wide distribution of histone H3K27ac and H3K27me3. (A) Genome browser shows distribution of histone H3K27ac and H3K27me3 in K562 cells. Public K562 ChIP-seq data were obtained from GEO database. Red and pale blue shadows represent H3K27ac- and H3K27me3-enriched regions, respectively. (B) Distribution of H3K27ac and H3K27me3 was visualized using density heatmap at transcription start sites (TSSs) ±5 Kb. TSSs were sorted by H3K27ac or H3K27me3 levels at TSSs ±5 Kb. Density heatmap was expanded to TSSs ±20 Kb with same ordering.

지금까지 다양한 종류의 히스톤 변형이 밝혀졌으며, 변형의 종류에 따라 유전체에서 분포하는 위치가 다른 것으로 나타났다. H3K27ac와 H3K27me3는 크로마틴 구조에 미치는 영향이 대비되기 때문에 이들 두 변형은 동일한 크로마틴 부위에서 함께 일어나지 않을 것으로 예상되며, 실제로 유전자와 인핸서에서 이들 히스톤 변형 분포의 역상관관계가 보고되었다(Wang et al., 2008; Katoh et al., 2018). 또한 히스톤 변형의 종류에 따라 분포 형태도 다른 것으로 보이며, 이는 Encyclopedia of DNA Elements (ENCODE) 프로젝트의 결과에서도 언급되었다(Landt et al., 2012). 따라서 H3K27ac와 H3K27me3의 경우처럼 그 분포 형태가 다르면 ChIP-seq 데이터의 생물정보학적 분석을 위해 히스톤 변형이 일어나는 지역, 즉 peak를 선별할 때, 분석 방법을 다르게 사용해야 할 것이며, 이는 유전체 수준에서 히스톤 변형을 제대로 해석하기 위해 중요한 과정으로 생각된다.

H3K27 변형 ChIP-seq 데이터에서 peak 선별 방법 비교

H3K27ac와 H3K27me3 ChIP-seq 데이터 분석 과정에서 peak를 적절하게 선별하기 위하여 총 세 가지 조건을 선정하고 비교하였다(Fig. 2A). 가장 일반적으로 사용되는 분석 도구인 MACS2와 넓게 분포하는 히스톤 변형에 주로 사용되는 SICER를 선택하였으며, MACS2에서는 narrow peaks (narrow) 옵션과 broad peaks (broad) 옵션을 비교하였다. 선별된 H3K27ac와 H3K27me3 peaks는 모두 각각의 히스톤 변형이 많이 분포하는 부위에 자리하고 있었다(Fig. 2B). 이 때 H3K27ac의 경우, narrow와 broad는 H3K27ac가 특히 많이 분포하는 밀집 부위를 우선적으로 peak로 선별하는 반면 SICER에서는 H3K27ac가 되어 있는 부위 전체가 넓게 선별되었다. 이는 고르게 분포하지 않는 H3K27ac의 분포 형태에 의해 나타나는 차이로 추정된다. 유전체 전체에서 선별된 peaks의 길이를 비교하면, narrow, broad, SICER 방법 순으로 그 길이가 증가하였다(Fig. 2C). Narrow와 broad로 선별한 H3K27ac peaks의 평균 길이는 491 bp와 1,915 bp였으며, 같은 방법으로 선별된 H3K27me3 peaks의 길이(433 bp와 1,686 bp)보다 조금 길었다. 하지만 SICER로 선별할 경우, H3K27me3 peaks의 평균 길이(6,036 bp)가 H3K27ac peaks의 평균 길이(5,170 bp)보다 길었으며, 이는 SICER에 의한 peak 선별이 H3K27me3의 넓은 분포 형태를 잘 반영하는 것으로 해석할 수 있다. 동일한 ChIP-seq 데이터로부터 길이가 다르게 선별된 peaks의 히스톤 변형 수준을 비교하기 위해 CPM 값을 계산하였으며, CPM 값은 peaks의 길이와 비례하여 증가하였다(Fig. 2D). CPM은 peak 내에 존재하는 read의 수를 전체 read 수로 보정한 값으로, peak 자체의 길이가 길수록 read 수가 증가하기 때문에 peak 길이와 CPM 모두 유사한 경향성을 보이는 것으로 생각된다.

Fig. 2. Peak calling using MACS2 and SICER tools in public histone H3K27ac and H3K27me3 ChIP-seq data. (A) H3K27ac and H3K27me3 peaks were identified using peak caller MACS2 and SICER with three analysis methods. (B) Peaks identified by each method were presented in genome browser of H3K27ac and H3K27me3 using different colors. Length of peaks (C) and counts per million (CPM) values of peaks (D) were graphed in boxplot. Box and whiskers represent interquartile range (IQR) and ±1.5×IQR values from edge of the box, respectively. (E) Signal-to-background ratio (S/B ratio) was calculated by analyzing H3K27ac or H3K27me3 level at each peak. (F) The S/B ratio of H3K27ac or H3K27me3 peaks was compared between three analysis methods. Difference of the S/B ratio was indicated in multiples.

Peak 선별 방법의 적절성을 파악하기 위해 signal-to-background ratio (S/B ratio)를 계산하였다(Fig. 2E). 히스톤 변형이 많이 분포하는 peak를 기준으로 하여 peak로 선별되지 않은 부위를 non-peak regions으로 구분한 후, peak와 non-peak region 내에 존재하는 read 수를 측정하였다. 측정한 read 수로 peak와 non-peak region의 CPM 값을 계산하였으며, peaks와 non-peak regions의 평균 CPM 값을 비교함으로써 S/B ratio를 도출하였다. 세 가지 분석 방법으로 선별한 peaks에서 S/B ratio를 비교한 결과, H3K27ac의 narrow와 broad, 그리고 narrow와 SICER 사이에서 그 비율은 각각 2배와 2.5배 증가하였지만, H3K27me3는 각각 7배와 18.8배 증가하였다(Fig. 2F). 이처럼 H3K27me3 ChIP-seq 데이터의 S/B ratio가 분석 방법에 따라 큰 차이를 보이는 것은 H3K27ac에 비해 H3K27me3 데이터가 peak 선별 방법에 의해 더 많이 영향 받을 수 있음을 시사한다. 따라서 넓게 분포하는 H3K27me3의 peak는 MACS2가 아닌 SICER를 이용하여 선별하는 것이 적절한 것으로 생각된다.

H3K27 변형 peak 선별 방법의 적절성 검증

Public K562 ChIP-seq 데이터를 이용하여 비교한 peak 선별 방법이 적절한지 확인하기 위해 본 연구실에서 수행한 K562 control (Con) H3K27ac ChIP-seq 데이터와 H3K27me3 ChIP-seq 데이터를 분석하였다. 2번의 반복 실험으로 얻은 Con ChIP-seq 데이터는 Spearman 상관관계 분석에서 나타난 것처럼 public ChIP-seq 데이터와 높은 유사성을 보였다(Fig. 3A, B). Con ChIP-seq 데이터에서 선별한 H3K27ac와 H3K27me3 peaks는 히스톤 변형이 많은 부위에서 나타났으며, peak의 길이도 narrow에서 SICER 순서로 증가하였다(Fig. 3C). Con ChIP-seq 데이터의 peaks에서 S/B ratio를 비교한 결과, H3K27ac는 모든 분석 방법에서 큰 차이가 나지 않은 반면(1.4~2.6배), H3K27me3는 narrow와 broad, 그리고 narrow와 SICER 사이에서 각각 5.4배와 9배 차이를 보여주었다(Fig. 3D).

Fig. 3. Peak calling in histone H3K27ac and H3K27me3 ChIP-seq data of control K562 cells. H3K27ac and H3K27me3 ChIP-seq were carried out in K562 control (Con) cells with two replicates (R1 and R2). (A) Scatter plots show correlation of public ChIP-seq data with Con R1 and R2 ChIP-seq data with Spearman correlation coefficients. (B) The correlation coefficients were also measured between H3K27ac and H3K27me3 ChIP-seq data and visualized by a heatmap. (C) Genome browser shows peaks identified by three analysis methods in green rectangles. (D) The S/B ratio of H3K27ac and H3K27me3 peaks was compared between three analysis methods with two replicates. Difference of the S/B ratio was indicated in multiples.

세 가지 방법으로 선별한 peaks의 S/B ratio는 H3K27ac에 비해 H3K27me3에서 큰 차이를 보였으며, 이러한 현상은 public 데이터와 Con 데이터 모두에서 나타났다. 이는 짧은 길이에 분포하는 H3K27ac에 비해 긴 거리에 넓게 분포하는 H3K27me3가 ChIP-seq 데이터 분석 방법에 의해 더 많이 영향 받음을 의미하므로 히스톤 변형의 분포 형태에 따라 적절한 peak 선별 방법을 사용하는 것이 필요한 것 같다. 아울러 public과 Con ChIP-seq 데이터 모두에서 H3K27ac의 S/B ratio가 H3K27me3보다 높았는데, 이는 ChIP 실험에 사용한 항체의 특이성과 관련된 것으로 추정된다. Fig. 1A와 Fig. 3C에서 보여지는 것처럼 H3K27me3 ChIP-seq 데이터의 noise signal이 H3K27ac 데이터에 비해 높게 나타나며, 높은 noise는 결과적으로 non-peak CPM 값을 증가시킨다.

Sequencing depth에 따른 H3K27 변형의 ChIP-seq 데이터 비교

ChIP-seq 라이브러리로 NGS 시퀀싱을 수행하면 일정한 길이의 reads를 얻게 되는데, read의 개수가 많을수록 sequencing depth는 증가한다. 하지만 시퀀싱 비용과 분석 시간이 증가하는 단점도 있다. 따라서 데이터를 효율적으로 분석하기 위해 적절한 길이와 개수의 read를 얻어야 한다. 분포 형태가 다른 H3K27ac와 H3K27me3의 ChIP-seq 데이터 분석에 적절한 sequencing depth를 알아보기 위해 Con ChIP-seq 데이터의 전체 read 수를 네 가지 기준(15 M, 30 M, 60 M, 120 M)으로 제한하고, 이후 분석을 수행하였다(Fig. 4A). Genome browser에서 보면 H3K27ac는 read 수와 상관없이 유사한 분포 형태를 나타내지만, H3K27me3는 60 M 이상의 read 수에서 길고 균일한 분포 형태를 보여준다(Fig. 4B). Sequencing depth와 peak 선별 사이의 관계를 조사하기 위해 H3K27ac의 peak는 broad 방법으로, 그리고 H3K27me3의 peak는 SICER 방법으로 선별하였다. H3K27ac는 60 M와 120 M reads에서 30 M에 비해 각각 1.9배와 3.6배 많은 peaks가 선별되었고, H3K27me3는 read 수와 상관없이 유사한 개수의 peaks가 선별되었다(Fig. 4C). 네 가지 기준의 read 수에 대하여 peak의 S/B ratio를 비교하였을 때 H3K27ac와 H3K27me3에서 15 M와 120 M 사이에 1.9배~2.6배의 차이를 보였으나, 이는 peak 선별 방법에 대한 S/B ratio의 격차보다 작다(Fig. 4D). H3K27ac에서 선별된 peaks의 위치를 비교하기 위해 중첩 여부를 분석한 결과, 60 M와 120 M에서 새로운 peaks가 많이 추가되었으며, 특히 120 M에서 700 bp 이하의 peaks의 증가가 두드러진다(Fig. 4E). 이는 많은 read 수에 의해 새로 선별된 peaks가 대부분 짧은 길이임을 알 수 있다. 반면 H3K27me3는 15 M에서 120 M까지 peaks의 대부분이 중첩되며, read 수가 증가할수록 1~4 Kb 사이의 peaks가 감소하고, 10~50 Kb 사이의 peaks가 증가한다(Fig. 4F). Fig. 4B에서 read 수에 따라 선별된 peaks를 H3K27ac와 H3K27me3의 분포 수준과 비교하면, H3K27ac 120 M에서 새로 선별된 peaks는 변형이 일어나는 부위를 제대로 반영하지 못하지만 H3K27me3가 넓게 분포하는 영역은 120 M에서 peak로 적절히 선별되었다. 따라서 히스톤 변형의 종류에 따라 높은 sequencing depth가 오히려 false-positive peak 선별을 초래할 수 있는 것으로 생각되며, H3K27me3처럼 넓은 부위에 분포하는 히스톤 변형은 sequencing depth의 증가가 적절한 peak 선별에 긍정적으로 작용하는 것으로 보인다.

Fig. 4. Analysis of histone H3K27ac and H3K27me3 ChIP-seq data in different sequencing depth. H3K27ac and H3K27me3 ChIP-seq data of control K562 cells (Con) were re-generated with 15, 30, 60 or 120 million (M) read counts. (A) The total read counts of Con ChIP-seq data were listed in a table. (B) Distribution of histone H3K27ac and H3K27me3 and peaks were visualized in genome browser with different read counts. Newly identified peaks in 120 M for H3K27ac were highlighted with red shadow. (C) The number of identified peaks was graphed in 15 M, 30 M, 60 M and 120 M. (D) The S/B ratio of peaks was compared with different read counts. Difference of the S/B ratio was indicated in multiples. Histone H3K27ac (E) and H3K27me3 peaks (F) identified in different read counts were overlapped and visualized by venn diagram (left). The number of peaks was graphed according to the length of them (right).

본 연구에서 분석한 Con ChIP-seq 데이터는 paired-end 시퀀싱을 통해 확보한 데이터로, 2개의 reads가 1개의 fragment에서 얻어지기 때문에 Fig. 4에서 언급한 15~ 120 M의 read 수는 7.5~60 M의 fragments를 반영한다. 그러나 single-end 시퀀싱으로 얻은 ChIP-seq 데이터는 1개의 read가 1개의 fragment를 반영하므로 NGS 데이터의 sequencing depth는 시퀀싱 방법에 따라 비교해야 한다. ENCODE 가이드라인에 따르면 mammalian cells에서 유의미한 ChIP-seq 데이터를 확보하기 위해 최소 1~2천만개(10~20 M) 이상의 uniquely mapped reads가 필요하다(Landt et al., 2012). 넓게 분포하는 히스톤 변형의 경우 sequencing depth가 높을수록 데이터의 질적 수준이 증가하는 경향을 보였는데, 이는 sequencing depth에 따라 다양한 히스톤 변형의 ChIP-seq 데이터를 비교한 연구 결과와 일치한다(Jeon et al., 2020). 따라서 좁게 분포하는 H3K27ac는 paired-end 기준 30 M의 reads만 얻어도 유의미한 ChIP-seq 데이터를 확보할 수 있지만, 넓게 분포하는 H3K27me3는 적어도 60 M 이상의 reads가 필요한 것으로 생각된다. 또한 read 수에 따른 S/B ratio 차이는 peak 선별 방법에 따른 차이보다 작게 나타났으며, H3K27ac 120 M에서 S/B ratio가 15 M보다 2.6배 높았지만 부적절한 peaks가 함께 선별되었던 것으로 보아 S/B ratio는 데이터 자체의 질적 수준 비교보다는 분석 방법의 비교에 활용하는 것이 적절하다고 사료된다.

본 연구 결과는 유전체 수준에서 히스톤 H3K27ac와 H3K27me3의 분포 길이가 다르며, 넓게 분포하는 H3K27me3는 좁게 일어나는 H3K27ac에 비해 peaks 선별 방법이 중요함을 보여준다. 또한 sequencing depth가 H3K27ac의 데이터 분석에 별 영향을 미치지 않지만, 너무 높은 sequencing depth는 peak 선별 과정에서 부정적인 영향을 미치는 것으로 나타났다. 반면 H3K27me3는 sequencing depth가 높을수록 ChIP-seq 데이터의 질적 수준이 증가하였다. 따라서 히스톤 변형에 대한 ChIP-seq 실험을 수행할 때 히스톤 변형이 분포하는 형태를 고려하여 sequencing depth를 결정하고, 분석 방법을 선택하는 것이 필요하다고 여겨진다. 이렇게 얻어진 결과는 정확한 생물학적 정보를 제공하고 나아가 질병의 진단 및 치료제 개발에 활용될 수 있을 것이다.

This study was supported by the 『BK21 Four Program』 of Pusan National University.

The authors have declared no conflict of interest.