Postępowanie doświadczalne
Patient xenografting and cell culture
Ludzkie linie komórkowe T-ALL obejmują CUTLL1 (dar od Adolfo Ferrando, Columbia University) i JURKAT (American Type Culture Collection (ATCC), Manassas, VA, #CCL-119). Komórki hodowano w podłożu RPMI1640 z l-glutaminą i 25 mM HEPES (Corning) uzupełnionym 10% inaktywowaną termicznie płodową surowicą bydlęcą (Sigma-Aldrich), 10 U/mL penicyliny-streptomycyny (Gibco) i 1× glutaMAX (Gibco) w nawilżanym inkubatorze w temperaturze 37 °C i 5% CO2. Komórki są okresowo testowane na obecność mykoplazmy przy użyciu zestawu Lonza Walkersville MycoAlert Mycoplasma Detection Kit (ostatni test w styczniu 2020). Linie komórkowe są utrzymywane w hodowli przez maksymalnie 20 pasaży i są uwierzytelniane przy użyciu profilowania krótkich powtórzeń tandemowych (JURKAT) lub przy użyciu PCR w celu wykrycia translokacji TCRb-NOTCH1 (TCRBJ2S4CUTLL1F:5′-GGACCCGGCTCTCAGTGCT-3′, NOTCH1CUTTL1R:5′-TCCCGCCCTCCAAAATAAGG-3′). Ostatnie uwierzytelnienie komórek przeprowadzono w lutym 2020 roku. Ludzkie komórki CD4+ T zostały zakupione od firmy AllCells. Pierwotne próbki ludzkie zostały pobrane za świadomą zgodą i analizowane pod nadzorem Institutional Review Board of Padova University, Associazone Italiana di Ematologia e Oncologia Pediatrica oraz pediatrycznych badań klinicznych Berlin-Frankfurt-Münster (AIEOP-BFM) ALL 2000/2006. Świadoma zgoda na wykorzystanie resztek materiału do celów badawczych została uzyskana od wszystkich pacjentów przy wejściu do badania zgodnie z Deklaracją Helsińską.
Przeciwciała i odczynniki
Western blots zostały wykonane przy użyciu następujących przeciwciał: Actin i CTCF firmy Millipore Sigma (klon C4; 07-729) oraz cleaved NOTCH1 (Val1744) firmy Cell Signaling Technology (4147). ChIP-seq przeprowadzono przy użyciu następujących przeciwciał: CTCF z Millipore Sigma (07-729), H3K27Ac (8173S), i H3K27me3 (9733S) z Cell Signaling Technology, i H3K4me1 (07-473) z Millipore.
In situ Hi-C
In situ Hi-C przeprowadzono na komórkach CD4+ T, Jurkat, CUTLL1, i ksenograftach pacjentów, jak wcześniej opisano . W skrócie, komórki zostały usieciowane 1% formaldehydem przez 10 min w temperaturze pokojowej. Dla każdej reakcji Hi-C, 5 milionów komórek było lizowanych, a jądra były permeabilizowane. DNA trawiono za pomocą MboI z New England Biolabs (R0147M). Strawione fragmenty znakowano biotynylowanym d-ATP firmy Jena Bioscience (NU-835-BIO14-S) i ligowano. Po potraktowaniu RNazą i Proteinazą K w celu odwrócenia wiązań krzyżowych, jądra poddano sonikacji przy użyciu aparatu Covaris E220 w celu uzyskania fragmentów o średniej długości 400 bp. Do ściągania fragmentów znakowanych biotyną użyto kulek streptawidyny firmy Thermo Fisher Scientific (65001). Po oczyszczeniu i izolacji DNA, ostateczne biblioteki przygotowano przy użyciu zestawu NEBNext® Ultra™ II DNA Library Prep Kit for Illumina® i sekwencjonowano za pomocą sekwencjonowania sparowanego końca przy długości odczytu 150 bp na Illumina HiSeq 2500, aby uzyskać średnio 400 milionów odczytów na próbkę.
ChIP-seq profiling
Komórki TCD4+, Jurkat, CUTLL1 i ksenograftów pacjentów sieciowano 1% formaldehydem i 1% płodową surowicą bydlęcą w PBS przez 10 min w temperaturze pokojowej. Reakcję wygaszano 0,2 M glicyną w temperaturze pokojowej przez 5 min. Komórki następnie przemywano PBS i granulowano.
Dla CTCF ChIPs, immunoprecypitację przeprowadzono w oparciu o protokół opisany wcześniej . Osad zawierający 50 milionów komórek lizowano za pomocą 5 mL buforu lizującego (50 mM HEPES-KOH, pH 7,5, 140 mM NaCl, 1 mM EDTA, 10% glicerol, 0,5% NP-40, 0,25% Triton X-100) przez 10 min w temperaturze 4 °C. Jądra granulowano przy 1350×g przez 7 min i ponownie zawieszano w 10 mM Tris pH 8, 1 mM EDTA i 0,1% SDS. Chromatynę ścinano przy użyciu systemu Covaris E220 do średniej długości fragmentu 400 bp i wirowano przy 15 000 obr/min przez 10 min w celu usunięcia nierozpuszczalnej chromatyny i resztek. Supernatant inkubowano z 20 μL Dynabeads Protein G przez 30 min przed odrzuceniem kulek. Jeden procent całkowitej objętości zapisywano jako wkład, a resztę inkubowano przez noc z przeciwciałem anty-CTCF. W sumie dodawano 100 μL Dynabeads Protein G przez 2 h. Związane fragmenty przemywano dwukrotnie 1 mL buforu niskosolnego (20 mM Tris-HCl pH 8.0, 150 mM NaCl, 2 mM EDTA, 1% w/v Triton X-100, i 0.1% w/v SDS), raz z buforem o wysokiej zawartości soli (20 mM Tris-HCl pH 8.0, 500 mM NaCl, 2 mM EDTA, 1% w/v Triton X-100, i 0.1% w/v SDS), raz z buforem chlorku litu (10 mM Tris-HCl pH 8.0, 250 mM LiCl, 1 mM EDTA, 1% w/v NP-40, and 1% w/v deoxycholic acid), and twice with TE (10 mM Tris pH 8, 1 mM EDTA).
For histone ChIPs, cells were lysed in 375 μL of nuclei incubation buffer (15 mM Tris pH 7.5, 60 mM KCl, 150 mM NaCl, 15 mM MgCl2, 1 mM CaCl2, 250 mM sacharozy, 0,3% NP-40, 1 mM NaV, 1 mM NaF, i 1 tabletka inhibitora proteazy bez EDTA (Roche)/10 mL w H2O) przez 10 min na lodzie. Jądra przemywano raz buforem trawiącym (10 mM NaCl, 10 mM Tris pH 7,5, 3 mM MgCl2, 1 mM CaCl2, 1 mM NaV, 1 mM NaF i 1 tabletka inhibitora proteazy bez EDTA (Roche)/10 mL w H2O) i ponownie zawieszano w 57-μL Digest Buffer zawierającym 4,5 jednostki MNase (USB) na 1 h w 37 °C. Aktywność MNazy była wygaszana przez 10 min na lodzie po dodaniu EDTA do końcowego stężenia 20 mM. Jądra były granulowane i ponownie zawieszane w 300-μL Nuclei Lysis Buffer (50 mM Tris-HCl pH 8,0, 10 mM EDTA pH 8,0, 1% SDS, 1 mM NaV, 1 mM NaF i 1 tabletka inhibitora proteazy bez EDTA (Roche)/10 mL w H2O) przed sonikacją przy użyciu Bioruptor Pico (Diagenode) przez 5 min (30 s włączone, 30 s wyłączone). Lizat odwirowywano przy maksymalnej prędkości przez 5 min w celu usunięcia resztek. Do supernatantu dodawano dziewięć objętości buforu do rozcieńczania IP (0,01% SDS, 1,1% Triton X-100, 1,2 mM EDTA pH 8,0, 16,7 mM Tris-HCl pH 8,0, 167 mM NaCl, 1 mM NaV, 1 mM NaF i 1 tabletka inhibitora proteazy bez EDTA (Roche)/10 mL w H2O). W sumie dodano 50 μl Dynabeads Protein G i próbkę inkubowano w temperaturze 4 °C przez 30 min, obracając. Jeden procent próbki pozostawiono jako wkład, a pozostałą próbkę rozdzielono do 3 probówek. W sumie do każdej probówki dodano 50 μl Dynabeads Protein G sprzężonych z 15 μl odpowiedniego przeciwciała przed inkubacją przez noc w temperaturze 4 °C, przy obracaniu. Kompleksy związane z kulkami przemywano po 5 min w 1 mL buforu o niskiej zawartości soli, buforu o wysokiej zawartości soli, buforu LiCl i dwukrotnie TE.
W celu elucji kompleksów związanych z kulkami, kulki zawieszano ponownie w 50 μL buforu elucyjnego (100 mM NaHCO3, 1% w/v SDS) i inkubowano w temperaturze 65 °C przez 15 min, wstrząsając przy 1000 obr/min w termomikserze (Thermo Scientific). Elucję powtórzono po raz drugi, a następnie do każdej próbki ChIP i wejściowej dodano 100 μL RNase Buffer (12 μL 5 M NaCl, 0.2 μL 30 mg/mL RNazy i 88 μL TE). Próbki inkubowano w 37 °C przez 20 min, a następnie dodawano 100 μL buforu z proteinazą K (2,5 μL 20 mg/mL proteinazą K, 5 μL 20% SDS i 92,5 μL TE) przez noc w 65 °C. Dodawano równą objętość roztworu fenol:chloroform i dokładnie mieszano. Mieszaninę przenoszono do probówek MaXtract High Density (Qiagen) i wirowano przez 8 min przy 15 000 rpm. Górną fazę przenoszono do nowych probówek i mieszano z 1,5 μL 20 mg/mL glikogenu, 30 μL 3M octanu sodu i 800 μL etanolu. Próbki inkubowano w temperaturze – 80 °C do momentu zamrożenia, a następnie odwirowywano przy 15 000 obr/min przez 30 min w temperaturze 4 °C. Supernatant usunięto, a osad przepłukano w 800 μL 70% lodowatego etanolu i wirowano przez 10 min w 4 °C przy 15 000 obr/min. Po dokładnym usunięciu etanolu, osady wysuszono na powietrzu i ponownie zawieszono w 30 μL 10 mM Tris o pH 8.
IP i wejściowe DNA były następnie oznaczane ilościowo przy użyciu fluorometru Qubit 3.0. Biblioteki przygotowywano przy użyciu zestawu KAPA HyperPrep Kit (KK8505) i sekwencjonowano przy użyciu Illumina NextSeq 500 do średniej głębokości 28 milionów odczytów na próbkę.
RNA-seq profiling
RNA izolowano z 3 milionów komórek na próbkę przy użyciu zestawu Bio-Rad Aurum™ Total RNA Mini Kit i kwantyfikowano przy użyciu zestawu Agilent RNA 6000 Nano Kit z bioanalizatorem Agilent. Biblioteki przygotowano przez zubożenie rRNA przy użyciu zestawu Illumina TruSeq® Stranded mRNA Library Prep Kit dla niskiego stężenia próbki wyjściowej i sekwencjonowano metodą sekwencjonowania pojedynczego końca na urządzeniu Illumina NextSeq 500 do średniej głębokości 18 milionów odczytów na próbkę.
Profilowanie metylacji DNA
Genomowe DNA wyizolowano przy użyciu zestawu AllPrep DNA/RNA Micro Kit (Qiagen). Aby ocenić status metylacji DNA w całym genomie, wykonaliśmy mRRBS . Po fluorometrycznej kwantyfikacji przy użyciu aparatu Qubit 3.0, trawiliśmy genomowe DNA enzymem restrykcyjnym MspI (New England Biolabs) i selekcjonowaliśmy fragmenty o długości około 100-250 par zasad przy użyciu kulek SPRI (MagBio Genomics). Otrzymane DNA poddano konwersji bisulfitowej przy użyciu zestawu EZ DNA Methylation-Lightning Kit (Zymo Research). Utworzyliśmy biblioteki z jednoniciowego DNA poddanego konwersji bisulfitowej przy użyciu zestawu Pico Methyl-Seq Library Prep Kit (Zymo Research), które następnie połączono do sekwencjonowania na aparacie Illumina NextSeq 500 przy użyciu zestawu odczynników NextSeq 500/550 V2 High Output (1 × 75 cykli) do minimalnej głębokości odczytu 50 milionów odczytów na próbkę.
Sekwencjonowanie całego genomu
Trzy miliony komórek z linii komórkowych lub próbek pobranych od pacjentów odmierzono i ponownie zawieszono w 1 mL roztworu Cell Lysis Solution (Qiagen) zmieszanego z 500 μg RNazy A. Reakcję lizy przeprowadzono w temperaturze 37 °C przez 15 min. W sumie do każdej próbki dodawano 333 μL Protein Precipitation Solution (Qiagen), którą następnie worteksowano, a następnie wirowano przy 2000×g przez 10 min. Supernatant mieszano z 1 mL izopropanolu do momentu wytrącenia nici DNA z roztworu. Po odrzuceniu supernatantu, osad DNA przemywano 1 mL 70% etanolu i odwirowywano przy 2000×g przez 1 min. Następnie wylano etanol, a osad suszono na powietrzu przez 15 min przed ponownym zawieszeniem w 50-100 μL roztworu do hydratacji DNA (Qiagen). DNA sekwencjonowano przy użyciu sekwencjonowania Illumina paired-end przy 30-krotnym pokryciu.
Immoprecypitacja
W sumie 100 milionów komórek dla każdej reakcji immunoprecypitacji granulowano i inkubowano w Buforze A (10 mM HEPES pH 8,0, 1,5 mM MgCl2, 10 mM KCl, 0,5 mM DTT) przez 10 min na lodzie. Następnie komórki lizowano po 12 uderzeniach 7-ml młynkiem do tkanek z luźnym tłuczkiem (Wheaton, 357542) i odwirowywano przy 2000 rpm przez 7 min. Osad jądrowy był ponownie zawieszony w 5 objętościach buforu TENT (50 mM Tris pH 7,5, 5 mM EDTA, 150 mM NaCl, 1% Triton X-100, 5 mM MgCl2) i traktowany benzonazą przez 30 minut przed 5 przejściami przez strzykawkę 25 g × 5/8 cala. Nierozpuszczalną frakcję usunięto po odwirowaniu przy 2000 obr/min przez 7 min i inkubowano przez noc z Dynabeads Protein G hybrydyzowanym z przeciwciałem. W sumie usunięto 2 miliony komórek dla danych wejściowych. Kulki i lizaty jąder były płukane 6 razy buforem TENT, a następnie eluowane w 0,1 M glicynie pH 2,5 z 100 mM Tris pH 8,0. Bufor do próbek NuPAGE LDS dodano do eluatów i danych wejściowych, które następnie inkubowano w temperaturze 70 °C przez 15 min przed analizą metodą western blot.
Gromadzenie danych publicznych
Publiczne dane CTCF ChIP-seq zostały pobrane z Cistrome Data Browser (dla plików szczytowych) i NCBI GEO (dla plików fastq, plik dodatkowy 2: Tabela S1). Dane Histone modification ChIP-seq zostały pobrane z NCBI GEO i ENCODE (dla plików bam). Publiczne dane RNA-seq w wielu typach komórek zostały pobrane z ENCODE (dla plików fastq). Dane dotyczące profilowania metylacji DNA zostały zebrane z ENCODE (dla plików bedMethyl) i NCBI GEO. Dane Hi-C zostały pobrane z NCBI GEO i ENCODE (dla plików fastq). Dane ATAC-seq zostały zebrane z NCBI GEO (dla plików fastq). Dane sekwencjonowania całego genomu dla próbek BRCA, COAD, LUAD i PRAD zostały pobrane z International Cancer Genome Consortium (ICGC) Data Portal. Szczegółowe informacje, w tym identyfikatory akcesji wszystkich publicznych zbiorów danych zebranych w tej pracy, można znaleźć w pliku dodatkowym 6: Tabela S5.
Opracowanie danych
Analiza danych ChIP-seq
Wyrównanie sekwencji dla danych ChIP-seq w plikach fastq przeprowadzono przy użyciu tego samego standardowego potoku analizy, który został użyty w Cistrome DB , w celu zapewnienia spójności i odtwarzalności. Wszystkie wyrównania genomowe danych sekwencyjnych zostały wykonane przy użyciu potoku Chilin z domyślnymi parametrami ($ chilin simple -p narrow -s hg38 –threads 8 -t IN.fq -i PRENAME -o OUTDIR). W skrócie, odczyty sekwencji zostały wyrównane do ludzkiego genomu referencyjnego (GRCH38/hg38) przy użyciu BWA ($ bwa aln -q 5 -l 32 -k 2 -t 8 INDEX IN.fq > PRENAME.sai $ bwa {samse | sampe} INDEX PRENAME.sai IN.fq > PRENAME.sam). Pliki sam były następnie konwertowane do plików bam za pomocą samtools ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Dla zestawów danych CTCF ChIP-seq zastosowano MACS2 do wywołania pików poniżej progu FDR 0,01 ($ macs2 callpeak –SPMR -B -q 0,01 –keep-dup 1 -g hs -t PRENAME.bam -n PRENAME –outidr OUTDIR). Zachowane zostały piki o wzbogaceniu co najmniej 4-krotnym. Pliki Bigwiggle zostały wygenerowane przy użyciu BEDTools i narzędzi UCSC ($ bedtools slop -i PRENAME.bdg -g CHROMSIZE -b 0|bedClip stdin CHROMSIZE PRENAME.bdg.clip $ LC_COLLATE=C sort -k1,1 -k2,2n PRENAME.bdg.clip > PRENAME.bdg.sort.clip $ bedGraphToBigWig PRENAME.bdg.sort.clip CHROMSIZE PRENAME.bw). Ostatecznie, tylko próbki CTCF ChIP-seq, które mają co najmniej 2000 pików, zostały włączone do dalszej analizy integracyjnej.
Analiza danychATAC-seq
Trim Galore został użyty do przycięcia surowych odczytów sekwencjonowania ($ trim_galore –nextera –phred33 –fastqc –paired R1.fq R2.fq -o OUTDIR). Odczyty zostały wyrównane do ludzkiego genomu referencyjnego (GRCH38/hg38) przy użyciu Bowtie2 ($ bowtie2 -p 10 -X 2000 -x INDEX -1 R1.fq -2 R2.fq -S PRENAME.sam). Pliki sam były następnie konwertowane do plików bam za pomocą samtools ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Bedtools posłużył do konwersji plików bam do formatu bed ($ bamToBed -i PRENAME.bam -bedpe > PRENAME_PE.bed). Odczyty zmapowane do mitochondrialnego DNA zostały odrzucone z dalszej analizy.
Analiza danych RNA-seq
Zestawy danych RNA-seq przetwarzano przy użyciu programu Salmon ($ salmon quant –gcBias -i INDEX -l A -p 8 {-1 R1.fq -2 R2.fq| -r IN.fq}. -o OUTDIR). Indeks transkryptomu został zbudowany na genomie referencyjnym człowieka (GRCH38/hg38). Szacunki obfitości transkryptów zostały podsumowane do poziomu genów przy użyciu pakietu „tximport” do analizy ekspresji różnicowej. DESeq2 został użyty do identyfikacji genów ulegających różnej ekspresji, a różne progi użyte w różnych analizach zostały odpowiednio wymienione w manuskrypcie.
Analiza danych Hi-C
Dane Hi-C zostały przetworzone przy użyciu HiC-Pro ($ HiC-Pro -i INDIR -o OUTDIR -c CONFIG -p). Mapy kontaktów zostały wygenerowane z rozdzielczością 5 kb. Surowe dane macierzowe zostały znormalizowane przy użyciu podejścia opisanego w Normalization of Chromatin Interactions.
Analiza danych metylacji DNA
Dane metylacji DNA (dla linii komórkowych T-ALL i pacjentów z T-ALL) zostały zdemultipleksowane przy użyciu bcl2fastq, a następnie przycięte o 10 par zasad z końca 5′ w celu usunięcia sekwencji startera i adaptora przy użyciu TrimGalore . Wyrównanie sekwencji do genomu referencyjnego GRCh38/hg38 i wywołania metylacji wykonano przy użyciu programu Bismark ($ bismark –multicore 8 –bowtie2 -q -N 1 INDEX INFILE.fq). Pliki pokrycia (counts) dla cytozyny w kontekście CpG zostały wygenerowane przy użyciu Bismark ($ bismark_methylation_extractor –multicore 8 –comprehensive –bedGraph INFILE_bismark_bt2.bam).
Analiza danych sekwencjonowania całego genomu
Mutacje zostały zidentyfikowane dla dwóch linii komórkowych T-ALL (Jurkat i CUTLL1) i dwóch próbek pacjentów T-ALL z danych sekwencjonowania całego genomu. Wyrównaliśmy sekwencje krótkich odczytów Illumina do ludzkiego genomu referencyjnego (GRCH38/hg38) używając BWA mem. Użyliśmy SAMBlaster do identyfikacji par niezgodnych, rozdzielenia odczytów i oznaczenia prawdopodobnych duplikatów PCR. Użyliśmy SAMBAMBA do konwersji SAM wyrównanych do formatu BAM, a samtools został użyty do sortowania tych wyrównanych w celu utworzenia pliku BAM odpowiadającego każdej próbce.
Użyliśmy VarDict do identyfikacji wariantów, które nakładały się na miejsca wiązania CTCF unii. Użyliśmy wszystkich domyślnych parametrów z wyjątkiem „-f 0.1”, który został użyty do identyfikacji wariantów, które były wspierane przez więcej niż 10% odczytów w tej lokalizacji. Dokonaliśmy anotacji wariantów za pomocą Variant Effect Predictor (VEP) i użyliśmy niestandardowych skryptów do identyfikacji wariantów, które wpływają na wiązanie TF.
Ponownie użyliśmy VarDict do identyfikacji wariantów w genach CTCF i NOTCH1 dla czterech próbek. Użyliśmy wszystkich domyślnych parametrów z wyjątkiem „-f 0.1”, który został użyty do identyfikacji wariantów, które były wspierane przez więcej niż 10% odczytów w tej lokalizacji. Dokonaliśmy anotacji wariantów przy użyciu Variant Effect Predictor (VEP), a następnie przefiltrowaliśmy je w celu zidentyfikowania mutacji, które (a) nie występowały w więcej niż 1% normalnej populacji ludzkiej, lub (b) posiadały wynik CADD dla deleteryjności > 20, lub (c) były obecne w bazie danych COSMIC.
Integracyjne modelowanie i analiza statystyczna
Identyfikacja repertuaru wiązania CTCF w ludzkim genomie
Dla CTCF ChIP-seq zebraliśmy łącznie 793 zestawy danych, w tym 787 publicznych zestawów danych i 6 zestawów danych wygenerowanych przez nas (Dodatkowy plik 2: Tabela S1). W sumie, 771 zestawów danych CTCF ChIP-seq z pikami powyżej 2000 zostało wykorzystanych w tym badaniu. Każdy zestaw danych może dostarczyć zidentyfikowanych przez MACS2 pików CTCF w zakresie od 2050 do 198 021, z medianą 46 451 i całkowitą liczbą 36 873 077 pików (Dodatkowy plik 1: Rys. S1a). Rozkład długości odstępów między sąsiednimi szczytami CTCF wszystkich 36 873 077 szczytów z 771 zestawów danych ma punkt przegięcia przy ~ 150 bp (plik dodatkowy 1: Fig. S1c) wskazujący granicę między tym samym miejscem wiązania a różnymi miejscami wiązania . Dlatego użyliśmy 150 bps jako punktu odcięcia do łączenia szczytów CTCF. W praktyce przedłużyliśmy ± 75 bps od każdego szczytu, aby wygenerować region 150-bp wyśrodkowany na szczycie, aby reprezentować każdy szczyt i połączyć wszystkie nakładające się regiony szczytowe, aby wygenerować zestaw unii miejsc wiążących CTCF, który zawiera 688 429 nienakładających się miejsc. Każdemu miejscu wiązania przypisano wynik zajętości CTCF, zdefiniowany jako suma zestawów danych ChIP-seq, które wykazują pik w danym miejscu. Odpowiednio, zdefiniowaliśmy częstotliwość zajętości jako stosunek wyniku zajętości do całkowitej liczby zestawów danych CTCF ChIP-seq. Aby dodatkowo zapewnić odporność zidentyfikowanych miejsc wiązania CTCF, wybraliśmy 285 467 miejsc o wysokiej wiarygodności z wynikiem zajętości ≥ 3 do dalszych analiz. Motywy CTCF w obrębie miejsc wiązania związków były przeszukiwane przez FIMO z macierzą Jaspar (ID: MA0139.1), z progiem wartości p wynoszącym 1e-4. Jeden motyw o najmniejszej wartości p został zachowany dla każdego miejsca wiążącego CTCF.
Identyfikacja konstytutywnych miejsc wiążących CTCF
Rozkład wyników zajętości wszystkich 285 467 miejsc wiążących CTCF (plik dodatkowy 1: Rys. S1d, niebieska krzywa) pokazuje, że większość miejsc wiążących CTCF występuje tylko w kilku zestawach danych, a liczba miejsc wiążących zmniejsza się wraz ze wzrostem wyniku zajętości, gdy wynik zajętości jest mały. Istnieją jednak miejsca wiążące CTCF, które są silnie konserwowane w prawie wszystkich zestawach danych (np. miejsca wiążące o zajętości większej niż 600). Używamy funkcji power law do dopasowania krzywej rozkładu (niebieskiej) pokazanej w pliku dodatkowym 1: Rys. S1d, aby określić wartość graniczną dla konstytutywnych miejsc CTCF. Oznaczamy Oi jako liczbę obserwowanych miejsc wiążących CTCF z wynikiem zajętości równym i, a Ei jako liczbę oczekiwanych miejsc CTCF z wynikiem zajętości równym i. Dopasowanie potęgowe do danych Oi można opisać jako (plik dodatkowy 1: Rys. S1d, kolor zielony):
Definiujemy cutoff A dla konstytutywnych miejsc wiązania CTCF jako:
Identyfikacja specyficznych dla nowotworu zyskanych/utraconych miejsc wiążących CTCF
Użyliśmy następujących 2 kryteriów do identyfikacji specyficznych dla nowotworu utraconych miejsc wiążących CTCF: (1) miejsce wiązania CTCF powinno mieć niższą częstotliwość zajmowania dla zbiorów danych tego typu nowotworu w porównaniu do częstotliwości zajmowania dla wszystkich zbiorów danych oraz (2) poziom wiązania CTCF (określony ilościowo jako znormalizowana liczba odczytów ChIP-seq) w danym miejscu jest niższy w zbiorach danych nowotworów niż w innych zbiorach danych. Dla pozyskanych miejsc CTCF zastosowaliśmy odwrotny zestaw kryteriów. W skrócie, dla każdego miejsca wiążącego CTCF w każdym typie nowotworu, obliczono wynik zajętości w zbiorach danych nowotworu wraz z jego wynikiem zajętości we wszystkich 771 zbiorach danych. Poziomy wiązania CTCF uzyskano z normalizowanej macierzy liczby odczytów, w której najpierw obliczono liczbę odczytów ChIP-seq (RPKM) dla miejsc wiązania CTCF we wszystkich zestawach danych, a następnie dokonano normalizacji kwantylowej. Użyliśmy niesparowanego, dwuogonowego testu t-Studenta do ilościowego określenia różnicy poziomów wiązania pomiędzy różnymi grupami zbiorów danych, a wartość p została następnie skorygowana przy użyciu procedury Benjamini-Hochberga. Ponadto, wyniki zajętości wiązania i poziomy wiązania zostały porównane między zestawami danych nowotworowych i zestawami danych z dopasowanych normalnych tkanek lub typów komórek, aby uwzględnić potencjalny czynnik zakłócający specyficzność tkanki, a nie specyficzność nowotworu. Szczegółowe kryteria identyfikacji specyficznych dla nowotworów miejsc wiązania CTCF opisano poniżej:
-
Swoiste dla nowotworów utracone miejsca wiązania CTCF: (1) częstość zajętości ≤ 0,2 w zbiorach danych dotyczących nowotworów; (2) częstość zajętości ≥ 0,7 w 771 zbiorach danych; (3) częstość zajętości ≥ 0.5 (z occupancy score ≥ 2) w dopasowanych zbiorach danych tkanki prawidłowej; (4) poziomy CTCF są niższe w raku w porównaniu do wszystkich innych zbiorów danych (wynik statystyczny < 0), (5) poziomy CTCF są niższe w raku w porównaniu do dopasowanych zbiorów danych tkanki prawidłowej (wynik statystyczny < 0), (6) uśrednione sygnały wiązania CTCF (RPKM) < 5 w zbiorach danych raka.
-
Cancer-specific gained CTCF binding sites: (1) częstość zajętości ≥ 0,5 (z occupancy score ≥ 2) w zbiorach danych nowotworowych, (2) częstość zajętości ≤ 0,2 w 771 zbiorach danych, (3) occupancy score = 0 w dopasowanych zbiorach danych tkanek prawidłowych, (4) poziomy CTCF są istotnie wyższe w nowotworach w porównaniu do wszystkich pozostałych zbiorów danych (FDR ≤ 0.01), (5) poziomy wiązania CTCF są znacząco wyższe w nowotworach w porównaniu z dopasowanymi zestawami danych tkanki prawidłowej (FDR ≤ 0.01), (6) uśrednione sygnały wiązania CTCF (RPKM) > 2 w zestawach danych nowotworowych.
Szczególne zyskane i utracone miejsca wiązania CTCF dla każdego typu nowotworu przedstawiono w pliku dodatkowym 4: Tabela S3.
Kwantyfikacja zróżnicowanej dostępności chromatyny
Użyliśmy przetworzonych danych z Ref. które zawierają macierz znormalizowanych liczebności insercji ATAC-seq w ramach zestawu szczytowego TCGA pan-cancer, aby ocenić zróżnicowaną dostępność chromatyny wokół miejsc wiązania CTCF. Dla każdego typu nowotworu spośród BRCA, CRC, LUAD i PRAD, pan-nowotworowe piki ATAC-seq, które pokrywają się ze zidentyfikowanymi, specyficznymi dla danego typu nowotworu utraconymi lub zyskanymi miejscami wiązania CTCF, zostały użyte do dalszych analiz. Wynik różnicowy ATAC-seq dla każdego piku został określony ilościowo jako zmiana fałdowa średniej znormalizowanej liczby insercji ATAC-seq z próbek pacjentów w odpowiednim typie nowotworu w porównaniu z pacjentami w innych typach nowotworów, a wynik różnicowy ATAC-seq został następnie przypisany do piku nakładającego się na miejsce wiązania CTCF.
Dla spójności, zastosowaliśmy to samo podejście stosowane do danych TCGA ATAC-seq do analizy zebranych danych ATAC-seq z linii komórkowej T-ALL Jurkat i normalnych komórek CD4+ T. Macierz danych została wygenerowana przy użyciu surowych odczytów ATAC-seq dla miejsc wiązania CTCF dla wszystkich zestawów danych Jurkat i komórek T. Normalizacja kwantylowa została zastosowana na matrycy skalowanej log2 (pseudo count = 5). Wynik różnicowy ATAC-seq był mierzony jako zmiana fałdowa uśrednionych znormalizowanych zliczeń ATAC-seq między zestawami danych Jurkat versus komórki CD4+ T w każdym miejscu wiązania CTCF.
Normalizacja oddziaływań chromatynowych
Gdy dana jest mapa kontaktów Hi-C A = {aij}, wynik aij odzwierciedla zmapowane odczyty pomiędzy dwoma regionami genomowymi i oraz j. Załóżmy, że rozmiar bloku wynosi 5 kb, regiony i oraz j będą miały odległość genomową ∣i – j ∣ × 5kb. Ponieważ prawdopodobieństwo kontaktu między dwoma binami maleje wraz ze wzrostem odległości genomowej, znormalizowaliśmy mapę kontaktów w następujący sposób: dla każdej danej odległości genomowej dk = k × 5kb, określamy współczynnik normalizacji \( {overline{S}}_{d_k} \) jako uśrednione interakcje między wszystkimi parami binów o tej samej odległości genomowej dk w tym samym chromosomie, np, \( {overline{S}}_{d_k}=left({sum}_{j-i=k}{a}_{ij}}right)/n \), gdzie n jest całkowitą liczbą par bin o odległości dk. Wynik interakcji aij pomiędzy dwoma binami o odległości dk został następnie znormalizowany przez \( {overline{S}_{d_k} \) jako \( {a}_{ij}^{prime }={a}_{ij}/{overline{S}_{d_k} \). Używając tego podejścia, znormalizowaliśmy macierz A do postaci ^( A^{ij}^{prime }={left}{a}^{prime}} \) w obrębie każdego chromosomu.
Detekcja zróżnicowanych interakcji chromatynowych
Zaznaczyliśmy znormalizowane mapy kontaktów Hi-C w zbiorze danych nowotworowych i normalnym zbiorze danych jako C = {cij} i N = {nij}, odpowiednio. Dla danego miejsca wiązania CTCF x (o współrzędnej xc) i wstępnie zdefiniowanej odległości genomowej L, interakcje chromatynowe pomiędzy x i jego pobliskimi, nienakładającymi się na siebie 5-kilobinami o odległości genomowej do L są zbierane odpowiednio z C i N. Konkretnie, wyniki interakcji pomiędzy x i jego pobliskimi 5-kb pinami w C są zbierane jako IC = {cij} , podczas gdy i lub j równa się ⌊xc/5kb⌋, a 0 < (j – i) × 5kb ≤ L. Podobnie, wyniki interakcji pomiędzy x i jego pobliskimi 5-kb binami w N zostały zebrane jako IN = {nij}. Następnie zastosowano sparowany dwuwargowy test t-Studenta na IC i IN, aby określić ilościowo zróżnicowanie interakcji między komórkami nowotworowymi i normalnymi otaczającymi miejsce wiązania CTCF x.
Asocjacja wiązania CTCF z ekspresją genów
W sumie wybrano 54 typy komórek, dla których zarówno dane CTCF ChIP-seq, jak i dane RNA-seq są publicznie dostępne (plik dodatkowy 6: Tabela S5) w celu zbadania asocjacji między wiązaniem CTCF a ekspresją genów dla każdej pary CTCF-genu w tym samym chromosomie. Aby uzyskać poziom wiązania CTCF, wygenerowano macierz liczby odczytów przy użyciu odczytów na kilobazę na milion (RPKM) w miejscach wiązania CTCF z danych ChIP-seq. Macierz liczby odczytów została przeskalowana pierwiastkiem kwadratowym z RPKM, a następnie poddana normalizacji kwantylowej. Poziom ekspresji genów był mierzony dla każdego genu przy użyciu pierwiastka kwadratowego transkryptów na milion (TPM) z danych RNA-seq. Dla każdej pary CTCF-gen określiliśmy ilościowo związek między miejscem CTCF a genem we wszystkich 54 typach komórek, używając współczynnika korelacji R między znormalizowanym poziomem wiązania CTCF a ekspresją genu (ryc. 3a). Pary CTCF-geny uznano za „wysoko skorelowane” z R2 większym niż 0,25, np, współczynnik korelacji większy niż 0,5 lub mniejszy niż – 0,5, a wysoko skorelowane pary CTCF-geny stanowią 1,3% wszystkich par CTCF-geny (plik dodatkowy 1: Fig. S8a).
Identyfikacja konstytutywnych domen chromatyny związanych z CTCF
Dla każdego miejsca wiązania CTCF zdefiniowaliśmy jego powiązaną domenę chromatyny jako region genomowy, który (1) obejmuje to konkretne miejsce wiązania CTCF, (2) jest ograniczony przez parę konstytutywnych miejsc wiązania CTCF z motywami o przeciwnych orientacjach i (3) zajmuje minimum 100 kb i maksimum 1 MB regionu po każdej stronie miejsca wiązania CTCF. Rysunek 3b zawiera schemat tego, jak zdefiniowano konstytutywne domeny chromatyny związane z CTCF.
Detekcja zmian metylacji DNA wokół miejsc wiązania CTCF
Zmiany metylacji DNA wykryto w regionie 300-bp wyśrodkowanym na każdym miejscu wiązania CTCF. Zachowano regiony z co najmniej 3 CpGs pokrytymi przez co najmniej 5 odczytów (≥ 5×) w obu liniach komórek nowotworowych i odpowiadających im prawidłowych tkankach. Region 300-bp został wykryty jako zróżnicowanie metylacji, jeśli uśrednione zróżnicowane poziomy metylacji wszystkich CpGs (≥ 5×) w tym regionie były większe niż 20% .
Detekcja wskaźnika mutacji i zróżnicowanego wyniku motywu
Dla każdego miejsca wiązania CTCF, surowa liczba mutacji została obliczona jako występowanie zdarzeń mutacji we wszystkich próbkach/pacjentach w każdej pojedynczej parze zasad w regionie 400-bp wyśrodkowanym na miejscu wiązania CTCF. Wskaźnik mutacji dla grupy miejsc wiążących CTCF został obliczony jako uśredniona liczba mutacji w stosunku do liczby miejsc wiążących CTCF dla każdej pary zasad w regionie 400-bp.
Wynik motywu został zmierzony przez ocenę macierzy wagowej pozycji CTCF (Jaspar , ID macierzy: MA0139.1) do 19-bp sekwencji DNA wyśrodkowanej na motywie CTCF lub miejscu wiążącym CTCF przy użyciu logicznych współczynników prawdopodobieństwa (z częstotliwością nukleotydów tła jak dla A,C,G,T). Różnicowy wynik motywu został obliczony przez porównanie wyników motywu dla sekwencji referencyjnej i zmutowanej.
Analiza motywów sekwencji DNA
Analizę wzbogacania motywów sekwencji DNA przeprowadzono przy użyciu MDSeqPos (wersja 1.0.0) na Cistrome z domyślnymi parametrami (-cisrome -Homo Sapien lub Mus musculus). Analizy motywów de novo przeprowadzono przy użyciu programu HOMER (wersja 4.10) z modułem findmotifs.pl oraz MEME (wersja 5.1.1) z następującymi parametrami: meme -dna -mod zoops -maxw 20 -evt -0.01.
Identyfikacja wewnątrzdomenowych regionów różnie oddziałujących z CTCF
Dla danego zestawu miejsc wiążących CTCF, dla każdego z miejsc wiążących CTCF zebrano zmiany interakcji chromatynowych między miejscem CTCF a każdym z jego wewnątrzdomenowych nienakładających się koszy, mierzonych na podstawie znormalizowanych map kontaktów Hi-C w komórkach nowotworowych w stosunku do dopasowanych komórek prawidłowych (plik dodatkowy 1: Fig. S14b). Regiony o zmniejszonych interakcjach (log2 FC < -1, uśredniony log2 interakcji > 0) ze specyficznymi dla raka utraconymi miejscami wiązania CTCF i regiony o zwiększonych interakcjach (log2 FC > 1, uśredniony log2 interakcji > 0) ze specyficznymi dla raka zyskanymi miejscami wiązania CTCF zostały wykorzystane do analizy wzbogacania czynników transkrypcyjnych (TF) w dół.
Analiza wzbogacania czynników transkrypcyjnych
Do analizy wzbogacania TF użyto poprawionej wersji algorytmu BART. W skrócie, kolekcja związkowych miejsc wrażliwych na nadwrażliwość DNase I (UDHS) została uprzednio opracowana jako repertuar wszystkich kandydujących elementów cis-regulacyjnych w ludzkim genomie, a 7032 zestawy danych ChIP-seq zostały zebrane dla 883 TF , przy czym każdy TF miał jeden lub więcej zestawów danych ChIP-seq z wielu typów komórek lub warunków. Dla każdego TF na UDHS wygenerowano profil binarny wskazujący, czy TF ma co najmniej jeden pik z któregokolwiek ze swoich zestawów danych ChIP-seq zlokalizowanych w każdym z UDHS. Analiza wzbogacenia wiązania została zastosowana dla każdego TF poprzez porównanie wiązania TF na podzbiorze UDHS pokrywających się z wybranymi regionami genomowymi w porównaniu do wiązania TF na UDHS. Wartość p została uzyskana przy użyciu dokładnego testu Fishera z dwoma ogonami.