Procédure expérimentale
Xénogreffe de patients et culture cellulaire
Les lignées cellulaires T-ALL humaines comprennent CUTLL1 (don d’Adolfo Ferrando, Columbia University) et JURKAT (American Type Culture Collection (ATCC), Manassas, VA, #CCL-119) . Les cellules ont été cultivées dans un milieu RPMI1640 avec de la l-glutamine et 25 mM HEPES (Corning) complété par 10 % de sérum fœtal bovin inactivé par la chaleur (Sigma-Aldrich), 10 U/mL de pénicilline-streptomycine (Gibco) et 1× glutaMAX (Gibco) dans un incubateur humidifié à 37 °C et 5 % de CO2. Les cellules sont périodiquement testées pour la présence de mycoplasmes à l’aide du kit de détection de mycoplasmes MycoAlert de Lonza Walkersville (dernier test en janvier 2020). Les lignées cellulaires sont maintenues en culture pour un maximum de 20 passages et sont authentifiées à l’aide du profilage des répétitions en tandem courtes (JURKAT) ou à l’aide de la PCR pour détecter la translocation TCRb-NOTCH1 (TCRBJ2S4CUTLL1F:5′-GGACCCGGCTCTCAGTGCT-3′, NOTCH1CUTTL1R:5′-TCCCGCCCTCCAAAATAAGG-3′). La dernière authentification des cellules a été effectuée en février 2020. Les cellules T CD4+ humaines ont été achetées auprès d’AllCells. Les échantillons humains primaires ont été collectés avec un consentement éclairé et analysés sous la supervision de l’Institutional Review Board de l’Université de Padoue, de l’Associazone Italiana di Ematologia e Oncologia Pediatrica et des essais cliniques pédiatriques ALL 2000/2006 de Berlin-Frankfurt-Münster (AIEOP-BFM). Le consentement éclairé pour utiliser le matériel restant à des fins de recherche a été obtenu de tous les patients à l’entrée de l’essai, conformément à la Déclaration d’Helsinki.
Anti-corps et réactifs
Des Western blots ont été réalisés en utilisant les anticorps suivants : Actin et CTCF de Millipore Sigma (clone C4 ; 07-729) et NOTCH1 clivé (Val1744) de Cell Signaling Technology (4147). Les ChIP-seq ont été réalisés en utilisant les anticorps suivants : CTCF de Millipore Sigma (07-729), H3K27Ac (8173S) et H3K27me3 (9733S) de Cell Signaling Technology, et H3K4me1 (07-473) de Millipore.
Hi-C in situ
La Hi-C in situ a été réalisée sur des cellules T CD4+, Jurkat, CUTLL1 et des xénogreffes de patients comme décrit précédemment . En bref, les cellules ont été réticulées avec du formaldéhyde à 1% pendant 10 minutes à température ambiante. Par réaction Hi-C, 5 millions de cellules ont été lysées et les noyaux ont été perméabilisés. L’ADN a été digéré avec MboI de New England Biolabs (R0147M). Les fragments digérés ont été marqués avec du d-ATP biotinylé de Jena Bioscience (NU-835-BIO14-S) et ligaturés. Après un traitement à la RNase et à la Protéinase K pour inverser les liaisons croisées, les noyaux ont été soniqués à l’aide d’un Covaris E220 pour produire une longueur moyenne de fragment de 400 pb. Les billes de streptavidine de Thermo Fisher Scientific (65001) ont été utilisées pour attirer les fragments marqués à la biotine. Après la purification et l’isolement de l’ADN, les bibliothèques finales ont été préparées à l’aide du kit de préparation de bibliothèque d’ADN NEBNext® Ultra™ II pour Illumina® et séquencées par séquençage à extrémité appariée à une longueur de lecture de 150 pb sur un Illumina HiSeq 2500 pour produire en moyenne 400 millions de lectures par échantillon.
Profilage CHIP-seq
Les cellules TCD4+, Jurkat, CUTLL1 et les xénogreffes de patients ont été réticulées avec 1% de formaldéhyde et 1% de sérum bovin fœtal dans du PBS pendant 10 min à température ambiante. La réaction a été éteinte avec de la glycine 0,2 M à température ambiante pendant 5 minutes. Les cellules ont ensuite été lavées avec du PBS et pelletées.
Pour les ChIPs de CTCF, l’immunoprécipitation a été réalisée sur la base d’un protocole décrit précédemment . Un culot contenant 50 millions de cellules a été lysé avec 5 mL de tampon de lyse (50 mM HEPES-KOH, pH 7,5, 140 mM NaCl, 1 mM EDTA, 10% de glycérol, 0,5% NP-40, 0,25% Triton X-100) pendant 10 min à 4 °C. Les noyaux ont été granulés à 1350×g pendant 7 minutes et remis en suspension dans 10 mM de Tris pH 8, 1 mM d’EDTA et 0,1% de SDS. La chromatine a été cisaillée avec un système Covaris E220 jusqu’à une longueur de fragment moyenne de 400 pb et centrifugée à 15 000 tr/min pendant 10 minutes pour éliminer la chromatine insoluble et les débris. Le surnageant a été incubé avec 20 μL de Dynabeads Protein G pendant 30 min avant de jeter les billes. Un pour cent du volume total a été conservé comme intrant et le reste a été incubé avec l’anticorps anti-CTCF pendant la nuit. Au total, 100 μL de Dynabeads Protein G ont été ajoutés pendant 2 h. Les fragments liés ont été lavés deux fois avec 1 mL de tampon à faible teneur en sel (20 mM Tris-HCl pH 8.0, 150 mM NaCl, 2 mM EDTA, 1% p/v Triton X-100, et 0.1 % p/v de SDS), une fois avec un tampon très salé (20 mM Tris-HCl pH 8.0, 500 mM NaCl, 2 mM EDTA, 1 % p/v de Triton X-100 et 0,1 % p/v de SDS), une fois avec un tampon au chlorure de lithium (10 mM Tris-HCl pH 8.0, 250 mM LiCl, 1 mM EDTA, 1 % p/v NP-40, et 1 % p/v acide désoxycholique), et deux fois avec TE (10 mM Tris pH 8, 1 mM EDTA).
Pour les ChIPs d’histones, les cellules ont été lysées dans 375 μL de tampon d’incubation des noyaux (15 mM Tris pH 7.5, 60 mM KCl, 150 mM NaCl, 15 mM MgCl2, 1 mM CaCl2, 250 mM sucrose, 0,3% NP-40, 1 mM NaV, 1 mM NaF et 1 comprimé d’inhibiteur de protéase sans EDTA (Roche)/10 mL dans H2O) pendant 10 min sur glace. Les noyaux ont été lavés une fois avec un tampon de digestion (10 mM NaCl, 10 mM Tris pH 7,5, 3 mM MgCl2, 1 mM CaCl2, 1 mM NaV, 1 mM NaF, et 1 comprimé d’inhibiteur de protéase sans EDTA (Roche)/10 mL dans H2O) et remis en suspension dans 57-μL de tampon de digestion contenant 4,5 unités de MNase (USB) pendant 1 h à 37 °C. L’activité MNase a été éteinte pendant 10 min sur la glace lors de l’ajout d’EDTA à une concentration finale de 20 mM. Les noyaux ont été granulés et remis en suspension dans 300-μL de tampon de lyse des noyaux (50 mM Tris-HCl pH 8,0, 10 mM EDTA pH 8,0, 1 % SDS, 1 mM NaV, 1 mM NaF et 1 comprimé d’inhibiteur de protéase sans EDTA (Roche)/10 mL dans H2O) avant sonication avec un Bioruptor Pico (Diagenode) pendant 5 min (30 s on, 30 s off). Le lysat a été centrifugé à la vitesse maximale pendant 5 min pour éliminer les débris. Neuf volumes de tampon de dilution IP (0,01 % de SDS, 1,1 % de Triton X-100, 1,2 mM d’EDTA pH 8,0, 16,7 mM de Tris-HCl pH 8,0, 167 mM de NaCl, 1 mM de NaV, 1 mM de NaF et 1 comprimé d’inhibiteur de protéase sans EDTA (Roche)/10 mL dans H2O) ont été ajoutés au surnageant. Au total, 50 μL de Dynabeads Protein G ont été ajoutés et l’échantillon a été incubé à 4 °C pendant 30 min, en rotation. Un pour cent de l’échantillon a été conservé en entrée, et l’échantillon restant a été divisé en 3 tubes. Au total, 50 μL de Dynabeads Protein G conjugués à 15 μL de l’anticorps approprié ont été ajoutés à chaque tube avant une incubation d’une nuit à 4 °C, en rotation. Les complexes liés aux billes ont été lavés pendant 5 min chacun dans 1 mL de tampon faiblement salé, de tampon fortement salé, de tampon LiCl et deux fois avec TE.
Pour éluer les complexes liés aux billes, les billes ont été remises en suspension dans 50 μL de tampon d’élution (100 mM NaHCO3, 1 % p/v de SDS) et incubées à 65 °C pendant 15 min, en les agitant à 1000 tours/minute sur un thermomixeur (Thermo Scientific). L’élution a été répétée une seconde fois, puis 100 μL de tampon RNase (12 μL de NaCl 5 M, 0,2 μL de RNase 30 mg/mL et 88 μL de TE) ont été ajoutés à chaque échantillon de ChIP et d’entrée. Les échantillons ont été incubés à 37 °C pendant 20 min, puis 100 μL de tampon de protéinase K (2,5 μL 20 mg/mL de protéinase K, 5 μL 20 % de SDS et 92,5 μL de TE) ont été ajoutés pendant la nuit à 65 °C. Un volume égal de solution de phénol:chloroforme a été ajouté et mélangé soigneusement. Le mélange a été transféré dans des tubes MaXtract haute densité (Qiagen) et centrifugé pendant 8 min à 15 000 rpm. La phase supérieure a été transférée dans de nouveaux tubes et mélangée avec 1,5 μL de glycogène à 20 mg/mL, 30 μL d’acétate de sodium 3M et 800 μL d’éthanol. Les échantillons ont été incubés à – 80 °C jusqu’à leur congélation, puis centrifugés à 15 000 tours/minute pendant 30 min à 4 °C. Le surnageant a été retiré et les culots ont été lavés dans 800 μL d’éthanol glacé à 70 % et centrifugés pendant 10 min à 4 °C à 15 000 tr/min. Après avoir soigneusement éliminé l’éthanol, les culots ont été séchés à l’air et remis en suspension dans 30 μL de Tris 10 mM à pH 8.
IP et l’ADN d’entrée ont ensuite été quantifiés à l’aide d’un fluoromètre Qubit 3.0. Les librairies ont été préparées à l’aide du kit KAPA HyperPrep (KK8505) et séquencées avec un Illumina NextSeq 500 à une profondeur moyenne de 28 millions de lectures par échantillon.
Profilage ARN-seq
L’ARN a été isolé à partir de 3 millions de cellules par échantillon à l’aide du mini kit d’ARN total Bio-Rad Aurum™ et quantifié avec le kit Agilent RNA 6000 Nano avec le bioanalyseur Agilent. Les librairies ont été préparées par déplétion de l’ARNr à l’aide de la trousse de préparation de librairies d’ARNm strandé TruSeq® d’Illumina pour une faible concentration d’échantillon de départ et séquencées par séquençage à une seule extrémité sur un NextSeq 500 d’Illumina à une profondeur moyenne de 18 millions de lectures par échantillon.
Profilage de la méthylation de l’ADN
L’ADN génomique a été isolé à l’aide de la trousse AllPrep DNA/RNA Micro Kit (Qiagen). Pour évaluer l’état de méthylation de l’ADN à l’échelle du génome, nous avons effectué un mRRBS . Après une quantification fluorométrique à l’aide d’un instrument Qubit 3.0, nous avons digéré l’ADN génomique avec l’enzyme de restriction MspI (New England Biolabs) et sélectionné la taille des fragments d’environ 100 à 250 paires de bases en utilisant des billes d’immobilisation réversible en phase solide (SPRI) (MagBio Genomics). L’ADN obtenu a subi une conversion au bisulfite à l’aide du kit EZ DNA Methylation-Lightning (Zymo Research). Nous avons créé des bibliothèques à partir d’ADN monocaténaire converti au bisulfite à l’aide du Pico Methyl-Seq Library Prep Kit (Zymo Research), qui ont ensuite été regroupées pour le séquençage sur un instrument NextSeq 500 d’Illumina à l’aide du kit de réactifs NextSeq 500/550 V2 High Output (1 × 75 cycles) jusqu’à une profondeur de lecture minimale de 50 millions de lectures par échantillon.
Séquençage du génome entier
Trois millions de cellules provenant de lignées cellulaires ou d’échantillons de patients ont été culottées et remises en suspension dans 1 mL de solution de lyse cellulaire (Qiagen) mélangée à 500 μg de RNase A. La réaction de lyse a été réalisée à 37 °C pendant 15 min. Au total, 333 μL de solution de précipitation des protéines (Qiagen) ont été ajoutés à chaque échantillon qui a ensuite été mélangé au vortex puis centrifugé à 2000×g pendant 10 min. Le surnageant a été mélangé avec 1 mL d’isopropanol jusqu’à ce que les brins d’ADN précipitent de la solution. Après avoir jeté le surnageant, le culot d’ADN a été lavé avec 1 ml d’éthanol à 70 % et centrifugé à 2000×g pendant 1 minute. L’éthanol a ensuite été versé et le culot a été séché à l’air pendant 15 min avant d’être remis en suspension dans 50 à 100 μL de solution d’hydratation de l’ADN (Qiagen). L’ADN a été séquencé par séquençage Illumina en paires avec une couverture de 30×.
Immunoprécipitation
Un total de 100 millions de cellules pour chaque réaction d’immunoprécipitation a été culotté et incubé dans le tampon A (10 mM HEPES pH 8,0, 1,5 mM MgCl2, 10 mM KCl, 0,5 mM DTT) pendant 10 min sur glace. Les cellules ont ensuite été lysées en 12 coups avec un broyeur de tissu à pilon libre de 7 ml (Wheaton, 357542) et centrifugées à 2000 rpm pendant 7 min. Les culots nucléaires ont été remis en suspension dans 5 volumes de tampon TENT (50 mM Tris pH 7,5, 5 mM EDTA, 150 mM NaCl, 1% Triton X-100, 5 mM MgCl2) et traités à la benzonase pendant 30 min avant 5 passages dans une seringue de 25 g × 5/8 in. La fraction insoluble a été éliminée après centrifugation à 2000 rpm pendant 7 min et incubée pendant la nuit avec des Dynabeads Protein G hybridés avec un anticorps. Un total de 2 millions de cellules a été prélevé pour l’entrée. Les lysats de billes et de noyaux ont été lavés 6 fois avec le tampon TENT, puis élués dans de la glycine 0,1 M pH 2,5 avec 100 mM Tris pH 8,0 avant. Le tampon d’échantillon NuPAGE LDS a été ajouté aux éluats et aux entrées, qui ont ensuite été incubés à 70 °C pendant 15 min avant d’être analysés par western blot.
Collecte de données publiques
Les données CTCF ChIP-seq publiques ont été collectées à partir de Cistrome Data Browser (pour les fichiers de pics) et NCBI GEO (pour les fichiers fastq, fichier supplémentaire 2 : tableau S1). Les données ChIP-seq de modification d’histone ont été collectées auprès de NCBI GEO et ENCODE (pour les fichiers bam). Les données publiques RNA-seq dans plusieurs types de cellules ont été collectées auprès d’ENCODE (pour les fichiers fastq). Les données de profilage de méthylation de l’ADN ont été recueillies auprès d’ENCODE (pour les fichiers bedMethyl) et de NCBI GEO. Les données Hi-C ont été collectées auprès de NCBI GEO et ENCODE (pour les fichiers fastq). Les données ATAC-seq ont été collectées auprès de NCBI GEO (pour les fichiers fastq). Les données de séquençage du génome entier pour les échantillons BRCA, COAD, LUAD et PRAD ont été recueillies sur le portail de données du Consortium international du génome du cancer (ICGC). Des informations détaillées, y compris les identifiants d’accession de tous les ensembles de données publics recueillis dans ce travail, peuvent être trouvées dans le fichier supplémentaire 6 : Tableau S5.
Traitement des données
Analyse des données ChIP-seq
L’alignement de séquence pour les données ChIP-seq dans les fichiers fastq a été effectué en utilisant le même pipeline d’analyse standard que celui utilisé dans Cistrome DB , pour la cohérence et la reproductibilité. Tous les alignements génomiques de données de séquence ont été effectués à l’aide du pipeline Chilin avec les paramètres par défaut ($ chilin simple -p narrow -s hg38 –threads 8 -t IN.fq -i PRENAME -o OUTDIR). En bref, les lectures de séquence ont été alignées sur le génome de référence humain (GRCH38/hg38) à l’aide de 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). Les fichiers Sam ont ensuite été convertis en fichiers bam à l’aide de samtools ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Pour les ensembles de données CTCF ChIP-seq, MACS2 a été utilisé pour appeler les pics sous le seuil FDR de 0,01 ($ macs2 callpeak –SPMR -B -q 0,01 –keep-dup 1 -g hs -t PRENAME.bam -n PRENAME –outidr OUTDIR). Les pics dont l’enrichissement en plis était d’au moins 4 ont été retenus. Les fichiers Bigwiggle ont été générés à l’aide de BEDTools et des outils 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). Enfin, seuls les échantillons CTCF ChIP-seq qui ont au moins 2000 pics ont été inclus dans l’analyse intégrative en aval.
Analyse des données ATAC-seq
Trim Galore a été utilisé pour tailler les lectures de séquençage brutes ($ trim_galore –nextera –phred33 –fastqc –paired R1.fq R2.fq -o OUTDIR). Les lectures ont été alignées sur le génome de référence humain (GRCH38/hg38) en utilisant Bowtie2 ($ bowtie2 -p 10 -X 2000 -x INDEX -1 R1.fq -2 R2.fq -S PRENAME.sam). Les fichiers Sam ont ensuite été convertis en fichiers bam à l’aide de samtools ($ samtools view -bS -q 1 -@ 8 PRENAME.sam > PRENAME.bam). Bedtools a été utilisé pour convertir les fichiers bam en format bed ($ bamToBed -i PRENAME.bam -bedpe > PRENAME_PE.bed). Les lectures mappées à l’ADN mitochondrial ont été écartées de l’analyse en aval.
Analyse des données ARN-seq
Les ensembles de données ARN-seq ont été traités à l’aide de Salmon ($ salmon quant –gcBias -i INDEX -l A -p 8 {-1 R1.fq -2 R2.fq| -r IN.fq} -o OUTDIR). L’index du transcriptome a été construit sur le génome de référence humain (GRCH38/hg38). Les estimations d’abondance au niveau des transcriptions ont été résumées au niveau des gènes en utilisant le paquet « tximport » pour l’analyse d’expression différentielle. DESeq2 a été utilisé pour identifier les gènes exprimés de manière différentielle, et les différents seuils utilisés dans les différentes analyses ont été listés de manière correspondante dans le manuscrit.
Analyse des données Hi-C
Les données Hi-C ont été traitées à l’aide de HiC-Pro ($ HiC-Pro -i INDIR -o OUTDIR -c CONFIG -p). Les cartes de contact ont été générées à une résolution de 5 kb. Les données matricielles brutes ont été normalisées à l’aide de l’approche décrite dans Normalization of Chromatin Interactions.
Analyse des données de méthylation de l’ADN
Les données de méthylation de l’ADN (pour les lignées cellulaires T-ALL et les patients T-ALL) ont été démultiplexées avec bcl2fastq, suivies d’un rognage de 10 paires de bases de l’extrémité 5′ pour éliminer les séquences d’amorce et d’adaptateur à l’aide de TrimGalore . L’alignement des séquences sur le génome de référence GRCh38/hg38 et les appels de méthylation ont été effectués avec Bismark ($ bismark –multicore 8 –bowtie2 -q -N 1 INDEX INFILE.fq). Les fichiers de couverture (comptes) pour les cytosines en contexte CpG ont été générés à l’aide de Bismark ($ bismark_methylation_extractor –multicore 8 –comprehensive –bedGraph INFILE_bismark_bt2.bam).
Analyse des données de séquençage du génome entier
Des mutations ont été identifiées pour deux lignées cellulaires T-ALL (Jurkat et CUTLL1) et deux échantillons de patients T-ALL à partir des données de séquençage du génome entier. Nous avons aligné les séquences à lecture courte d’Illumina sur le génome de référence humain (GRCH38/hg38) en utilisant BWA mem. Nous avons utilisé SAMBlaster pour identifier les paires discordantes, diviser les lectures et signaler les duplicatas PCR putatifs. Nous avons utilisé SAMBAMBA pour convertir les alignements SAM au format BAM, et samtools a été utilisé pour trier ces alignements afin de créer un fichier BAM correspondant à chaque échantillon.
Nous avons utilisé VarDict pour identifier les variants qui chevauchent les sites de liaison CTCF de l’union. Nous avons utilisé tous les paramètres par défaut à l’exception de « -f 0,1 » qui a été utilisé pour identifier les variants qui étaient soutenus par plus de 10% des lectures à cet endroit. Nous avons annoté les variants à l’aide du Variant Effect Predictor (VEP) et utilisé des scripts personnalisés pour identifier les variants qui influencent la liaison TF.
Nous avons à nouveau utilisé VarDict pour identifier les variants dans les gènes CTCF et NOTCH1 pour les quatre échantillons. Nous avons utilisé tous les paramètres par défaut à l’exception de « -f 0,1 » qui a été utilisé pour identifier les variants qui étaient soutenus par plus de 10% des lectures à cet endroit. Nous avons annoté les variants à l’aide de Variant Effect Predictor (VEP) , puis nous les avons filtrés pour identifier les mutations qui (a) n’étaient pas observées dans plus de 1 % d’une population humaine normale, ou (b) avaient un score CADD de délétérité > 20, ou (c) étaient présentes dans la base de données COSMIC.
Modélisation intégrative et analyse statistique
Identification du répertoire de liaison CTCF dans le génome humain
Pour CTCF ChIP-seq, nous avons collecté un total de 793 jeux de données, dont 787 jeux de données publics et 6 jeux de données que nous avons générés (fichier supplémentaire 2 : tableau S1). Au total, 771 jeux de données CTCF ChIP-seq avec des pics supérieurs à 2000 ont été utilisés dans cette étude. Chaque ensemble de données peut produire des pics CTCF identifiés par MACS2 dans une fourchette comprise entre 2050 et 198 021, avec une médiane de 46 451 et un total de 36 873 077 pics (Additional file 1 : Fig. S1a). La distribution des longueurs d’intervalle entre les sommets adjacents des pics CTCF des 36 873 077 pics des 771 ensembles de données présente un point d’inflexion à ~ 150 pb (fichier supplémentaire 1 : figure S1c), ce qui indique la limite entre le même site de liaison et des sites de liaison différents. Par conséquent, nous avons utilisé 150 pb comme seuil pour fusionner les pics CTCF. En pratique, nous avons étendu ± 75 pb de chaque sommet de pic pour générer une région de 150 pb centrée sur le sommet pour représenter chaque pic et avons fusionné toutes les régions de pic qui se chevauchent pour générer un ensemble d’union de sites de liaison CTCF, qui contient 688 429 sites non chevauchants. Chaque site de liaison s’est vu attribuer un score d’occupation de CTCF, défini comme le décompte des ensembles de données ChIP-seq qui présentent un pic dans le site. En conséquence, nous avons défini la fréquence d’occupation comme le rapport entre le score d’occupation et le nombre total d’ensembles de données ChIP-seq CTCF. Pour garantir la robustesse des sites de liaison CTCF identifiés, nous avons sélectionné 285 467 sites de haute confiance avec un score d’occupation ≥ 3 pour les analyses en aval. Les motifs CTCF au sein des sites de liaison syndicaux ont été recherchés par FIMO avec la matrice Jaspar (ID : MA0139.1), avec un seuil de valeur p de 1e-4. Un motif avec la plus petite valeur p a été retenu pour chaque site de liaison CTCF.
Identification des sites de liaison CTCF constitutifs
La distribution des scores d’occupation de l’ensemble des 285 467 sites de liaison CTCF (fichier additionnel 1 : Fig. S1d, courbe bleue) montre que la majorité des sites de liaison CTCF ne se trouvent que dans quelques ensembles de données, et que le nombre de sites de liaison diminue avec l’augmentation du score d’occupation lorsque celui-ci est faible. Cependant, certains sites de liaison CTCF sont très conservés dans presque tous les ensembles de données (par exemple, les sites de liaison avec un score d’occupation supérieur à 600). Nous utilisons une fonction de loi de puissance pour ajuster la courbe de distribution (bleue) présentée dans le fichier supplémentaire 1 : Fig. S1d pour déterminer le seuil des sites CTCF constitutifs. Nous désignons Oi comme le nombre de sites de liaison CTCF observés avec un score d’occupation égal à i, et Ei comme le nombre de sites CTCF attendus avec un score d’occupation égal à i. L’ajustement de la loi de puissance aux données Oi peut être décrit comme (fichier supplémentaire 1 : Fig. S1d, vert):
Nous définissons le seuil A pour les sites de liaison CTCF constitutifs comme:
En d’autres termes, le total des sites CTCF observés avec un score d’occupation supérieur à A devrait être 6 fois plus important que prévu. Nous avons alors déterminé A = 615, et utilisé un seuil de fréquence d’occupation de 80 % pour définir 22 097 sites de liaison CTCF constitutifs, ce qui correspond au score d’occupation ≥ 616 dans les 771 ensembles de données ChIP-seq CTCF.
Identification des sites de liaison CTCF gagnés/perdus spécifiques au cancer
Nous avons utilisé les 2 critères suivants pour identifier les sites de liaison CTCF perdus spécifiques au cancer : (1) Le site de liaison CTCF doit avoir une fréquence d’occupation plus faible pour les ensembles de données de ce type de cancer par rapport à la fréquence d’occupation de tous les ensembles de données et (2) le niveau de liaison CTCF (quantifié en tant que nombre de lectures ChIP-seq normalisées) sur le site est plus faible dans les ensembles de données du cancer que dans les autres ensembles de données. Pour les sites CTCF gagnés, nous avons utilisé l’ensemble de critères inverses. En bref, pour chaque site de liaison CTCF dans chaque type de cancer, le score d’occupation dans les ensembles de données sur le cancer a été calculé ainsi que son score d’occupation dans les 771 ensembles de données. Les niveaux de liaison CTCF ont été obtenus à partir d’une matrice de comptage de lecture normalisée dans laquelle les comptages de lecture ChIP-seq (RPKM) ont d’abord été calculés pour les sites de liaison CTCF de l’union dans tous les ensembles de données, puis suivis d’une normalisation par quantile. Nous avons utilisé le test t de Student non apparié et à double sens pour quantifier la différence des niveaux de liaison entre les différents groupes d’ensembles de données, et la valeur p a ensuite été ajustée en utilisant la procédure de Benjamini-Hochberg. En outre, les scores d’occupation des liaisons et les niveaux de liaison ont été comparés entre les ensembles de données sur le cancer et les ensembles de données sur les types de tissus ou de cellules normaux appariés, afin de prendre en compte le facteur de confusion potentiel de la spécificité du tissu plutôt que du cancer. Les critères détaillés d’identification des sites de liaison CTCF spécifiques du cancer sont décrits ci-dessous :
-
Sites de liaison CTCF perdus spécifiques du cancer : (1) fréquence d’occupation ≤ 0,2 dans les ensembles de données sur le cancer ; (2) fréquence d’occupation ≥ 0,7 dans 771 ensembles de données ; (3) fréquence d’occupation ≥ 0.5 (avec un score d’occupation ≥ 2) dans les ensembles de données de tissus normaux appariés ; (4) les niveaux de CTCF sont plus faibles dans le cancer par rapport à tous les autres ensembles de données (score statistique < 0), (5) les niveaux de CTCF sont plus faibles dans le cancer par rapport aux ensembles de données de tissus normaux appariés (score statistique < 0), (6) signaux de liaison CTCF moyens (RPKM) < 5 dans les ensembles de données de cancer.
-
Sites de liaison CTCF gagnés spécifiques au cancer : (1) fréquence d’occupation ≥ 0,5 (avec score d’occupation ≥ 2) dans les ensembles de données sur le cancer, (2) fréquence d’occupation ≤ 0,2 dans 771 ensembles de données, (3) score d’occupation = 0 dans les ensembles de données sur les tissus normaux appariés, (4) les niveaux de CTCF sont significativement plus élevés dans le cancer par rapport à tous les autres ensembles de données (FDR ≤ 0.01), (5) les niveaux de liaison CTCF sont significativement plus élevés dans le cancer par rapport aux ensembles de données de tissus normaux appariés (FDR ≤ 0,01), (6) signaux de liaison CTCF moyens (RPKM) > 2 dans les ensembles de données du cancer.
Les sites de liaison CTCF spécifiques gagnés et perdus pour chaque type de cancer sont indiqués dans le fichier additionnel 4 : tableau S3.
Quantification de l’accessibilité chromatinienne différentielle
Nous avons utilisé les données traitées de la Réf. qui incluent une matrice des comptes d’insertion ATAC-seq normalisés dans l’ensemble de pics du pan-cancer TCGA pour évaluer l’accessibilité chromatinienne différentielle autour des sites de liaison CTCF. Pour chaque type de cancer parmi BRCA, CRC, LUAD et PRAD, les pics ATAC-seq de l’ensemble des cancers qui chevauchent les sites de liaison CTCF perdus ou gagnés spécifiques au type de cancer ont été utilisés pour les analyses en aval. Le score différentiel ATAC-seq de chaque pic a été quantifié comme le changement de pli de la moyenne des comptes d’insertion ATAC-seq normalisés des échantillons de patients dans le type de cancer correspondant par rapport aux patients dans d’autres types de cancer, et le score différentiel ATAC-seq a ensuite été attribué au pic chevauchant le site de liaison CTCF.
Pour la cohérence, nous avons appliqué la même approche utilisée pour les données ATAC-seq TCGA pour analyser les données ATAC-seq collectées de la lignée cellulaire T-ALL Jurkat et des cellules T CD4+ normales. Une matrice de données a été générée en utilisant le nombre de lectures brutes ATAC-seq sur les sites de liaison CTCF de l’union pour tous les ensembles de données de Jurkat et de cellules T. Une normalisation par quantile a été appliquée à la matrice mise à l’échelle log2 (pseudo compte = 5). Le score différentiel ATAC-seq a été mesuré comme le changement de pli des comptes ATAC-seq normalisés moyens entre les ensembles de données de Jurkat par rapport aux cellules T CD4+ à chaque site de liaison CTCF.
Normalisation des interactions chromatiniennes
Du fait d’une carte de contact Hi-C A = {aij}, le score aij reflète les lectures cartographiées entre deux régions génomiques i et j. Supposons que la taille de la case soit de 5 kb, les régions i et j auront une distance génomique de ∣i – j ∣ × 5kb. Puisque la probabilité de contact entre deux bacs diminue avec l’augmentation de la distance génomique , nous avons normalisé la carte de contact comme suit : pour toute distance génomique donnée dk = k × 5kb, nous quantifions un facteur de normalisation \( {\overline{S}}_{d_k} \) comme la moyenne des interactions entre toutes les paires de bacs ayant la même distance génomique dk dans un même chromosome, par ex, \( {\overline{S}_{d_k}=\left({\sum}_{j-i=k}{a}_{ij}\right)/n \), où n est le nombre total de paires de cases avec la distance dk. Le score d’interaction aij entre deux cases de distance dk a ensuite été normalisé par \( {\overline{S}}_{d_k} \) comme \( {a}_{ij}^{\prime }={a}_{ij}/{\overline{S}}_{d_k} \). En utilisant cette approche, nous avons normalisé la matrice A en \( A^{\prime }=\left\{a}_{ij}^{\prime}\right\} \) au sein de chaque chromosome.
Détection des interactions chromatiniennes différentielles
Nous avons désigné les cartes de contact Hi-C normalisées dans l’ensemble de données sur le cancer et l’ensemble de données normales par C = {cij} et N = {nij}, respectivement. Pour un site de liaison CTCF donné x (avec la coordonnée xc) et une distance génomique prédéfinie L, les interactions chromatiniennes entre x et ses bins voisins non chevauchés de 5 kb avec une distance génomique jusqu’à L sont collectées à partir de C et N respectivement. Plus précisément, les scores d’interaction entre x et ses cases de 5 kb proches dans C sont collectés sous la forme IC = {cij}. , alors que i ou j est égal à ⌊xc/5kb⌋, et 0 < (j – i) × 5kb ≤ L. De même, les scores d’interaction entre x et ses bins 5-kb proches dans N ont été collectés comme IN = {nij}. Un test t de Student bilatéral apparié a ensuite été appliqué sur IC et IN pour quantifier l’interaction différentielle entre le cancer et les cellules normales entourant le site de liaison CTCF x.
Association de la liaison CTCF avec l’expression génique
Au total, 54 types de cellules pour lesquelles les données ChIP-seq CTCF et les données RNA-seq sont publiquement disponibles ont été sélectionnées (fichier additionnel 6 : tableau S5) pour étudier l’association entre la liaison CTCF et l’expression génique pour chaque paire CTCF-gène dans le même chromosome. Pour obtenir le niveau de liaison du CTCF, une matrice de comptage des lectures a été générée en utilisant les lectures par kilobase par million (RPKM) sur les sites de liaison du CTCF de l’union à partir des données ChIP-seq. La matrice de comptage des lectures a été mise à l’échelle par la racine carrée de RPKM, puis normalisée par quantile. Le niveau d’expression des gènes a été mesuré pour chaque gène en utilisant la racine carrée des transcriptions par million (TPM) à partir des données RNA-seq. Pour chaque paire CTCF-gène, nous avons quantifié l’association entre le site CTCF et le gène dans les 54 types de cellules en utilisant le coefficient de corrélation R entre le niveau de liaison CTCF normalisé et l’expression du gène (Fig. 3a). Les paires CTCF-gène ont été considérées comme « hautement corrélées » avec un R2 supérieur à 0,25, par ex, coefficient de corrélation supérieur à 0,5 ou inférieur à – 0,5, et les paires CTCF-gène fortement corrélées contribuent à 1,3 % de toutes les paires CTCF-gène (fichier supplémentaire 1 : Fig. S8a).
Identification des domaines chromatiniens constitutifs limités par CTCF
Pour chaque site de liaison CTCF, nous avons défini son domaine chromatinien associé comme la région génomique qui (1) comprend ce site de liaison CTCF spécifique, (2) est limitée par une paire de sites de liaison CTCF constitutifs avec des motifs d’orientations opposées, et (3) occupe un minimum de 100 kb et un maximum de 1 MB de région de chaque côté du site de liaison CTCF. La figure 3b contient un schéma de la façon dont les domaines chromatiniens constitutifs limités par le CTCF ont été définis.
Détection des changements de méthylation de l’ADN entourant les sites de liaison du CTCF
Les changements de méthylation de l’ADN ont été détectés dans une région de 300 pb centrée sur chaque site de liaison du CTCF. Les régions comportant au moins 3 CpGs couverts par au moins 5 lectures (≥ 5×) dans les deux lignées cellulaires cancéreuses et les tissus normaux correspondants ont été retenues. Une région de 300 pb a été détectée comme étant différentiellement méthylée si les niveaux de méthylation différentiels moyens de tous les CpG (≥ 5×) dans cette région étaient supérieurs à 20 % .
Détection du taux de mutation et du score de motif différentiel
Pour chaque site de liaison CTCF, le nombre brut de mutations a été calculé comme l’occurrence d’événements de mutation dans tous les échantillons/patients à chaque paire de bases unique dans une région de 400 pb centrée sur le site de liaison CTCF. Le taux de mutation pour un groupe de sites de liaison CTCF a été calculé comme le nombre moyen de mutations sur le nombre de sites de liaison CTCF pour chaque paire de bases dans la région de 400 pb.
Le score du motif a été mesuré en évaluant la matrice de poids de position CTCF (Jaspar , Matrix ID : MA0139.1) à une séquence d’ADN de 19 pb centrée sur le motif CTCF ou le site de liaison CTCF en utilisant les rapports de vraisemblance logarithmique (avec la fréquence des nucléotides de fond comme pour A,C,G,T). Le score différentiel de motifs a été calculé en comparant les scores de motifs pour les séquences de référence et les séquences mutées.
Analyse des motifs de séquences d’ADN
L’analyse d’enrichissement des motifs de séquences d’ADN a été réalisée en utilisant MDSeqPos (version 1.0.0) sur Cistrome avec les paramètres par défaut (-cisrome -Homo Sapien ou Mus musculus). Les analyses de novo des motifs ont été réalisées en utilisant HOMER (version 4.10) avec le module findmotifs.pl et MEME (version 5.1.1) avec les paramètres suivants : meme -dna -mod zoops -maxw 20 -evt -0,01.
Identification des régions à interaction différentielle intra-domaine de CTCF
Pour un ensemble donné de sites de liaison CTCF, les changements d’interaction chromatinienne entre un site CTCF et chacune de ses bins intra-domaine non chevauchés, mesurés à partir des cartes de contact Hi-C normalisées dans les cellules cancéreuses par rapport aux cellules normales appariées, ont été collectés pour chacun des sites de liaison CTCF (fichier additionnel 1 : figure S14b). Les régions présentant des interactions réduites (log2 FC < -1, interaction log2 moyenne > 0) avec les sites de liaison CTCF perdus spécifiques au cancer, et les régions présentant des interactions accrues (log2 FC > 1, interaction log2 moyenne > 0) avec les sites de liaison CTCF gagnés spécifiques au cancer ont été utilisées pour l’analyse d’enrichissement des facteurs de transcription (TF) en aval.
Analyse d’enrichissement des facteurs de transcription
Une version révisée de l’algorithme BART a été utilisée pour l’analyse d’enrichissement des TF. En bref, une collection de sites hypersensibles à la DNase I de l’union (UDHS) a été précédemment curée comme un répertoire de tous les éléments cis-régulateurs candidats dans le génome humain, et 7032 ensembles de données ChIP-seq ont été collectés pour 883 TF , chaque TF ayant un ou plusieurs ensembles de données ChIP-seq provenant de plusieurs types de cellules ou conditions. Un profil binaire a été généré pour chaque TF sur l’UDHS, indiquant si la TF a au moins un pic de l’un de ses ensembles de données ChIP-seq situé dans chacun des UDHS. Une analyse d’enrichissement de la liaison a été appliquée pour chaque TF en comparant la liaison de la TF sur un sous-ensemble d’UDHS chevauchant les régions génomiques sélectionnées par rapport à la liaison de la TF sur les UDHS. La valeur p a été obtenue en utilisant le test exact de Fisher à deux intervalles.