dataesr/scientific-paragraphs-categorization huggingface dataset:
Spoligotyping and database comparison Spoligotyping was carried out according to the protocol previously described by Kamerbeek et al  [9] , using a commercially available kit (Ocimum Biosolutions Inc. , Hyderabad, India). M. tuberculosis H37Rv and M. bovis BCG DNAs were used as parallel positive controls and distilled water as a negative control. The patterns obtained were compared by using the SITVIT2 proprietary database of the Pasteur Institute of Guadeloupe, which is an updated in-house version of the publicly released SITVITWEB database  [10] ; available online at http://www. pasteur-guadeloupe. fr:8081/SITVIT_ONLINE/]. In this database, Spoligotype International Type (SIT) designates an identical pattern shared by two or more patient isolates, whereas ''orphan'' designates patterns reported for a single isolate that does not correspond to any of the strains recorded in the database repository.
The BioNumerics software Version 6. 6 (Applied Maths NV, Sint-Martens-Latem, Belgium) was used to build minimum spanning trees (MST) based on spoligotyping data. MST is an undirected network in which all of the isolates are linked together with the fewest possible linkages between nearest neighbors, represented by branches (continuous vs. dashed and dotted lines). Note that the length of the branches represents the distance between patterns while the complexity of the lines (continuous, gray dashed and gray dotted) denotes the number of allele/spacer changes between two patterns: solid lines, 1 or 2 or 3 changes (thicker ones indicate a single change, while the thinner ones indicate 2 or 3 changes); gray dashed lines represent 4 changes; and gray dotted lines represent 5 or more changes. Spoligoforest tree based on Fruchterman-Reingold algorithm was drawn using the SpolTools software  [14, 15].
Contrary to the MST, the Spoligoforest trees are directed graphs which only evolve by loss of spacers, represented by directed edges between nodes, with the arrowheads pointing to descendant spoligotypes. In this representation, the heuristic approach used selects a single inbound edge with a maximum weight using a Zipf model. Solid black lines link patterns that are very similar, i. e. , loss of one spacer only (maximum weight being 1. 0), while dashed lines represent links of weight comprised between 0. 5 and 1, and dotted lines a weight less than 0. 5. In a Spoligoforest, nodes are not necessarily all connected (indeed, in case of too many changes between two strains, there are no edges linking them. GraphViz software (http://www. graphviz. org) was used to color the strains based on their lineages on the Spoligoforest trees  [16].
Worldwide distribution analysis The predominant MTBC SITs/lineages obtained in our study (SITs containing. 1% i. e. , 3 or more isolates) were screened for their worldwide distribution analysis using the SITVIT2 database country-wise as well as macro-geographically as United Nations sub-regions (http://unstats. un. org/unsd/methods/m49/ m49regin. htm): AFRI (Africa), AMER (Americas), ASIA (Asia), EURO (Europe), and OCE (Oceania), subdivided in: E (Eastern), M (Middle), C (Central), N (Northern), S (Southern), SE (South-Eastern), and W (Western). In this classification scheme, CARIB (Caribbean) belongs to Americas, while Oceania is subdivided in 4 sub-regions, AUST (Australasia), MEL (Melanesia), MIC (Micronesia), and POLY (Polynesia). Note that Russia was attributed a new subregion by itself (Northern Asia) instead of including it among the rest of Eastern Europe.
Mapping in neighboring countries We also compared MTBC population structure in our study vs. results retrieved from the SITVIT2 database for neighboring countries in South Asia n54271 strains (India n52647, Bangladesh n51090, Pakistan n5385, Nepal n549, Sri-Lanka n5100), East Asia n55216 strains (China n53072, Japan n51569, South Korea n5152, Taiwan n5423), and South East Asia n53059 strains (Indonesia n5422, Myanmar n5310, Malaysia without our present study n5749, Philippines n5155, Singapore n5333, Thailand n5307, Vietnam n5783).
Statistical analysis Statistical analyses were performed using the Statistical Package for the Social Sciences (SPSS) version 18. 0 and STATA version 10. 1. Associations between M. tuberculosis clades with patient demographic and drug resistant criteria were assessed using Chi-square analysis and Fisher's exact test. P-values of #0. 05 were considered statistically significant and the 95% confidence interval (95% CI) was calculated. Mean and standard deviation (SD) were calculated for age of the patients.
0%) of isolates were not available due to several factors including data loss and contamination of the culture during the process. Spoligotyping of the study sample (n5220 isolates) as well as the detailed results obtained including demographic, epidemiologic, drug-resistance, and genotyping information are summarized in S1 Table. Spoligotyping revealed 77 different patterns: 22 corresponded to orphan patterns not yet reported in the SITVIT2 database, while 55 patterns containing 198 isolates were assigned a Spoligo International Type (SIT) designation in the database (the latter included 6 ''newly created'' SITs containing 11 isolates from this study). Refer to S1 Table for a listing of orphan patterns (highlighted in blue), and Table  1  for a description of 55 SITs. Note that SITs followed by an asterisk indicate "newly created'' SITs due to 2 or more strains belonging to an identical new pattern within this study or after a match with an orphan in the database (SITs followed by number of strains: SIT3993* this study n51, MYS n51; SIT3994* this study n52, MYS n51; SIT3994* this study n52; SIT3995* this study n52; SIT3996* this study n52; SIT3997* this study n52; 3998* this study n52.
82%, and clustered isolates to 161/220 or 73. 18%). If the study sample was classified by clades (and not SITs), the predominant M. tuberculosis clades belonged to the ancestral EAI clade (56. 4%) followed by Beijing (28. 6%), Haarlem (4. 0%) and T (4. 0%). The single M. bovis clade identified in this study was isolated from a 14 year old Malay boy with pulmonary TB; unfortunately, DST patterns were not available for this strain. Table  3  summarizes available demographic, epidemiologic and DST results in function of clades, and one may notice that the proportion of drug resistant strains was higher among isolates belonging to the Beijing lineage (8 drug resistant isolates vs. 29 pansusceptible isolates) as compared to the isolates belonging to EAI lineage (9 Lineage designations according to SITVITWEB rules; ''Unknown'' designates patterns with signatures that do not belong to any of the major lineages described in the database.
c Clustered strains correspond to a similar spoligotype pattern shared by 2 or more strains ''within this study''; as opposed to unique strains harboring a spoligotype pattern that does not match with another strain from this study. Unique strains matching a preexisting pattern in the SITVIT2 database are classified as SITs, whereas in case of no match, they are designated as ''orphan'' (see Table  1 ). doi:10. 1371/journal. pone. 0114832. t001 drug resistant isolates vs. 83 pansusceptible isolates). Apart from that, almost all clades were reported higher in male group than female group except for H clade. Nevertheless, statistical analyses of M. tuberculosis clades by patient demographic and drug resistance criteria showed that neither differences between and clades nor those between drug resistance and clades were statistically significant (p.
0. 5). On the other hand, the association between city of isolation and clades was statistically significant; the proportion of Haarlem and T lineages being higher in Kuala Lumpur than Kelantan (p,0. 01). Investigation of phylogenetic relationships between MTBC strains by Spoligoforest and minimum spanning tree (MST) is shown in Figs. 1A  and 1B , Worldwide distribution is reported for regions with more than 3% of a given SITs as compared to their total number in the SITVIT2 database. The definition of macro-geographical regions and sub-regions (http://unstats. un. org/unsd/methods/m49/m49regin. htm) is according to the United Nations; Regions: AFRI (Africa), AMER (Americas), ASIA (Asia), EURO (Europe), and OCE (Oceania), subdivided in: E (Eastern), M (Middle), C (Central), N (Northern), S (Southern), SE (South-Eastern), and W (Western).
Furthermore, CARIB (Caribbean) belongs to Americas, while Oceania is subdivided in 4 sub-regions, AUST (Australasia), MEL (Melanesia), MIC (Micronesia), and POLY (Polynesia). Note that in our classification scheme, Russia has been attributed a new sub-region by itself (Northern Asia) instead of including it among rest of the Eastern Europe. It reflects its geographical localization as well as due to the similarity of specific TB genotypes circulating in Russia (a majority of Beijing genotypes) with those prevalent in Central, Eastern and South-Eastern Asia. b The 3 letter country codes are according to http://en. wikipedia. org/wiki/ISO_3166-1_alpha-3; countrywide distribution is only shown for SITs with $3% of a given SITs as compared to their total number in the SITVIT2 database. (A) Spoligoforest tree drawn using the SpolTools software (available through http://www.
2 ) highlighted that EAI lineage is the most prevalent lineage in Malaysia, Philippines, Singapore, Myanmar, Thailand, Bangladesh and Sri Lanka, while the Beijing lineage is the most prevalent in China, Japan, South Korea, Taiwan, Vietnam and Thailand. Interestingly, the MTBC population structure observed for Malaysia in the present study is almost similar to the one recorded for a total of 749 MTBC strains recorded in the SITVIT2 database. Last but not least, particularly noticeable was the ongoing evolution of the SIT745 in Malaysia which was closely connected to a network of 8 other derived spoligotype patterns (Fig. 3 ), all showing a remarkable phylogeographical specificity for Malaysia. The parental SIT745 pattern and the locally evolved descendant strains shared a specific signature characterized by absence of spacers 37, 38, and 40.
g. , it accounted for 70% of TB cases in Hong Kong  [19] , 80. 4% in Jiangsu province, China  [20] , and 77. 7% among homeless TB patients in Osaka City, Japan  [21]. A similar distribution pattern was also confirmed for China, Japan and Vietnam in the SITVIT2 database (Fig. 2 ). In contrast, the prevalence of the Beijing genotype is low in Latin America as compared to Asian and East-European countries  [22] [23] [24] [25] [26] [27]. Although the Beijing genotype strains are strongly associated with drug resistance and specifically multiple drug resistant (MDR) and extensively drug resistant (XDR) TB in various geographical settings  [20, [28] [29] [30] , this association holds true particularly in areas where the proportion of the Beijing genotype strain is actively increasing (reviewed in  [31]  and references therein).
Interestingly, SIT745/EAI1-MYS sublineage is tentatively identifiable by the absence of spacers 37 and 38 (a signature common to EAI3-IND sublineage) in conjunction with the absence of spacer 40 (a signature common to EAI1-SOM sublineage). In SITVIT2 database, such strains have only been reported from Malaysia and India, although majorities (. 95%) are reported in Malaysia (Table  2 ). Now such strains should be ideally reconfirmed using supplementary genotyping markers and compared to the prototypes of EAI1-SOM and EAI3-IND sublineages to corroborate our findings. This should not come as a surprise since the endemic nature of EAI lineage in Southeast Asia has been well proven with the existence of specific sublineages, e. g. EAI2-Manila with phylogeographical specificity to the Philippines and EAI4-VNM to Vietnam  [10].
1A ). EAI is the predominant lineage in Malaysia, identified in more than half of overall strains followed by Beijing lineage (Table  3 ). This finding is in line with the previous data from Malaysia reported in the SITVIT2 database (Fig. 2 ). It shows that EAI lineage has high transmission of TB in our study population. One may notice the majority or orphan strains in this study were emerging from this lineage (Fig. 1 ). Even though EAI lineage was the most predominant strain in this study, there was no significant association with gender and drug susceptibility pattern. However, a study in Montreal reported that EAI lineage was associated with lower rates of TB transmission and less likely to cause severe forms of TB, though it was associated with a higher proportion of TB lymphadenitis and extra-pulmonary TB  [5].
Figure1. Phylogenetical trees illustrating evolutionary relationships between M. tuberculosis spoligotypes in Malaysia. (A) Spoligoforest tree drawn using the SpolTools software (available through http://www. emi. unsw. edu. au/spolTools [14, 15] ), and shown as Fruchterman Reingold tree. The Figure was drawn on all patterns including orphan patterns (n5220). In both cases, each spoligotype pattern from the study is represented by a node with area size being proportional to the total number of isolates with that specific pattern. Changes (loss of spacers) are represented by directed edges between nodes, with the arrowheads pointing to descendant spoligotypes. In this representation, the heuristic used selects a single inbound edge with a maximum weight using a Zipf model. Solid black lines link patterns that are very similar, i.
doi: 10. Figure 1. 101 Figure1. Phylogenetical trees illustrating evolutionary relationships between M. tuberculosis spoligotypes in Malaysia. (A) Spoligoforest tree drawn using the SpolTools software (available through http://www. emi. unsw. edu. au/spolTools [14, 15] ), and shown as Fruchterman Reingold tree. The Figure was drawn on all patterns including orphan patterns (n5220). In both cases, each spoligotype pattern from the study is represented by a node with area size being proportional to the total number of isolates with that specific pattern. Changes (loss of spacers) are represented by directed edges between nodes, with the arrowheads pointing to descendant spoligotypes. In this representation, the heuristic used selects a single inbound edge with a maximum weight using a Zipf model.
A total of 18/55 SITs containing 161 isolates were clustered within this study (2 to 56 isolates per cluster) while 37/55 SITs containing 37 strains were unique (for total unique strains, one should add to this number the 22 orphan strains, which brings the number of unclustered isolates in this study to 59/220 or 26. 82%, and clustered isolates to 161/220 or 73. 18%). Note that SITs followed by an asterisk indicates "newly created'' SITs due to 2 or more strains belonging to an identical new pattern within this study or after a match with an orphan in the database; SIT designations followed by number of strains: 3993* this study n51, MYS n51; 3994* this study n52, MYS n51; 3994* this study n52; 3995* this study n52; 3996* this study n52; 3997* this study n52; 3998* this study n52. b.
2007; Mikkelsen et al. 2007 ). We further collected data for mRNA expression levels in five biological replicates of these progenitor rat strains. Subsequently, we collected genome-wide H3K4me3, H3K27me3, and RNA-seq profiles across a panel of 30 BXH/HXB RI rats. We then combined these comprehensive chromatin and expression profiles and analyzed them computationally (see Supplemental Fig. 1 ). We defined regions of interest based on peak calling algorithms or known gene models and quantified histone modification levels as normalized ChIP-seq read counts. The high quality of all ChIP-seq data is in line with recently published standard measurements  (Landt et al. 2012 ; Supplemental Table  1 ), and typical histone occupancy profiles  (Barski et al. 2007)  were observed for expressed and nonexpressed genes (Supplemental Fig.
2010). For each motif of the TRANSFAC database  (Wingender et al. 1996 ; http://www. gene-regulation. com), we tested whether high scores for differential TF binding occurred more often in cis-regulated regions compared with the remaining regions. In heart tissue we found 14 motifs of TF, whose binding affinity to cis-regulated regions was significantly altered (FDR < 0. 1; Supplemental Table  4 ). Among these factors was T-box 5 (TBX5), which is involved in heart development (Fig. 3C, D ). In liver there were eight motifs of differentially occupied TF. One of the TFs was hepatocyte nuclear factor 4, alpha (HNF4A), which is essential for liver development and maintenance of liver-specific gene expression. These results suggest a novel role of the identified transcription factors to directly or indirectly recruit H3K4 methyltransferase or demethylase complexes to the promoters of the differentially modified genes in heart and liver.
3E ). To validate our findings in vivo, we performed ChIP-seq of SHR. BN-chr3 congenic rats, in which the region of the QTL hotspot at chromosome 3 was transferred from the BN donor strain into the genetic background of the SHR recipient strain. We then compared differential data from congenic rats and QTLs from BXH/HXB RI rats and detected a significant overlap of differentially regulated histone marks in the congenic rats with those found to be coregulated by the QTL hotspot at chromosome 3 (Fig. 3F ). Using the congenic rats, we confirmed 768 (85%) of the trans-linked H3K4me3 marks, and interestingly, higher histone modification levels were always associated with the BN allele, indicating a common regulation of these marks. Genes corresponding to trans-linked promoter regions (<4 kb from TSS) were enriched for GO terms chromatin organization (P = 3.
mRNA preparation and RNA-seq data generation We performed 80 RNA-seq experiments to determine gene expression and alternative promoter usage of 30 RI strains and five biological replicates of their progenitor strains. Total RNA was extracted from heart and liver tissue and RNA quality was assessed (2100 Bioanalyzer, Agilent) to assure a RIN score >9. 5 for all samples. After heat fragmentation, rRNA-depleted poly(A) + mRNA was used to generate random-primed double-stranded cDNA. Adapters were ligated using the SPRI bead system and the final library was amplified with 15 PCR cycles. We then sequenced one sample per lane on a HiSeq 2000 instrument from Illumina with TruSeq 23 100 bp PE chemistry. Reads were mapped to the BN reference genome RGSC 3. 4 using TopHat v1. 2. 0  (Trapnell et al. 2009 ).
This approach allows alignment of reads across known and predicted splice junctions. We supplied known splice junctions annotated in the Ensembl reference database  (Flicek et al. 2012)  to the pipeline and also enabled de novo splice junction detection. Gene expression levels were estimated using read counts within gene bodies. We counted all alignments that were overlapping at least one exon of the gene and normalized the counts to the length of the union of all exonic regions of the gene. Gene expression values were normalized across samples using a quantilebased scaling method  (Schulte et al. 2010). Differential analysis and QTL mapping of gene expression levels were performed using the methods described in the sections ''Differential testing of histone marks in progenitor strains'' and ''Differential testing of histone marks in RI strains and genetic mapping.
Quantification and data normalization For each sample, we counted aligned reads that fall in a region with an identified peak. In order to be able to compare these values across samples, we have applied quantile normalization  (Kasowski et al. 2010). We selected only reproducible regions that had a minimum average coverage of 13 in all three replicates of one of the parental strains. For H3K4me3 and H3K27me3 we also require regions to have more than 13 average coverage in at least 15% of the recombinant inbred strains.
Quality assessment of ChIP-seq data We assessed the quality of our ChIP-seq data using metrics developed by the ENCODE Consortium  (Landt et al. 2012) , which is summarized in Supplemental Table  1. For H3K4me3 we compared the MACS peak regions to genome annotations from Ensembl release 56 (Supplemental Fig. 10 ). In addition we generated read coverage plots (Supplemental Fig. 2 ) similar to  Barski et al. (2007). For H3K4me3 and H3K4me1 we plotted the average read coverage 2000 bp upstream of and downstream from the TSS for Ensemblannotated expressed and nonexpressed genes. For H3K27me3 and H4K20me1 we plotted the average read coverage 2000 bp upstream of and downstream from the gene combined with the average coverage per base in 10 bins across the length of the gene body for expressed and nonexpressed genes.
Differential transcription factor binding In order to determine whether the binding sites of particular transcription factors were preferentially altered by SNPs in cisregulated H3K4me3 regions, we first obtained 50 bp of flanking sequences for all SNPs located in all H3K4me3 regions subject to QTL mapping in each respective tissue. For each PWM from the TRANSFAC database set  (Wingender et al. 1996 ; http://www. generegulation. com) of nonredundant PWMs, we scored the difference of binding affinities of the reference and alternative alleles using the sTRAP method  (Manke et al. 2010). We then classified the SNPs into two classes: SNPs in cis-regulated regions (FDR < 0. 1) and SNPs in other regions. We tested for each TF whether the absolute value of the log ratio of binding P-values was larger in the cis-regulated group compared with the other group using a one-sided Wilcoxon-Mann-Whitney test.
Integration with gene expression data For integration of ChIP-seq and RNA-seq data sets, we used a genebased analysis of Ensembl-annotated genes, as this is the more conservative approach, to minimize false-positive findings in our data. Since H3K4me3 is located in promoter regions, we quantified ChIP-seq read coverage in regions 2-kb upstream of and downstream from annotated TSS (Supplemental Fig. 2 ; Supplemental Table  7 ). For H3K27me3 we quantified reads in the entire gene body of Ensembl annotated genes. In order to assess the relation between genetic variants, histone modification levels, and gene expression, we followed a two-step procedure. Each gene is represented by three random variables, which correspond to each of the three traits. First we performed a QTL analysis for each of the three traits separately.
Genotypes are encoded as discrete variables with levels 0 and 1. The log-transformed histone and expression values are modeled as continuous variables. The like-lihood function of the model can be factorized as L(G, P) = P i f[x i |pa(i)]. All models were fitted using maximum likelihood estimators for the parameters. We have used the statistical software R in order to implement the graphical model analysis. The source code of the R package is available in the Supplemental Material and on our website http://www. molgen. mpg. de/;heinig/ histoneQTL/. In order to select the best model, we compute the Akaike information criterion [AIC = À2 log L(G,P) + 2 p, where p is the number of parameters of the distributions in P] and select the model with the smallest AIC. The robustness of the results is assessed by drawing 100 bootstrap samples and recording how many times a given model had the minimal AIC.
Quantitative real-time PCR The findings of ChIP-seq data were confirmed by ChIP followed by regular qPCR on the ABI 7900HT detection system (Applied Biosystems). PCR primers were designed to amplify designated genomic regions using Primer Express software (Applied Biosystems). qPCR assays were carried out in 384-well plates with a final volume of 20 mL each for 40 cycles. We used Power SYBR Green PCR Master Mix (Applied Biosystems) with 1/10 diluted ChIPed DNA or unenriched input DNA as template. Enrichment ratios were calculated according to the 2 ÀnnCt method with endogenous controls (Gapdh, Myod1). Primer and probe sequences can be found in Supplemental Table  8.