Experimentelles Vorgehen
Patienten-Xenotransplantation und Zellkultur
Die menschlichen T-ALL-Zelllinien sind CUTLL1 (Geschenk von Adolfo Ferrando, Columbia University) und JURKAT (American Type Culture Collection (ATCC), Manassas, VA, #CCL-119). Die Zellen wurden in RPMI1640-Medium mit L-Glutamin und 25 mM HEPES (Corning) kultiviert, das mit 10 % hitzeinaktiviertem fötalem Rinderserum (Sigma-Aldrich), 10 U/mL Penicillin-Streptomycin (Gibco) und 1× GlutaMAX (Gibco) in einem befeuchteten Inkubator bei 37 °C und 5 % CO2 ergänzt wurde. Die Zellen werden regelmäßig mit dem Lonza Walkersville MycoAlert Mycoplasma Detection Kit auf das Vorhandensein von Mykoplasmen getestet (letzter Test im Januar 2020). Die Zelllinien werden für maximal 20 Passagen in Kultur gehalten und mittels Short-Tandem-Repeats-Profiling (JURKAT) oder mittels PCR zum Nachweis der TCRb-NOTCH1-Translokation (TCRBJ2S4CUTLL1F:5′-GGACCCGGCTCTCAGTGCT-3′, NOTCH1CUTTL1R:5′-TCCCGCCCTCCAAAATAAGG-3′) authentifiziert. Die letzte Zellauthentifizierung wurde im Februar 2020 durchgeführt. Humane CD4+ T-Zellen wurden von AllCells erworben. Humane Primärproben wurden mit informierter Zustimmung entnommen und unter der Aufsicht des Institutional Review Board der Universität Padua, der Associazone Italiana di Ematologia e Oncologia Pediatrica und der Berlin-Frankfurt-Münster (AIEOP-BFM) ALL 2000/2006 pädiatrischen klinischen Studien analysiert. Von allen Patienten wurde bei Studienbeginn gemäß der Deklaration von Helsinki eine informierte Zustimmung zur Verwendung von Restmaterial für Forschungszwecke eingeholt.
Antikörper und Reagenzien
Western Blots wurden mit den folgenden Antikörpern durchgeführt: Actin und CTCF von Millipore Sigma (Klon C4; 07-729) und gespaltenes NOTCH1 (Val1744) von Cell Signaling Technology (4147). ChIP-seq wurden mit den folgenden Antikörpern durchgeführt: CTCF von Millipore Sigma (07-729), H3K27Ac (8173S) und H3K27me3 (9733S) von Cell Signaling Technology und H3K4me1 (07-473) von Millipore.
In situ Hi-C
In situ Hi-C wurde an CD4+ T-Zellen, Jurkat, CUTLL1 und Patienten-Xenografts wie zuvor beschrieben durchgeführt. Die Zellen wurden 10 Minuten lang mit 1% Formaldehyd bei Raumtemperatur vernetzt. Pro Hi-C-Reaktion wurden 5 Millionen Zellen lysiert und die Zellkerne permeabilisiert. Die DNA wurde mit MboI von New England Biolabs (R0147M) verdaut. Die verdauten Fragmente wurden mit biotinyliertem d-ATP von Jena Bioscience (NU-835-BIO14-S) markiert und ligiert. Nach der RNase-Behandlung und der Behandlung mit Proteinase K zur Umkehrung der Querverbindungen wurden die Kerne mit einem Covaris E220 beschallt, um eine durchschnittliche Fragmentlänge von 400 bp zu erhalten. Streptavidin-Beads von Thermo Fisher Scientific (65001) wurden verwendet, um Biotin-markierte Fragmente herunterzuziehen. Nach der Aufreinigung und Isolierung der DNA wurden die endgültigen Bibliotheken mit dem NEBNext® Ultra™ II DNA Library Prep Kit für Illumina® vorbereitet und mittels Paired-End-Sequenzierung mit einer Leselänge von 150 bp auf einem Illumina HiSeq 2500 sequenziert, um durchschnittlich 400 Millionen Reads pro Probe zu erzeugen.
ChIP-seq-Profiling
CD4+ T-Zellen, Jurkat, CUTLL1 und Xenotransplantate von Patienten wurden mit 1% Formaldehyd und 1% fötalem Rinderserum in PBS für 10 Minuten bei Raumtemperatur vernetzt. Die Reaktion wurde mit 0,2 M Glycin bei Raumtemperatur für 5 Minuten gequencht. Die Zellen wurden dann mit PBS gewaschen und pelletiert.
Für CTCF-ChIPs wurde die Immunpräzipitation nach einem zuvor beschriebenen Protokoll durchgeführt. Ein Pellet mit 50 Millionen Zellen wurde mit 5 mL Lysepuffer (50 mM HEPES-KOH, pH 7,5, 140 mM NaCl, 1 mM EDTA, 10 % Glycerin, 0,5 % NP-40, 0,25 % Triton X-100) für 10 Minuten bei 4 °C lysiert. Die Kerne wurden 7 Minuten lang bei 1350×g pelletiert und in 10 mM Tris pH 8, 1 mM EDTA und 0,1% SDS resuspendiert. Das Chromatin wurde mit einem Covaris E220-System auf eine durchschnittliche Fragmentlänge von 400 bp geschert und 10 Minuten lang bei 15.000 U/min geschleudert, um unlösliches Chromatin und Trümmer zu entfernen. Der Überstand wurde mit 20 μl Dynabeads Protein G 30 Minuten lang inkubiert, bevor die Beads verworfen wurden. Ein Prozent des Gesamtvolumens wurde als Input aufbewahrt, der Rest wurde über Nacht mit dem Anti-CTCF-Antikörper inkubiert. Insgesamt wurden 100 μL Dynabeads Protein G für 2 h zugegeben. Die gebundenen Fragmente wurden zweimal mit 1 mL eines salzarmen Puffers (20 mM Tris-HCl pH 8,0, 150 mM NaCl, 2 mM EDTA, 1 % w/v Triton X-100 und 0.1% w/v SDS), einmal mit Puffer mit hohem Salzgehalt (20 mM Tris-HCl pH 8,0, 500 mM NaCl, 2 mM EDTA, 1% w/v Triton X-100 und 0,1% w/v SDS), einmal mit Lithiumchloridpuffer (10 mM Tris-HCl pH 8.0, 250 mM LiCl, 1 mM EDTA, 1% Gew./Vol. NP-40 und 1% Gew./Vol. Desoxycholsäure) und zweimal mit TE (10 mM Tris pH 8, 1 mM EDTA).
Für Histon-ChIPs wurden die Zellen in 375 μL Kern-Inkubationspuffer (15 mM Tris pH 7.5, 60 mM KCl, 150 mM NaCl, 15 mM MgCl2, 1 mM CaCl2, 250 mM Saccharose, 0,3% NP-40, 1 mM NaV, 1 mM NaF und 1 EDTA-freie Proteaseinhibitor-Tablette (Roche)/10 mL in H2O) für 10 min auf Eis. Die Kerne wurden einmal mit Verdauungspuffer (10 mM NaCl, 10 mM Tris pH 7,5, 3 mM MgCl2, 1 mM CaCl2, 1 mM NaV, 1 mM NaF und 1 EDTA-freie Proteaseinhibitor-Tablette (Roche)/10 mL in H2O) gewaschen und in 57-μL Verdauungspuffer mit 4,5 Einheiten MNase (USB) für 1 h bei 37 °C resuspendiert. Die MNase-Aktivität wurde 10 Minuten lang auf Eis durch Zugabe von EDTA bis zu einer Endkonzentration von 20 mM gequencht. Die Zellkerne wurden pelletiert und in 300 μl Zellkern-Lysepuffer (50 mM Tris-HCl pH 8,0, 10 mM EDTA pH 8,0, 1 % SDS, 1 mM NaV, 1 mM NaF und 1 EDTA-freie Proteaseinhibitor-Tablette (Roche)/10 mL in H2O) resuspendiert, bevor sie 5 Minuten lang mit einem Bioruptor Pico (Diagenode) beschallt wurden (30 s an, 30 s aus). Das Lysat wurde 5 Minuten lang bei maximaler Geschwindigkeit zentrifugiert, um Trümmer zu entfernen. Neun Volumina IP-Verdünnungspuffer (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 und 1 EDTA-freie Proteaseinhibitor-Tablette (Roche)/10 mL in H2O) wurden zum Überstand hinzugefügt. Insgesamt wurden 50 μl Dynabeads Protein G zugegeben und die Probe bei 4 °C 30 Minuten lang unter Rotation inkubiert. Ein Prozent der Probe wurde als Input aufbewahrt, und die restliche Probe wurde auf 3 Röhrchen aufgeteilt. Insgesamt wurden 50 μl Dynabeads Protein G, die mit 15 μl des entsprechenden Antikörpers konjugiert waren, in jedes Röhrchen gegeben, bevor die Inkubation über Nacht bei 4 °C unter Rotation erfolgte. Die an die Beads gebundenen Komplexe wurden jeweils 5 Minuten lang in 1 ml Puffer mit niedrigem Salzgehalt, Puffer mit hohem Salzgehalt, LiCl-Puffer und zweimal mit TE gewaschen.
Um die an die Beads gebundenen Komplexe zu eluieren, wurden die Beads in 50 μl Elutionspuffer (100 mM NaHCO3, 1 % w/v SDS) resuspendiert und 15 Minuten lang bei 65 °C unter Schütteln bei 1000 U/min in einem Thermomixer (Thermo Scientific) inkubiert. Die Elution wurde ein zweites Mal wiederholt, und dann wurden 100 μL RNase-Puffer (12 μL 5 M NaCl, 0,2 μL 30 mg/mL RNase und 88 μL TE) zu jeder ChIP- und Input-Probe gegeben. Die Proben wurden 20 Minuten lang bei 37 °C inkubiert, gefolgt von der Zugabe von 100 μL Proteinase-K-Puffer (2,5 μL 20 mg/mL Proteinase K, 5 μL 20% SDS und 92,5 μL TE) über Nacht bei 65 °C. Ein gleiches Volumen an Phenol:Chloroform-Lösung wurde hinzugefügt und gründlich gemischt. Das Gemisch wurde in MaXtract High Density-Röhrchen (Qiagen) überführt und 8 Minuten lang bei 15.000 U/min zentrifugiert. Die obere Phase wurde in neue Röhrchen überführt und mit 1,5 μL 20 mg/mL Glykogen, 30 μL 3M Natriumacetat und 800 μL Ethanol gemischt. Die Proben wurden bei – 80 °C bebrütet, bis sie gefroren waren, und dann bei 15.000 U/min für 30 Minuten bei 4 °C zentrifugiert. Der Überstand wurde entfernt und die Pellets wurden in 800 μL 70%igem, eiskaltem Ethanol gewaschen und 10 min bei 4 °C und 15.000 U/min geschleudert. Nach sorgfältiger Entfernung des Ethanols wurden die Pellets an der Luft getrocknet und in 30 μl 10 mM Tris bei pH 8 resuspendiert.
IP und Input-DNA wurden dann mit einem Qubit 3.0 Fluorometer quantifiziert. Bibliotheken wurden mit dem KAPA HyperPrep Kit (KK8505) vorbereitet und mit einem Illumina NextSeq 500 bis zu einer durchschnittlichen Tiefe von 28 Millionen Reads pro Probe sequenziert.
RNA-seq profiling
RNA wurde aus 3 Millionen Zellen pro Probe mit dem Bio-Rad Aurum™ Total RNA Mini Kit isoliert und mit dem Agilent RNA 6000 Nano Kit mit dem Agilent Bioanalyzer quantifiziert. Bibliotheken wurden durch rRNA-Abreicherung mit dem Illumina TruSeq® Stranded mRNA Library Prep Kit für eine niedrige Konzentration der Ausgangsprobe vorbereitet und durch Single-End-Sequenzierung auf einem Illumina NextSeq 500 mit einer durchschnittlichen Tiefe von 18 Millionen Reads pro Probe sequenziert.
DNA-Methylierungsprofilierung
Genomische DNA wurde mit dem AllPrep DNA/RNA Micro Kit (Qiagen) isoliert. Zur Bewertung des genomweiten DNA-Methylierungsstatus führten wir mRRBS durch. Nach der fluorometrischen Quantifizierung mit einem Qubit 3.0-Instrument verdauten wir genomische DNA mit dem Restriktionsenzym MspI (New England Biolabs) und führten eine Größenselektion für Fragmente mit einer Länge von etwa 100-250 Basenpaaren unter Verwendung von SPRI-Beads (Solid Phase Reversible Immobilization) (MagBio Genomics) durch. Die resultierende DNA wurde mit dem EZ DNA Methylation-Lightning Kit (Zymo Research) bisulfitiert. Wir erstellten Bibliotheken aus Bisulfit-konvertierter einzelsträngiger DNA mit dem Pico Methyl-Seq Library Prep Kit (Zymo Research), die dann für die Sequenzierung auf einem Illumina NextSeq 500-Gerät unter Verwendung des NextSeq 500/550 V2 High Output-Reagenzienkits (1 × 75 Zyklen) mit einer Mindestlesetiefe von 50 Millionen Reads pro Probe gepoolt wurden.
Ganzgenomsequenzierung
Drei Millionen Zellen aus Zelllinien oder Patientenproben wurden pelletiert und in 1 mL Cell Lysis Solution (Qiagen), gemischt mit 500 μg RNase A, resuspendiert. Die Lyse-Reaktion wurde 15 Minuten lang bei 37 °C durchgeführt. Insgesamt wurden 333 μl Proteinpräzipitationslösung (Qiagen) zu jeder Probe gegeben, die dann mit dem Schüttler geschüttelt und anschließend 10 Minuten lang bei 2000×g zentrifugiert wurde. Der Überstand wurde mit 1 mL Isopropanol gemischt, bis die DNA-Stränge aus der Lösung ausfielen. Nach Verwerfen des Überstands wurde das DNA-Pellet mit 1 mL 70%igem Ethanol gewaschen und 1 Minute lang bei 2000×g zentrifugiert. Anschließend wurde das Ethanol abgegossen und das Pellet 15 Minuten lang an der Luft getrocknet, bevor es in 50 bis 100 μl DNA-Hydrierungslösung (Qiagen) resuspendiert wurde. Die DNA wurde mit Paired-End-Sequenzierung von Illumina mit 30-facher Abdeckung sequenziert.
Immunpräzipitation
Insgesamt 100 Millionen Zellen für jede Immunpräzipitationsreaktion wurden pelletiert und in Puffer A (10 mM HEPES pH 8,0, 1,5 mM MgCl2, 10 mM KCl, 0,5 mM DTT) für 10 Minuten auf Eis inkubiert. Die Zellen wurden dann nach 12 Schlägen mit einem 7-mL-Gewebezerkleinerer mit losem Stößel (Wheaton, 357542) lysiert und 7 Minuten lang bei 2000 rpm zentrifugiert. Die Kernpellets wurden in 5 Volumina TENT-Puffer (50 mM Tris pH 7,5, 5 mM EDTA, 150 mM NaCl, 1 % Triton X-100, 5 mM MgCl2) resuspendiert und 30 Minuten lang mit Benzonase behandelt, bevor sie fünfmal durch eine 25 g × 5/8 Zoll-Spritze gepresst wurden. Die unlösliche Fraktion wurde nach 7-minütiger Zentrifugation bei 2000 rpm abgetrennt und über Nacht mit Dynabeads Protein G, die mit einem Antikörper hybridisiert waren, inkubiert. Insgesamt wurden 2 Millionen Zellen für die Eingabe entnommen. Beads und Zellkernlysate wurden 6-mal mit TENT-Puffer gewaschen und anschließend in 0,1 M Glycin pH 2,5 mit 100 mM Tris pH 8,0 eluiert. NuPAGE LDS-Probenpuffer wurde zu Eluaten und Inputs hinzugefügt, die dann bei 70 °C für 15 Minuten vor der Analyse durch Western Blot inkubiert wurden.
Öffentliche Datensammlung
Öffentliche CTCF ChIP-seq-Daten wurden von Cistrome Data Browser (für Peak-Dateien) und NCBI GEO (für fastq-Dateien, Additional file 2: Table S1) gesammelt. Histon-Modifikation ChIP-seq-Daten wurden von NCBI GEO und ENCODE (für bam-Dateien) gesammelt. Öffentliche RNA-seq-Daten in mehreren Zelltypen wurden von ENCODE gesammelt (für fastq-Dateien). DNA-Methylierungsprofilierungsdaten wurden von ENCODE (für bedMethyl-Dateien) und NCBI GEO gesammelt. Hi-C-Daten wurden von NCBI GEO und ENCODE (für fastq-Dateien) gesammelt. ATAC-seq-Daten wurden von NCBI GEO (für fastq-Dateien) gesammelt. Whole Genome Sequencing-Daten für BRCA-, COAD-, LUAD- und PRAD-Proben wurden vom Datenportal des International Cancer Genome Consortium (ICGC) gesammelt. Detaillierte Informationen, einschließlich der Zugangs-IDs aller öffentlichen Datensätze, die in dieser Arbeit gesammelt wurden, finden Sie in Zusatzdatei 6: Tabelle S5.
Datenverarbeitung
ChIP-seq-Datenanalyse
Das Sequenz-Alignment für ChIP-seq-Daten in fastq-Dateien wurde mit der gleichen Standard-Analyse-Pipeline durchgeführt, die auch in Cistrome DB verwendet wird, um Konsistenz und Reproduzierbarkeit zu gewährleisten. Alle Sequenzdaten wurden mit der Chilin-Pipeline mit Standardparametern ausgerichtet ($ chilin simple -p narrow -s hg38 –threads 8 -t IN.fq -i PRENAME -o OUTDIR). Die Sequenzlesungen wurden mit Hilfe von BWA an das menschliche Referenzgenom (GRCH38/hg38) angeglichen ($ 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). Die Sam-Dateien wurden dann mit samtools in bam-Dateien umgewandelt ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Für CTCF ChIP-seq-Datensätze wurde MACS2 verwendet, um Peaks unter dem FDR-Schwellenwert von 0,01 zu ermitteln ($ macs2 callpeak –SPMR -B -q 0,01 –keep-dup 1 -g hs -t PRENAME.bam -n PRENAME –outidr OUTDIR). Peaks mit einer Faltenanreicherung von mindestens 4 wurden beibehalten. Bigwiggle-Dateien wurden mit BEDTools und UCSC-Tools erzeugt ($ 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). Schließlich wurden nur die CTCF-ChIP-seq-Proben mit mindestens 2000 Peaks in die nachgelagerte integrative Analyse einbezogen.
ATAC-seq-Datenanalyse
Trim Galore wurde zum Trimmen der Rohsequenzierungs-Reads verwendet ($ trim_galore –nextera –phred33 –fastqc –paired R1.fq R2.fq -o OUTDIR). Die Reads wurden mit Hilfe von Bowtie2 am menschlichen Referenzgenom (GRCH38/hg38) ausgerichtet ($ bowtie2 -p 10 -X 2000 -x INDEX -1 R1.fq -2 R2.fq -S PRENAME.sam). Die Sam-Dateien wurden dann mit samtools in bam-Dateien umgewandelt ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Bedtools wurde zur Konvertierung von bam-Dateien in das bed-Format verwendet ($ bamToBed -i PRENAME.bam -bedpe > PRENAME_PE.bed). Reads, die der Mitochondrien-DNA zugeordnet waren, wurden von der nachgeschalteten Analyse ausgeschlossen.
RNA-seq-Datenanalyse
RNA-seq-Datensätze wurden mit Salmon verarbeitet ($ salmon quant –gcBias -i INDEX -l A -p 8 {-1 R1.fq -2 R2.fq| -r IN.fq} -o OUTDIR). Der Transkriptom-Index wurde auf dem menschlichen Referenzgenom (GRCH38/hg38) aufgebaut. Die Abundanzschätzungen auf Transkriptionsebene wurden mit Hilfe des „tximport“-Pakets für die differenzielle Expressionsanalyse auf Genebene zusammengefasst. DESeq2 wurde verwendet, um differenziell exprimierte Gene zu identifizieren, und die in den verschiedenen Analysen verwendeten Schwellenwerte wurden im Manuskript entsprechend aufgeführt.
Hi-C-Datenanalyse
Hi-C-Daten wurden mit HiC-Pro verarbeitet ($ HiC-Pro -i INDIR -o OUTDIR -c CONFIG -p). Es wurden Kontaktkarten mit einer Auflösung von 5 kb erstellt.
DNA-Methylierungsdatenanalyse
DNA-Methylierungsdaten (für T-ALL-Zelllinien und T-ALL-Patienten) wurden mit bcl2fastq demultiplexiert, gefolgt von einem Trimmen von 10 Basenpaaren vom 5′-Ende zur Entfernung von Primer- und Adaptorsequenzen mit TrimGalore. Das Sequenz-Alignment zum GRCh38/hg38-Referenzgenom und die Methylierungs-Calls wurden mit Bismark ($ bismark –multicore 8 –bowtie2 -q -N 1 INDEX INFILE.fq) durchgeführt. Coverage-Dateien (Zählungen) für Cytosine im CpG-Kontext wurden mit Bismark ($ bismark_methylation_extractor –multicore 8 –comprehensive –bedGraph INFILE_bismark_bt2.bam) erstellt.
Ganzgenom-Sequenzierungsdatenanalyse
Mutationen wurden für zwei T-ALL-Zelllinien (Jurkat und CUTLL1) und zwei T-ALL-Patientenproben aus den Ganzgenom-Sequenzierungsdaten identifiziert. Die Illumina Short-Read-Sequenzen wurden mit Hilfe von BWA mem an das menschliche Referenzgenom (GRCH38/hg38) angeglichen. Mit SAMBlaster identifizierten wir die diskordanten Paare, teilten die Reads auf und markierten die mutmaßlichen PCR-Duplikate. Mit SAMBAMBA konvertierten wir die ausgerichteten SAM-Daten in das BAM-Format, und mit samtools sortierten wir die ausgerichteten Daten, um eine BAM-Datei für jede Probe zu erstellen.
Wir verwendeten VarDict, um die Varianten zu identifizieren, die sich mit den CTCF-Bindungsstellen der Union überschnitten. Wir verwendeten alle Standardparameter mit Ausnahme von „-f 0.1“, das zur Identifizierung von Varianten verwendet wurde, die von mehr als 10% der Reads an dieser Stelle unterstützt wurden. Wir haben die Varianten mit Hilfe von Variant Effect Predictor (VEP) annotiert und benutzerdefinierte Skripte verwendet, um die Varianten zu identifizieren, die die TF-Bindung beeinflussen.
Wir haben erneut VarDict verwendet, um die Varianten in den CTCF- und NOTCH1-Genen für die vier Proben zu identifizieren. Wir verwendeten alle Standardparameter mit Ausnahme von „-f 0.1“, das zur Identifizierung von Varianten verwendet wurde, die von mehr als 10 % der Reads an dieser Stelle unterstützt wurden. Wir haben die Varianten mit Hilfe von Variant Effect Predictor (VEP) annotiert und dann gefiltert, um die Mutationen zu identifizieren, die entweder (a) nicht in mehr als 1 % der normalen menschlichen Bevölkerung vorkommen oder (b) einen CADD-Score der Deleteriosität > 20 haben oder (c) in der COSMIC-Datenbank vorhanden sind.
Integrative Modellierung und statistische Analyse
Identifizierung des CTCF-Bindungsrepertoires im menschlichen Genom
Für CTCF ChIP-seq haben wir insgesamt 793 Datensätze gesammelt, darunter 787 öffentliche Datensätze und 6 von uns selbst erstellte Datensätze (Zusatzdatei 2: Tabelle S1). Insgesamt wurden in dieser Studie 771 CTCF ChIP-seq-Datensätze mit mehr als 2000 Peaks verwendet. Jeder Datensatz kann MACS2-identifizierte CTCF-Peaks im Bereich zwischen 2050 und 198.021 liefern, mit einem Median von 46.451 und einer Gesamtzahl von 36.873.077 Peaks (Zusätzliche Datei 1: Abb. S1a). Die Verteilung der Intervalllängen zwischen benachbarten CTCF-Peaks aller 36.873.077 Peaks aus den 771 Datensätzen hat einen Wendepunkt bei ~ 150 bp (Zusatzdatei 1: Abb. S1c), der die Grenze zwischen derselben Bindungsstelle und verschiedenen Bindungsstellen anzeigt. Daher verwendeten wir 150 bps als Grenzwert für die Zusammenführung von CTCF-Peaks. In der Praxis verlängerten wir ± 75 bps von jedem Peak-Gipfel, um eine 150-bp-Region in der Mitte des Gipfels zu erzeugen, um jeden Peak zu repräsentieren, und fügten alle überlappenden Peak-Regionen zusammen, um einen Vereinigungssatz von CTCF-Bindungsstellen zu erzeugen, der 688.429 nicht überlappende Stellen enthält. Jeder Bindungsstelle wurde ein CTCF-Besetzungswert zugewiesen, definiert als die Anzahl der ChIP-seq-Datensätze, die einen Peak innerhalb der Stelle aufweisen. Dementsprechend definierten wir die Belegungshäufigkeit als das Verhältnis der Belegungsbewertung zur Gesamtzahl der CTCF-ChIP-seq-Datensätze. Um die Robustheit der identifizierten CTCF-Bindungsstellen weiter zu gewährleisten, wählten wir 285.467 Stellen mit hohem Vertrauen und einer Belegungszahl ≥ 3 für nachgeschaltete Analysen aus. CTCF-Motive innerhalb der Union Bindungsstellen wurden von FIMO mit Jaspar-Matrix (ID: MA0139.1), mit einem p-Wert Schwelle von 1e-4 gesucht. Ein Motiv mit dem kleinsten p-Wert wurde für jede CTCF-Bindungsstelle beibehalten.
Identifizierung von konstitutiven CTCF-Bindungsstellen
Die Verteilung der Belegungsscores aller 285.467 CTCF-Bindungsstellen (Zusatzdatei 1: Abb. S1d, blaue Kurve) zeigt, dass die Mehrzahl der CTCF-Bindungsstellen nur in wenigen Datensätzen vorkommt und die Anzahl der Bindungsstellen mit zunehmendem Belegungsscore abnimmt, wenn der Belegungsscore klein ist. Es gibt jedoch CTCF-Bindungsstellen, die in fast allen Datensätzen hoch konserviert sind (z. B. Bindungsstellen mit einer Belegungszahl von mehr als 600). Wir verwenden eine Potenzgesetzfunktion zur Anpassung der Verteilungskurve (blau), die in Zusatzdatei 1: Abb. S1d dargestellte Verteilungskurve (blau) an, um den Grenzwert für konstitutive CTCF-Stellen zu bestimmen. Wir bezeichnen Oi als die Anzahl der beobachteten CTCF-Bindungsstellen mit einem Besetzungsgrad gleich i und Ei als die Anzahl der erwarteten CTCF-Stellen mit einem Besetzungsgrad gleich i. Die Potenzgesetz-Anpassung an die Daten Oi kann wie folgt beschrieben werden (Zusätzliche Datei 1: Abb. S1d, grün):
Wir definieren den Cutoff A für konstitutive CTCF-Bindungsstellen als:
Mit anderen Worten sollte die Gesamtzahl der beobachteten CTCF-Stellen mit einem Besetzungswert größer als A 6-mal höher sein als erwartet. Wir haben dann A = 615 bestimmt und einen Besetzungshäufigkeits-Cutoff von 80 % verwendet, um 22.097 konstitutive CTCF-Bindungsstellen zu definieren, was der Besetzungszahl ≥ 616 in allen 771 CTCF-ChIP-seq-Datensätzen entspricht.
Identifizierung krebsspezifischer gewonnener/verlorener CTCF-Bindungsstellen
Wir haben die folgenden zwei Kriterien verwendet, um krebsspezifische verlorene CTCF-Bindungsstellen zu identifizieren: (1) Die CTCF-Bindungsstelle sollte eine geringere Belegungshäufigkeit für Datensätze dieses Krebstyps im Vergleich zur Belegungshäufigkeit für alle Datensätze aufweisen und (2) das CTCF-Bindungsniveau (quantifiziert als normalisierte ChIP-seq-Lesezahlen) an der Stelle ist in Krebsdatensätzen geringer als in anderen Datensätzen. Für die Gewinnung von CTCF-Stellen haben wir die umgekehrten Kriterien verwendet. Kurz gesagt, für jede CTCF-Bindungsstelle in jedem Krebstyp wurde der Belegungswert in den Krebsdatensätzen zusammen mit dem Belegungswert in allen 771 Datensätzen berechnet. Die CTCF-Bindungsniveaus wurden aus einer normalisierten Read-Count-Matrix gewonnen, in der die ChIP-seq-Read-Counts (RPKM) zunächst für die Vereinigung der CTCF-Bindungsstellen in allen Datensätzen berechnet wurden und dann eine Quantil-Normalisierung folgte. Wir verwendeten einen ungepaarten zweiseitigen Student’s t-Test, um den Unterschied der Bindungsniveaus zwischen verschiedenen Gruppen von Datensätzen zu quantifizieren, und der p-Wert wurde dann mit dem Benjamini-Hochberg-Verfahren angepasst. Darüber hinaus wurden die Werte für die Bindungsbelegung und die Bindungsniveaus zwischen Krebsdatensätzen und Datensätzen aus den entsprechenden normalen Geweben oder Zelltypen verglichen, um den potenziellen Störfaktor der Gewebespezifität und nicht der Krebsspezifität zu berücksichtigen. Die detaillierten Kriterien für die Identifizierung krebsspezifischer CTCF-Bindungsstellen sind im Folgenden beschrieben:
-
Krebsspezifische verlorene CTCF-Bindungsstellen: (1) Belegungshäufigkeit ≤ 0,2 in Krebsdatensätzen; (2) Belegungshäufigkeit ≥ 0,7 in 771 Datensätzen; (3) Belegungshäufigkeit ≥ 0.5 (mit occupancy score ≥ 2) in angepassten normalen Gewebedatensätzen; (4) CTCF-Levels sind bei Krebs im Vergleich zu allen anderen Datensätzen niedriger (statistic score < 0), (5) CTCF-Levels sind bei Krebs im Vergleich zu angepassten normalen Gewebedatensätzen niedriger (statistic score < 0), (6) gemittelte CTCF-Bindungssignale (RPKM) < 5 in Krebsdatensätzen.
-
Krebsspezifische gewonnene CTCF-Bindungsstellen: (1) Belegungshäufigkeit ≥ 0,5 (mit Belegungsscore ≥ 2) in Krebsdatensätzen, (2) Belegungshäufigkeit ≤ 0,2 in 771 Datensätzen, (3) Belegungsscore = 0 in angepassten Normalgewebedatensätzen, (4) CTCF-Niveaus sind bei Krebs im Vergleich zu allen anderen Datensätzen signifikant höher (FDR ≤ 0.01), (5) CTCF-Bindungsniveaus sind signifikant höher bei Krebs im Vergleich zu abgeglichenen Normalgewebedatensätzen (FDR ≤ 0,01), (6) gemittelte CTCF-Bindungssignale (RPKM) > 2 in Krebsdatensätzen.
Die spezifischen gewonnenen und verlorenen CTCF-Bindungsstellen für jeden Krebstyp sind in Zusatzdatei 4: Tabelle S3.
Quantifizierung der differentiellen Chromatin-Zugänglichkeit
Wir haben die verarbeiteten Daten aus Ref. verwendet, die eine Matrix der normalisierten ATAC-seq-Insertionszahlen innerhalb des TCGA-Pan-Krebs-Peaksatzes enthalten, um die differentielle Chromatin-Zugänglichkeit um CTCF-Bindungsstellen zu bewerten. Für jeden Krebstyp unter BRCA, CRC, LUAD und PRAD wurden die pan-Krebs-ATAC-seq-Peaks, die sich mit identifizierten krebstypspezifischen verlorenen oder gewonnenen CTCF-Bindungsstellen überschneiden, für nachgeschaltete Analysen verwendet. Der ATAC-seq-Differential-Score für jeden Peak wurde als fache Änderung des Durchschnitts der normalisierten ATAC-seq-Insertionszahlen von Patientenproben des entsprechenden Krebstyps im Vergleich zu Patienten anderer Krebstypen quantifiziert, und der ATAC-seq-Differential-Score wurde dann der Peak-überlappenden CTCF-Bindungsstelle zugewiesen.
Aus Gründen der Konsistenz haben wir den gleichen Ansatz wie bei den TCGA-ATAC-seq-Daten angewandt, um die gesammelten ATAC-seq-Daten der T-ALL-Zelllinie Jurkat und der normalen CD4+ T-Zellen zu analysieren. Aus den ATAC-seq-Rohdaten der CTCF-Bindungsstellen der Union wurde eine Datenmatrix für alle Jurkat- und T-Zell-Datensätze erstellt. Auf die log2-skalierte Matrix wurde eine Quantil-Normalisierung angewandt (Pseudoanzahl = 5). Der ATAC-seq-Differential-Score wurde als fache Änderung der gemittelten normalisierten ATAC-seq-Zahlen zwischen den Datensätzen von Jurkat und CD4+ T-Zellen an jeder CTCF-Bindungsstelle gemessen.
Normalisierung von Chromatin-Interaktionen
Bei einer Hi-C-Kontaktkarte A = {aij} spiegelt der Score aij die gemappten Reads zwischen zwei genomischen Regionen i und j wider. Angenommen, die Bin-Größe beträgt 5 kb, dann haben die Regionen i und j einen genomischen Abstand von ∣i – j ∣ × 5kb. Da die Kontaktwahrscheinlichkeit zwischen zwei Bins mit zunehmendem genomischen Abstand abnimmt, haben wir die Kontaktkarte wie folgt normalisiert: Für jeden gegebenen genomischen Abstand dk = k × 5kb quantifizieren wir einen Normalisierungsfaktor \( {\overline{S}}_{d_k} \) als die gemittelten Interaktionen zwischen allen Bin-Paaren mit demselben genomischen Abstand dk im selben Chromosom, z. B., \( {\overline{S}}_{d_k}=\left({\sum}_{j-i=k}{a}_{ij}\right)/n \), wobei n die Gesamtzahl der Bin-Paare mit dem Abstand dk ist. Der Interaktionswert aij zwischen zwei Bins mit dem Abstand dk wurde dann durch \( {\overline{S}}_{d_k} \) als \( {a}_{ij}^{\prime }={a}_{ij}/{\overline{S}}_{d_k} \) normalisiert. Mit diesem Ansatz normalisierten wir die Matrix A in \( A^{\prime }=\left\{{a}_{ij}^{\prime}\right\} \) innerhalb jedes Chromosoms.
Detektion von differentiellen Chromatin-Interaktionen
Wir bezeichneten die normalisierten Hi-C-Kontaktkarten im Krebsdatensatz und im normalen Datensatz als C = {cij} bzw. N = {nij}. Für eine gegebene CTCF-Bindungsstelle x (mit der Koordinate xc) und einem vordefinierten genomischen Abstand L werden die Chromatin-Interaktionen zwischen x und den benachbarten, nicht überlappenden 5-kb-Bins mit einem genomischen Abstand von bis zu L in C bzw. N gesammelt. Konkret werden die Interaktionswerte zwischen x und seinen benachbarten 5-kb-Bins in C als IC = {cij} , wobei entweder i oder j gleich ⌊xc/5kb⌋ ist und 0 < (j – i) × 5kb ≤ L. In ähnlicher Weise wurden die Interaktionswerte zwischen x und den benachbarten 5-kb-Bins in N als IN = {nij} erfasst. Ein gepaarter zweiseitiger Student’s t-Test wurde dann auf IC und IN angewendet, um die unterschiedliche Interaktion zwischen Krebs- und normalen Zellen in der Umgebung der CTCF-Bindungsstelle x zu quantifizieren.
Assoziation der CTCF-Bindung mit der Genexpression
Insgesamt wurden 54 Zelltypen, für die sowohl CTCF-ChIP-seq-Daten als auch RNA-seq-Daten öffentlich verfügbar sind, ausgewählt (Additional file 6: Table S5), um die Assoziation zwischen CTCF-Bindung und Genexpression für jedes CTCF-Gen-Paar auf demselben Chromosom zu untersuchen. Um das CTCF-Bindungsniveau zu erhalten, wurde eine Read-Count-Matrix unter Verwendung von Reads per Kilobase per Million (RPKM) auf Union-CTCF-Bindungsstellen aus ChIP-seq-Daten erstellt. Die Read-Count-Matrix wurde mit der Quadratwurzel aus RPKM skaliert, gefolgt von einer Quantilsnormalisierung. Das Niveau der Genexpression wurde für jedes Gen anhand der Quadratwurzel der Transkripte pro Million (TPM) aus RNA-seq-Daten gemessen. Für jedes CTCF-Gen-Paar quantifizierten wir die Assoziation zwischen der CTCF-Stelle und dem Gen in allen 54 Zelltypen anhand des Korrelationskoeffizienten R zwischen der normalisierten CTCF-Bindungsstärke und der Genexpression (Abb. 3a). CTCF-Gene-Paare wurden als „hoch korreliert“ eingestuft, wenn R2 größer als 0,25 war, z. B, Korrelationskoeffizient größer als 0,5 oder kleiner als – 0,5, und die hoch korrelierten CTCF-Gen-Paare tragen zu 1,3 % aller CTCF-Gen-Paare bei (Additional file 1: Abb. S8a).
Identifizierung konstitutiver CTCF-gebundener Chromatindomänen
Für jede CTCF-Bindungsstelle definierten wir die zugehörige Chromatindomäne als die genomische Region, die (1) diese spezifische CTCF-Bindungsstelle einschließt, (2) von einem Paar konstitutiver CTCF-Bindungsstellen mit Motiven entgegengesetzter Orientierung begrenzt wird und (3) mindestens 100 kb und maximal 1 MB Region auf jeder Seite der CTCF-Bindungsstelle einnimmt. Abbildung 3b zeigt schematisch, wie konstitutive CTCF-begrenzte Chromatindomänen definiert wurden.
Detektion von DNA-Methylierungsänderungen in der Umgebung von CTCF-Bindungsstellen
DNA-Methylierungsänderungen wurden innerhalb einer 300-bp-Region in der Mitte jeder CTCF-Bindungsstelle festgestellt. Regionen mit mindestens 3 CpGs, die von mindestens 5 Reads (≥ 5×) sowohl in Krebszelllinien als auch in den entsprechenden Normalgeweben abgedeckt wurden, wurden ausgewählt. Eine 300-bp-Region wurde als differentiell methyliert erkannt, wenn die durchschnittlichen differentiellen Methylierungsniveaus aller CpGs (≥ 5×) innerhalb dieser Region mehr als 20 % betrugen.
Ermittlung der Mutationsrate und des differentiellen Motiv-Scores
Für jede CTCF-Bindungsstelle wurde die rohe Mutationszahl als das Auftreten von Mutationsereignissen in allen Proben/Patienten an jedem einzelnen Basenpaar innerhalb einer 400-bp-Region berechnet, die an der CTCF-Bindungsstelle zentriert war. Die Mutationsrate für eine Gruppe von CTCF-Bindungsstellen wurde als gemittelte Mutationszahl über die Anzahl der CTCF-Bindungsstellen für jedes Basenpaar innerhalb der 400-bp-Region berechnet.
Der Motiv-Score wurde gemessen, indem die CTCF-Positionsgewichtsmatrix (Jaspar , Matrix ID: MA0139.1) mit einer 19-bp-DNA-Sequenz, die auf das CTCF-Motiv oder die CTCF-Bindungsstelle zentriert ist, unter Verwendung von Log-Likelihood-Verhältnissen (mit einer Hintergrundnukleotidhäufigkeit wie für A,C,G,T) verglichen wurde. Der differentielle Motiv-Score wurde durch Vergleich der Motiv-Scores für die Referenz- und die mutierten Sequenzen berechnet.
DNA-Sequenz-Motiv-Analyse
Die Analyse der Anreicherung von DNA-Sequenz-Motiven wurde mit MDSeqPos (Version 1.0.0) auf Cistrome mit Standardparametern (-cisrome -Homo Sapien oder Mus musculus) durchgeführt. De-novo-Motivanalysen wurden mit HOMER (Version 4.10) mit findmotifs.pl-Modul und MEME (Version 5.1.1) mit den folgenden Parametern: meme -dna -mod zoops -maxw 20 -evt -0.01.
Identifizierung von CTCF-Intra-Domain-Regionen, die unterschiedlich interagieren
Für einen gegebenen Satz von CTCF-Bindungsstellen wurden die Chromatin-Interaktionsänderungen zwischen einer CTCF-Stelle und jedem ihrer nicht überlappenden Intra-Domain-Bins, gemessen aus normalisierten Hi-C-Kontaktkarten in Krebszellen im Vergleich zu angepassten normalen Zellen, für jede der CTCF-Bindungsstellen gesammelt (Additional file 1: Abb. S14b). Regionen mit verminderten Interaktionen (log2 FC < -1, gemittelt log2 Interaktion > 0) mit krebsspezifischen verlorenen CTCF-Bindungsstellen und Regionen mit erhöhten Interaktionen (log2 FC > 1, gemittelt log2 Interaktion > 0) mit krebsspezifischen gewonnenen CTCF-Bindungsstellen wurden für die nachgeschaltete Transkriptionsfaktor (TF) Anreicherungsanalyse verwendet.
Transkriptionsfaktor-Anreicherungsanalyse
Eine überarbeitete Version des BART-Algorithmus wurde für die TF-Anreicherungsanalyse verwendet. Kurz gesagt, wurde eine Sammlung von union DNase I hypersensitive sites (UDHS) als Repertoire aller cis-regulierenden Elemente im menschlichen Genom kuratiert, und 7032 ChIP-seq-Datensätze wurden für 883 TFs gesammelt, wobei jeder TF einen oder mehrere ChIP-seq-Datensätze aus mehreren Zelltypen oder Bedingungen hatte. Für jeden TF auf UDHS wurde ein binäres Profil erstellt, das angibt, ob der TF mindestens einen Peak aus einem seiner ChIP-seq-Datensätze in jedem der UDHS aufweist. Die Analyse der Bindungsanreicherung wurde für jeden TF durchgeführt, indem die TF-Bindung auf einer Untergruppe von UDHS, die sich mit den ausgewählten genomischen Regionen überschneidet, mit der TF-Bindung auf UDHS verglichen wurde. Der p-Wert wurde mit dem exakten Test von Fisher mit zwei Schwellenwerten ermittelt.