dataesr/scientific-paragraphs-categorization huggingface dataset:
1d  and 2b, c ). Next, we assessed the effect of DNA copy number alterations on TBX2 expression levels using an ANOVA analysis in the NRC tumor dataset (n = 218, GSE85047). Only in high-stage disease    4 ) Ewingsvsarcoma (  9 ) Glioma (43)  Haematopoietic_neoplasm (  50 ) Lymphoid_neoplasm (122) Melanoma (58)  Mesothelioma (  9 ) Neuroblastoma (17)  Osteosarcoma (  9 ) Medulloblastoma (  4   17 ) Medulloblastoma(  4 ) Urinary_tract(  28 ) Ewings_sarcoma(  12 ) Prostate(  8 ) Mesothelioma(  11 ) Melanoma(  63 ) Glioma(  66 ) Chondrosarcoma(  4 ) Ovary(  55 ) Lung_small_cell(  54 ) Upper_aerodigestive(  33 ) Liver(  29 ) Soft_tissue(  20 ) Osteosarcoma(  10 ) Other(  8 ) Lung_NSC(  136 ) Kidney(  37 ) Breast(  60 ) Thyroid(  12 ) Etsophagus(  27 ) AML(  39 ) Endometrium(  28 ) Pancreas(  46 ) CML(15) B-cell_ALL(  13 ) Leukemia_other(  5 ) Colorectal(  63 ) Multiple_myeloma(29) B-cell_lymphoma_other(  16 ) Lymphoma_Hodgkin(  13 ) Bile_duct(  8 ) T-cell_lymphoma_other(  11 ) Stomach(  39 ) Lymphoma_DLBCL(  18 ) T-cell_ALL(  16 ) Giant_cell_tumour   16 ) Meningioma(  3 ) Lymphoma_Burkitt(  11 ) giant_cell_tumour(  3 ) Leukemia_other(  5 ) Urinary_tract(  28 ) Lymphoma_Hodgkin(  13 ) Colorectal(  63 ) Thyroid(  12 ) B-cell_ALL(  13 ) NA(387) Lymphoma_DLBCL(  18 ) AML(  39 ) Multiple_myeloma(  29 ) Other(  8 ) Prostate(  8 ) Esophagus(  27 ) Upper_aerodigestive(  33 ) stomach(  39 ) Ewings_sarcoma(  12 ) Endomet Glioma(  66 ) Lung_NSC(  136 ) soft_tissue(  20 ) osteosarcoma(  10 ) Ovary(  55 ) CML(  15 ) T-cell_ALL(  16 ) Pancreas(  46 ) Lung_small_cell(  54 ) T-cell_lymphoma_other(  11 ) Liver(  29 ) Kidney(  37 ) Melanoma(  63 ) Mesothelioma(  11 ) Medulloblastoma(  4 ) Breast(  60 ) Bile_duct(  8 ) Neuroblastoma(  17   ) and increased TBX2 copy number due to 17q gain (ANOVA, n = 218, p = 6.
To gain further insight into the TBX2-controlled regulatory network, we assessed TBX2 DNA occupancy by ChIP-sequencing and ATACsequencing in the NB cell line IMR-32. A total of 557 significant (adj. P. val < 0. 05) TBX2 binding sites were identified and motif analysis confirmed enrichment for a TBX motif (AGGTGTGA, p = 1e-41), supporting the validity of our ChIP-seq data (Supplementary Data 1). In total, 81, 28, and 94% of TBX2 binding sites in IMR-32 respectively overlap H3K27ac, H3K4me3, and ATAC-sequencing peaks (Fisher test p < 2. 2e-16, Fig. 3a, b ), which confirms the binding of TBX2 to active promotor and enhancer regions. Moreover, respectively 41 and 30% of the TBX2 ChIP-seq peaks are found intergenic or are annotated to lncRNAs (Supplementary Fig. 3a, b ), and 19% overlap with the SEs annotated in the cell line IMR-32 (Fisher test p < 2.
4b , Supplementary Data 4). In addition to FOXM1 target genes, FOXM1 itself is also downregulated upon TBX2 KD in the IMR-5/75 and CLB-GA cell lines (Fig. 4d ). Using a similar approach as for MYCN amplified IMR-5/75 cells, we also investigated the TBX2 transcriptional regulated network in the non-MYCN amplified CLB-GA NB cells and observed similar downstream targets and enriched genes sets as those observed in IMR-5/75 upon TBX2 KD (Supplementary Fig. 4a-d , Supplementary Data 3 and 5), with the E2F/FOXM1 axis being most prominent. To further validate these results, we performed correlation analysis in publicly available transcriptome data of 283 primary NB tumors (NRC, GSE85047). Supporting the TBX2 KD data described above, GSEA showed enrichment for cell cycle, DNA repair and DNA replication as well as chromatin architecture among the genes positively correlated with TBX2 expression levels in the NB tumor data set (FDR < 0.
A recent report suggested that TWIST and HAND2 would cooperate with MYCN through binding of enhancers of a set of developmental genes  16. In view of these findings, we also decided to further study the TBX2-MYCN interrelationship. We first obtained indirect evidence for a functional relationship between MYCN and TBX2 from several datasets: (1) significant enrichment for publicly available MYC (N) target/signature gene sets in the MYCN amplified IMR-5/75 and MYCN non-amplified cell line CLB-GA upon TBX2 KD (Figs. 4a, 6a and Supplementary Fig. 4d);  (2) significant correlation of TBX2 and MYCN mRNA and protein expression in NB tumors (Supplementary Fig. 3f ) and NB cell lines (Supplementary Fig. 6a ), respectively; (3) correlation between MYCN activity and TBX2 shRNA signature scores (established in both IMR5/75 and sets and finally (5) significant dynamical downregulation of the shTBX2 signature was noted during Tg(TH-MYCN) driven NB formation in transgenic mice, at one, two and six weeks after birth  38  (Supplementary Fig.
6c ). To further experimentally explore this presumed cooperation between TBX2 and MYCN, we assessed the effects of TBX2 KD in the presence of high versus low MYCN levels. Using this approach, we observed a stronger decrease in cell proliferation and increased G1-phase arrest when combining TBX2 and MYCN KD together (Fig. 6b, c ). Next, we performed RNA-sequencing to further explore the transcriptional effects of individual versus combined MYCN and TBX2 KD and confirmed synergistic effects on expression levels of gene sets implicated in cell cycle and the DREAM-E2F-FOXM1 complex, as well as publicly available MYCN signatures (Fig. 6d ). Using iRegulon analysis on the enforced affected genes, motif enrichment of DREAM complex core members, such as FOXM1, E2F4, and MYBL2 was observed in the additively downregulated genes, while ChIP-seq targets for EP300 and NANOG (public datasets) were enriched in the upregulated genes (Supplementary Fig.
Real-time quantitative PCR. Total RNA was extracted using miRNeasy kit (Qiagen) according to the manufacturer's instructions, including on-column DNase treatment, and concentration was determined with the Nanodrop (Thermo Scientific). cDNA synthesis was performed using the iScript Advanced cDNA synthesis kit from BioRad. PCR mix contained 5 ng of cDNA, 2. 5 ul SsoAdvanced SYBR qPCR supermix (Bio-Rad) and 0. 25 µl forward and reverse primer (to a final concentration of 250 nM, IDT) and was analysed on the LC-480 device (Roche) for RT-qPCR cycling. Expression levels were normalized using expression data of 3 stable reference genes out of 5 reference genes tested (SDHA, YWHAZ, TBP, B2M, HPRT1) and analyzed using qBasePlus software (http://www. biogazelle. com). All primer pair sequences can be found in Supplementary Table  2.
Genes were only retained if they were expressed at counts per million (cpm) above one in at least four, eight or twelve samples for the drugging, shTBX2 in IMR-5/75 and shTBX2 in CLB-GA datasets respectively. Counts were normalized with the TMM method (Rpackage edgeR), followed by voom transformation and differential expression analysis using limma (R-package limma). A general linear model was built with the treatment groups (drugging or KD) and the replicates as a batch effect. Statistical testing was done using the empirical Bayes quasi-likelihood F-test. GSEA  59  was performed on the genes ordered according to differential expression statistic value (t). Signature scores were conducted using a rank-scoring algorithm  65. Limma voom barcodeplots were used for visualization of gene set enrichment.
As secondary antibody, we used HRP-labeled anti-rabbit (7074 S, Cell Signalling, 1:10,000 dilution) and anti-mouse (7076P2, Cell Signalling, 1:10,000 dilution) antibodies. Antibodies against Vinculin (V9131; Sigma-Aldrich, 1:10,000 dilution), alpha-Tubulin (T5168, Sigma-Aldrich, 1:10,000 dilution) or bèta actine (A2228; Sigma-Aldrich, 1:10,000 dilution) were used as loading control. All antibodies were diluted in milk/TBST (5 % non-fat dry milk in TBS with 0. 1 % Tween-20). Binding of the antibodies with the membrane was evaluated using the SuperSignal West Dura Extended Duration Substrate (Thermo Scientific). Pictures were taken with the ChemiDoc-It Imaging System (UVP) using the VisionWorks analysis software (UVP), quantification of the blots were performed using ImageJ. Uncropped scans of the blots used in the main figures can be found in the Supplementary Figure  9.
Masking was done using the IncuCyte® ZOOM Software. For cell cycle analysis, 7 × 10 5 cells were seeded in a T25 in complete medium and transduced with TBX2 shRNAs and controls and selected with puromycin, as described above. Cells were trypsinized and washed with PBS. The cells were resuspended in 300 µl cold PBS and while vortexing, 700 µl of 70% ice-cold ethanol was added dropwise to fix the cells. Following incubation of the sample for minimum 1 h at -20 °C, cells were washed in PBS and resuspended in 500 µl PBS with RNase A to a final concentration of 0. 25 mg/ml. Upon 1 h incubation at 37 °C, 20 µl Propidium Iodide solution was added to a final concentration of 40 µg/ml. Samples were loaded on a BioRad S3 TM Cell sorter and analysed with the Dean-Jett-Fox algorithm for cell-cycle analysis using the FlowJo ® software package.
Cell viabilities of DMSO-treated cells were set to 100%. Excess over Bliss scores (score >0 indicates drug synergy, <0 indicates drug antagonism) were calculated and subsequently visualized using the R package synergyfinder. Values more than 100% were replaced by 100% for the concentration range experiments of SH-SY5Y and SK-N-BE(2c) to fulfill the criteria for Bliss score calculation. IncuCyte® assay for assessment of drug synergism. Neuroblastoma cells (CLB-GA, IMR-32, IMR-5/75, SH-SY5Y, Kelly, SK-N-BE(2)-C and SK-N-AS) and MCF-7 cells were seeded in complete medium and a 1:200 dilution of IncuCyte® Annexin V Red Reagent for apoptosis (Essen Bioscience) in 96-well tissue culture plates (Corning costar 3596) in triplicate at 30% confluency and allowed to recover overnight. The next day, cells were treated with the indicated concentrations of JQ1 (BPS Bioscience), THZ1 hydrochloride salt (Medchem express) or Panobinostat (LBH589, Selleck Chemicals), or combinations thereof.
Cell proliferation and cell death (AnnexinV positivity) were assessed continuously for 68-96 h after treatment by using the IncuCyte® Live Cell imaging system (Essen BioScience), by taking photos every 3-4 h of the whole well. Masking was done using the IncuCyte® ZOOM Software. Possible synergism was calculated using the Bliss method. Chromatine immunoprecipitation (ChIP) assay. A total of 3-5 × 10 7 cells were crosslinked in complete medium with 1% formaldehyde while shaking for 7 min at room temperature. Crosslinking was quenched with 125 mM glycine (Sigma-Aldrich), cells were washed twice with PBS and stored at -80 °C. Cells were lysed and sonicated with the M220 Covaris for 15 min to obtain 200-300 bp long fragments. Chromatin fragments were immunoprecipitated overnight using 1 µg for every 10 7 cells of the following antibodies anti-TBX2 (SC-17880), anti-H3K4Me3 (Ab8580), anti-H3K4me1 (Ab8895), anti-H3K27ac (Ab4729) and anti-MYCN (SC-53993).
2. 3. 1), quality of the raw sequencing data was evaluated using FastQC and adapter dimers were removed using cutadapt when necessary. Peak calling was performed using MACS2 taking a q value of 0. 05 as threshold. IGV was used for visualization purposes with the by IGVtools generated tdf files. Bedgraph files generated by MACS2 were converted to bigwig files and used as input for the ChIP-tracks in the figures, further processed using the R package Sushiplot. Homer 68 was used to perform motif enrichment analysis, taking 200 bp around the peak summit as input. Gene annotation tables, tag density plots, heatmaps and overlap of peaks maps were generated using Homer and the R packages EnrichedHeatmap and ChIPpeakAnno, with default parameters. Super-enhancer analysis. ChIP-sequencing data of three different datasets were used for SE calling.
(a) H3K27ac ChIP-sequencing for the cell lines CLB-GA (rep1), NGP, IMR-5/75 and IMR-32 (rep1) was performed as reported above with the H3K27ac antibody (ab4729, Abcam). (b) H3K27ac ChIP-sequencing for the cell lines SK-N-AS (rep1), CHP-212 (rep1), GI-M-EN (rep1) and IMR-5 was performed by Zymo Research using between 5 × 10 6 and 10 × 10 6 frozen cross-linked cells and pull down with the anti-H3K27ac antibody, as previously reported by Henssen et al  69. (c) The public available H3K27ac ChIP-seq fastq files for CHP-212, CLB-MA, CLB-PE, CLB-GA (rep2 and rep3), GI-C-AN, GI-M-EN (rep2), IMR-32 (rep2), LA-N-1, N206, NB-EBc1, SH-SY5Y, SJNB-1, SJNB-12, SJNB-6, SJNB-8, SK-N-AS (rep2), SK-N-BE(2)-C, SK-N-DZ, SK-N-FI, TR14 (GSE90683  15  ) and MCF-7 (GSE69112 70 ) were obtained using the SRAtoolkit and mapping to ref genome hg19 was done using bowtie2 (v.
1 Mb-qter). Among the top 500 SEs ranked according to the median rank, six regions were prioritized based on the presence of an SE cluster in min 20 NB cell lines. Ensembl genes (bioMart) in a 500 Kb range around the SE cluster were annotated to the different SE clusters. Of the annotated transcription factors  71, 72  , only those with a H3K27ac activity called peak at their TSS (Macs2, q < 0. 05) were selected, as SEs are supposed to highly activate their target genes  13. Genes were only selected when they had a H3K27ac marked TSS in all cell lines harboring the SE. Both IGV and the R package Sushi were used to visualize the SE regions in the different cell lines. 4C-sequencing and data analysis. Preparation of 4C samples was performed as described previously  73  , starting with 10 7 separated single cells that were crosslinked with 2% formaldehyde.
DNA was digested using DpnII as a primary and Csp6I as a secondary restriction enzyme. 50 ng of the final 4C template of two NB cell lines (CLB-GA and SK-N-AS) was used for the amplification step with three different primer sets (Supplementary Table  2 ) designed for the TBX2 locus. PCR products were purified and pooled for multiplexed sequencing after adding 40% PhiX to increase sequence complexity. Sequencing was performed on the HiSeq 2000 V4 platform (50 bps single reads). Sequencing reads from the 4C library were demultiplexed using the demultiplex. py script (https://gist. github. com/meren/ 7632184) and aligned using bowtie2 (v2. 3. 1) with parameters-sensitive-timeend-to-end -q -p 2. Aligned data were stored together with the metadata into a FourCSeq object  74. The aligned libraries were subjected to a QC according to criteria published elsewhere  73.
Cis-and trans-interactions were determined using a custom script implementing the method described in  73, 74. Briefly, a z-score is computed on the binarized signal (i. e. , counting all fragments covered by reads) using sliding windows (w = 100 fragments, W = 3000 fragments), and an FDR is computed from this z-score to identify significant interactions. We use a FDR of 1% for cis-interactions and 0. 5% for trans-interactions. Domainograms showing the signal intensity of different ranges are plotted using a custom R-script. Bigwig files for visualization are generated using a smoothing over 21 consecutive fragments, and normalized to the total library size. TBX2 expression analysis in NB tumors and cell lines. TBX2 expression analysis was performed on 283 NB tumors for which copy number (n = 218), mRNA expression (n = 283) and patient survival (n = 276) data were available from the Neuroblastoma Research Consortium (NRC, GSE85047), which is a collaboration between several European NB research groups.
Additionally, the NB dataset from Su et al. (n = 489, GSE45547) was used as validation cohort  25. All statistical analyses were performed using R (version 3. 3. 0). Genes were ranked according to the correlation of their expression (Spearman correlation) with TBX2 expression levels. Next, pre-ranked GSEA was performed as previously described  59  with the MsigDB c5. bp. v4. 0 geneset and clustering of enriched genesets (only when FDR < 0. 01, and R > 4/R < -3) was performed using the tool EnrichmentMap in Cytoscape (version 3. 2. 1). TBX2 expression in normal tissue, tumor, and cell line datasets was visualized on the R2 genomics analysis and visualization platform (http://r2. amc. nl/). DNA copy number analysis. The presence of focal gains encompassing the TBX2 locus was investigated in published copy number profiling data of 556 high-risk tumors  29.
Acknowledgements The authors would like to thank  W. Van Loocke ,  L. Mus ,  G. Dewyn ,  E. Sanders ,  J. Van Laere ,  E. De Smet ,  S. Claeys  and  E. Janssens  for technical assistance and  Prof. Stein Aerts  for providing the transcription factor list used in iRegulon and iCistarget. This research was supported by the following funding agencies: the  Belgian Foundation against Cancer  (project  2015-146 ) to F. S. P. ,  Ghent University  ( BOF10/GOA/019  and  BOF16/GOA/23 ) to F. S. P. , the  Belgian Program of Interuniversity Poles of Attraction (IUAP Phase   VII-P7/03 ) to F. S. P. , the  Fund for Scientific Research Flanders  (Research projects  G053012N ,  G050712N , and  G051516N  to F. S. P. ,  G021415N  to K. D. P and F. S. P. ), ' Kom op tegen Kanker ' (Stand up to Cancer) the  Flemish cancer society  (Research grant to F.
S. P. ), the  European Union   H2020  ( OPTIMIZE-NB GOD9415N  and  TRANSCAN-ON THE TRAC GOD8815N  to F. S. P. ), 'Kinderkankerfonds' (Research grant to F. S. P. ),  Olivia Fund  to F. S. P. and  Villa Joep  to F. S. P. FWO grants were supporting the work of B. D. C. , S. L. N. , S. C. L. , K. D. R. , C. N. N. , S. V. H. , and C. V. N ; K. V. B. and K. D. N. are supported by a  BOF  grant. D. R. M. is supported by a  IWT grant. A. G. H. is supported by the  Deutsche Forschungsgemeinschaft (DFG, Germ and Research Foundation -398299703 ).
Data availability The RNA-sequencing, ChIP-sequencing and ATA-sequencing datasets generated during this study were deposited in the ArrayExpress database at EMBL-EBI (www. ebi. ac. uk/ arrayexpress) with accession numbers: E-MTAB-6570, E-MTAB-6562, E-MTAB-6568, E-MTAB-6567 and E-MTAB-7025. Data that supports the findings of this study are available from the Neuroblastoma Research Consortium (NRC, GSE85047), Su et al. [GSE45547]  25  and Depuydt et al. [GSE103123] 29.
Annotation of the >116 million unigenes (Methods and Supplementary Fig. 1a ) revealed that we could assign a taxonomy level (from "cellular organism" to species name) to only 48. 3% of the unigenes (Fig. 2a  and Supplementary Fig. 1b ). By mapping the unigenes onto known gene annotations from marine genomes, we found a mean value of 2. 20 (s. d. = 0. 47) unigenes per gene (Methods and Supplementary Data 2). We then estimated the number of distinct transcriptomes (originating from different species) that were present in the catalog by counting the mean number of copies of conserved ribosomal protein genes, which indicated that the catalog contains genes from 8823 (s. d. = 1532) different organisms (Supplementary Data 3). These values indicate that the unigenes are derived from around 53 (44-68)  million genes, with a mean of 6014 (4226-9223) genes per sampled organism (Supplementary Data 4).
2b ), confirming the efficiency of the fractionation-based approach. We observed 1. 13% of unigenes that are affiliated to prokaryotes. These were not removed from the catalog, as they can be true nonpolyadenylated transcripts from this group, or alternatively to the low level of eukaryotic annotations with respect to prokaryotes in reference databases, or to horizontal gene transfers. Our metatranscriptomic data also captured transcripts (or RNA genomes) of viruses actively infecting their eukaryotic hosts. Their activities were found to be pervasive across the geographic and organismal size ranges examined in this study. Of the taxonomically assignable unigenes, 33,870 (0. 06%) were predicted to be of eukaryotic virus origin, the vast majority of which (86%) originated from nucleocytoplasmic large dsDNA viruses (NCLDVs)  21  (Fig.
We also observed that the relative expression of photosynthesis genes (seen through chlorophyll-binding proteins) vs. proteorhodopsins (Bac_rhodopsin Pfam domain corresponding to type-I rhodopsins  22, 23  ) showed a strong preference for photosynthesis in groups dominated by autotrophs, supporting that rhodopsin is not a major way of using light energy in these groups in natural conditions (Supplementary Fig. 5a ). To further investigate the distribution of the expression of the rhodopsins present in the catalog, we isolated all the unigenes bearing a Bac_rhodopsin Pfam domain. We added to the dataset 2112 proteins-mainly from fungi (40%), bacteria (35%), and archaea (18%)-from public databases and 2538 eukaryotic protein sequences from MMETSP  14. Protein sequences from the 71,576 unigenes carrying the Bac_rhodopsin Pfam domain were aligned and clustered with reference sequences to study their diversity (Methods section).
3c ). Gene novelty. The majority (51. 2%) of unigenes currently have no matches in public sequence databases, which limits the insights that can be derived from the gene catalog. Some sequences may be derived from non-coding genes or non-coding portions of coding genes, very short open reading frames, parts of genes where only another region is functionally known, or completely new open reading frames. To distinguish between these possibilities and better classify the catalog, we clustered all the unigenes according to a nucleic acid similarity threshold of >70% (Methods; Supplementary Fig. 6a ). Despite its size, the gene catalog is not saturated, and accordingly we observed that 59. 6% of unknown unigenes (UU) and 39. 8% of known unigenes are represented by singletons (Fig. 4a, b ). The clusters may thus be considered as being representative of gene family (GF) content of the catalog, with most singletons likely being derived from  Among the highly expressed Pfam domains (Supplementary Fig.
4 ), those with specific patterns are shown. The relative expression of Pfam domains in the four filter sizes (left panel) and the contribution of each taxonomic group to the total expression of the Pfam domain (right panel) are shown as an average of all Tara Oceans SRF and DCM samples. O/U, other or unassigned. b Unrooted phylogenetic tree of type-I rhodopsin subfamilies (PF01036) obtained using sampling of 300 sequences of the three largest MCL clusters (see details in Supplementary Fig. 5b ). The vertical size of the triangles represents the number of unigenes in each cluster (explicitly indicated in white) and their width represents the maximum branch length of 95% of sequences in the cluster. Taxonomic assignments of reference sequences (inner ring) and unigenes (outer ring) are indicated for each cluster with the color code of a.
4d  and Supplementary Fig. 7a ). This pattern was not due to higher mean expression levels of ftGFs or tGFs that would render them more detectable than nGFs (Fig. 4e ). We conclude that the gene novelty detected corresponds to families that are present in fewer environments, yet are not less expressed than known families. Moreover, nGFs generally represent smaller GFs (6. 3 unigenes per cluster) than fGFs (8. 9) and ftGFs (11. 4), suggesting that nGFs are conserved in a smaller range of species than characterized GFs (Fig. 4a  and Supplementary Fig. 6b ), or that they are present in less abundant taxonomic groups. It has been previously suggested that newly discovered genes are either biased taxonomically (which restrains their presence in databases), or that they correspond to genes that are necessary only in some conditions, potentially related to the adaptation of organisms to specific environments  26.
The current database of protein domains such as Pfam 27 contains 16,712 different domains of known and unknown functions, whereas we detected 11,431 protein groups of nGFs, and 22,072 groups of tGFs based only on Clustering of the largest families, indicating the high discovery rate of new conserved domains that could be used to derive a more exhaustive list of conserved domains within eukaryotes. In summary, we have found that UUs can be part of known GFs but that a large proportion are predicted to be novel proteincoding genes. As they are distributed less globally than known functions, their extent remains to be evaluated, although we have shown here that they represent a highly significant portion of the gene repertoire of eukaryotic plankton. The environmental footprint of gene expression in phytoplankton.