retroviruses , transposons , and other mobile dna elements directly integrate their dna into the chromosomes of host cells.1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 distributions of newly integrated dna elements can be characterized using next - generation sequencing , as is described in the companion paper in this issue of moleculartherapy : methods & clinical development and in many previous studies ( e.g. , bushman , craig et  al . , maldarelli et  al . , cohn et  al . , wu et  al . , contemporary methods take advantage of illumina paired - end sequencing to report the location of newly integrated dna.13 , 14 , 15 , 16 multiple reports have described methods for quantifying and analyzing such data , but optimal methods for statistical analysis , data reduction , and data visualization are the topics of ongoing development.13 , 14 , 16 , 17 , 18 , 19 here , we describe a suite of tools for integration site analysis and visualization that are useful for characterizing samples from human gene therapy and other applications . in the case of gene modification in circulating blood cells , it is possible to sample cell populations from blood longitudinally and sequence sites of integration of the gene - correcting vector.20 , 21 , 22 , 23 , 24 , 25 , 26 , 27 an important question centers on how best to quantify the numbers and types of gene - modified cell clones contributing blood cells to the periphery . for example , adverse events have been reported where expanded cell clones in blood became frank leukemia,28 , 29 , 30 , 31 and this can be tracked using quantitative integration site data . complicating the analysis , simply counting the number of integration site sequence reads does not accurately report clonal abundance because of distortions resulting from pcr steps in the integration site recovery procedure.32 , 33 we have previously described a method for abundance estimation based on paired - end sequencing of pcr products containing integration site sequences that allows for accurate quantification of gene - modified cells . dna is sheared using sonication , dna linkers are ligated to free dna ends , and then samples are amplified using primers complementary to the integrated vectors and ligated linkers . genomic sequence information is acquired from both the linker end and the integrated vector end . for the case of an expanded clone , many different dna breaks and sites of linker ligation are associated with the unique integration site from the expanded clone . this allows for the estimation of abundance using the number of different linker ligation sites as a surrogate for the number of cells sampled . we have published statistical tools for analysis of such data and applied these tools to track several gene therapy trials.20 , 21 , 25 , 26 , 34 , 35 , 36 , 37 , 38 , 39 other groups have also used related methods.6 , 8 , 24 , 40 , 41 , 42 , 43 , 44 here , we present tools for integration site sequence analysis and the quantification of clonal abundance , and describe applications in human gene therapy . these methods can also be used for tracking latently infected cells in hiv - positive subjects and monitoring experiments using insertional mutagens , and in mechanistic studies of dna integration . we describe a heatmap format for the analysis of relationships among integration site distributions , genomic features , and sites of epigenetic modification . these analyses allow users to carry out numerous custom statistical comparisons with annotations , random distributions , and other datasets by simply pointing and clicking . we also present a series of analytical tools for use with patient samples to characterize integration site population structure and possible adverse events . results are packaged into reproducible reports ( html or pdf file format ) , allowing for version tracking of the code , datasets used , and external datasets queried . using these tools , we describe examples of tracking a subject from the first gene therapy trial to treat severe combined immunodeficiency - x1 ( scid - x1 ) deficiency . the data demonstrate durable reconstitution accompanied by a clonal expansion of cells harboring an integrated vector near the cancer - associated gene ccnd2 . the first steps involve the generation of a sequence library and sequence data acquisition , genomic alignment , analysis of viral integration sites in repeated sequences , and quality controls ( described in the accompanying paper ) . after sequence quality filtering , trimming of dna sequences added during library construction , and sequence alignment , unique sites and integration sites in repeated sequences (  multihits  ) are saved in an intermediary binary format ( rdata file format ) . all sites from each run of the sequencing instrument are uploaded into a mysql database ( intsitedb ) for storage and for use in downstream analysis using a utility script ( intsiteuploader ) . alternatively , the inspiired pipeline also supports use of a sqlite database , which is provided with the software . the intsitedb stores genomic locations of integration sites together with pcr break points and their counts , which is used for estimation of abundance . we use intsiteretreiver , a key component of the inspiired software package , to retrieve unique sites and multihits for chosen samples from the intsitedb for analysis . integration sites are annotated with the hiannotator r software package ( http://bioconductor.org/packages/release/bioc/html/hiannotator.html ) , which makes use of genomic features compiled by the ucsc genome bioinformatics group . for studies of gene therapy subjects , patient metadata and specimen information are stored in an accompanying gene therapy specimen management database that includes anonymized patient identifiers , cell types analyzed , and time point data . the separation of the pipeline into several output products ( patient report and heatmaps ) and databases ( integration site , patient metadata , and annotation ) provides flexibility in development and use . we use heatmaps to summarize the relationships between integration site distributions and genomic annotations ( figures 2 , 3 , and 4 ) . these maps were introduced in berry et  al . , which presents more background and examples of their uses . the heatmaps summarize information on integration site datasets in columns and different genomic features in rows . for each comparison , integration site distributions are compared with distributions of randomly selected sites in the human genome.4 , 35 early integration site recovery methods involved use of cleavage by restriction enzymes , which are unevenly distributed in the human genome.32 , 33 for this reason , random controls were matched based on proximity to restriction enzyme cleavage sites . the inspiired pipeline uses random cleavage by sonication , so purely random control sites are generated in  silico and used in the analysis described here . in the heatmaps , colored tiles indicate the intensity and direction of any departures from the distributions of random controls for each genomic feature in each integration site dataset . the coincidence of genomic feature  j  with each integration site and random control site is measured . the nonparametric method of estimating receiver operating characteristic ( roc ) curve areas and their covariance structure of delong et  al . is used . each integration site is compared in a pairwise fashion with random control sites , and a number is assigned indicating the relative rank of the integration site : 1 if the measurement of j is higher at the integration site than at a random control site , 0 if the measurement of j is lower at the integration site than at a random control site , and 0.5 if the measurement of j is equal for the two sites . all such values are calculated for a dataset of integration sites and averaged to obtain the overall roc area for the feature measured ( https://github.com/bushmanlab/hotrocs ) . this is equivalent to comparing the ranks of the sites with those of the controls . in older datasets with integration sites recovered by cleavage with restriction enzymes , matched random controls based on proximity to restriction enzyme cleavage sites were used . in that setting roc curve areas were based on comparing each site only with its matched controls . an roc area between 0 and 0.5 indicates the genomic feature occurs less frequently at or near integration sites than at or near random sites in the genome and is therefore disfavored . an roc area between 0.5  and 1 indicates the genomic feature is enriched at integration sites . an roc area of exactly 0.5 indicates that integration sites in the dataset are neither enriched nor depleted with respect to the feature of interest . the roc area is converted to a color tile according to the colorimetric scale shown at the bottom of the heatmap . in figure  2 , positive associations ( enrichment compared with random ) are shown as increasingly intense shades of blue , negative associations ( depletion compared with random ) as increasing intense shades of yellow , and no difference from random as black . each tile represents a comparison of integration sites with the randomly sampled controls for one genomic feature ( row ) in one experimental dataset ( column ) . note that we do not present the magnitude of effect in terms of the original units of measurement . we simply ask whether the average integration site has a higher rank for a given type of feature than its matched random control sites . the color indicates the average quantile of each integration site relative to its random controls . this removes skewing effects contributed by non - normal distributions of the data and also reduces the effect of data points with extreme values for a feature . statistical tests are carried out to determine whether the roc areas calculated are significantly different from one another or from 0.5 ( indistinguishable from random controls ; methods are further explained in supplemental materials and methods).35 , 46 the heatmap shown in figure  2 compares the integration site distribution at two time points from a gene - corrected patient ( patient 1  [ p1 ] ) with epigenetic marks mapped in cd133 progenitor cells . patient 1 ( p1 ) was treated with an early gammaretroviral vector , used to deliver the missing il2rg gene to treat scid - x1.47 , 48 the two samples were isolated from peripheral blood mononuclear cells ( pbmcs ) taken at 177 and 189.5  months after gene therapy . for comparison , another sample using a lentiviral - vector - infected human - derived hap-1 cell line has been included to illustrate differences with the lentiviral integration pattern . the distributions of integration sites datasets ( table s1 ) were compared with the distributions of 10 different epigenetic marks or bound dna binding proteins . each of these was mapped by  chromatin immunoprecipitation sequencing ( chip - seq ) , in which each protein was covalently cross - linked to dna , and bound dna fragments were recovered by immunoprecipitation , sequenced , and then mapped to the human genome to identify relative density . densities of mapped chip - seq annotations were compared with distributions of integration sites within 10 kb windows , and the collection of values was used to generate roc areas . for the gene therapy specimens made by infection of stem cells with a gammaretroviral vector ( patient 1 [ p1 ] ) , the distribution mostly favored marks associated with active transcription ( h3k9me1 , h3k4me1 , h4k20me1 , and h3k27me1 ) . integration was disfavored near marks associated with repressive chromatin ( h3k27me3 and h3k9me3 ) . however , the gammaretroviruses favor integration near transcription start sites , and h2az and h3k4me3 were positively associated , and rna polymerase ii ( pol ii ) more strongly than for lentiviruses . for the lentiviral infection in hap1 cells , integration is favored near sites of covalent modifications of histones associated with active transcription , including h3k9me1 , h3k4me1 , h4k20me1 , h3k27me1 , and h3k6me3 . lentiviral integration is favored within transcription units,1 , 3 , 4 where the h3k36me3 mark is found , so integration was favored near this mark for lentiviral infection , but not gammaretroviral infection . marks found near transcription start sites ( h3k4me3 and h2az ) were either slightly disfavored or neither favored nor disfavored . these patterns parallel those seen previously for lentiviral vectors in diverse cell types.4 , 8 , 11 , 52 , 53 an added feature of these heatmaps is that they have been engineered to allow interactive statistical tests ( figure  3 ; interactive heatmaps are available in supplemental information and data s1 ) . heatmaps are generated as scalable vector graphics ( svg ) , which can be opened in an internet browser . users can click on a row or column , and statistical results appear on the heatmap tiles documenting whether results in other rows or columns differ from the query . users can also click on a button to the right of the maps to allow comparison of all tiles with the random control . some examples are shown in figure  3 , illustrating comparisons among random ( figure  3a ) , the leftmost hap1 dataset ( figure  3b ) , or the pol ii chip - seq distribution ( figure  3c ) . figure  4 presents another form of the heatmap that queries the results of multiple additional features , including mapped dnase i cleavage site ( which reports dna accessibility ) , cpg islands ( important in gene regulation ) , guanine / cytosine ( gc ) percentage , gene counts as documented in the refseq dataset , and proximity to gene boundaries . for those features that are mapped in intervals ( gc percentage over 1  mb , 100 kb , 10 kb , and so on ) , it is often unknown a priori what width is the most relevant to the biological question at hand . thus , for these features , results for a number of different interval sizes are shown . all three datasets are compared over their four replicates against these features ( figure  4 ) . high densities of dnase i hypersensitive sites and high densities of cpg islands are associated with favored integration for both lentiviral and gammaretroviral vectors ( red coloration ) . high gc content is also favored , likely because high gc content is characteristic of gene - rich regions , although for lentiviruses , the preference switches to local high adenine / thymine ( at ) , possibly associated with binding of ledgf / p75 , the tethering cofactor , or wrapping of integration target site dna on nucleosomes . paralleling the favored high gc content , direct measures of gene richness integration is disfavored relative to regions with long gene widths or long intergenic distances , because these are indicative of gene - sparse regions . the gammaretroviral vector sites are disfavored relative to long gene boundary distances because they favor integration near transcription start sites . integration is favored within genes ( as annotated by the refseq dataset ) for lentiviral vectors , but only weakly favored for gammaretroviral vectors . thus , numerous relationships between integration site datasets and genomic features can be explored statistically using these interactive heatmaps . a question of interest in many therapeutic applications centers on whether integration sites accumulate near the transcription start sites of cancer - related genes . a complication is that there are many ways of defining cancer - associated genes , and most such genes are important only in specific types of human cancers . for annotating gene therapy results , we have thus generated multiple lists of cancer - associated genes that can be queried as appropriate for integration site analysis ( http://www.bushmanlab.org/links/genelists ) . in one approach , we created a maximally comprehensive list ( allonco ) for use in first - pass screening based on the idea that we can not predict what cancer - associated genes are most important in the novel clinical setting of human gene modification . the list incorporates known human cancer genes and human homologs of cancer genes in model organisms , and so includes 2,125 total genes , or roughly 8.5% of all human genes ( assuming 25,000 total ) . comparison with oncogene annotation is summarized using the heatmap format ( figure  4 , bottom row ) , which scores the frequency of integration sites within 100 kb of cancer gene transcription start sites in integration sites versus random sites . for this we have developed a standardized patient report template that rests on top of the inspiired pipeline . use of a reproducible report format allows tracking of datasets used in each study and version control of code ( which are specified by dates ) . the report software takes in integration site and break point positional information , annotates the sites using hiannotator , and outputs targeted analyses of integration site distributions . an example of a patient report is provided in data s2 , showing two recent time points monitored for patient 1 ( p1 ) treated for scid - x1 . earlier time points were analyzed by 454 roche pyrosequencing and were previously reported . some excerpts from the report the software generates a summary table ( figure  5a ) that reports the patient , time point , cell type , patient metadata , and summary statistics for each sample . among these are the total numbers of reads , the number of cells inferred to have been sampled ( sum of break points captured ) , and the number of unique integration sites after dereplication . the minimum population size is inferred from a chao1 analysis with jackknife correction , which takes advantage of the four replicate analyses typically run for each sample . skewing in proportional abundance is calculated using the gini index , where 0 indicates a perfectly even distribution of integration sites over the cells sampled and increases up to 1 with increasing oligoclonality . diversity is calculated using the shannon index , which summarizes both the number of different unique integration sites and the evenness of distribution of cells sampled ( sonicabundance ) among integration sites . here , we introduce a new metric , called the uc50 ( unique cell progenitors contributing the most to the expanded 50% of progeny cell clones ) . to generate the number , progenitor cell clones ( reported as unique integration sites ) are first ranked by the relative abundance of progeny cells using sonicabundance ( reported by linker ligation site data ) . the uc50 then reports the number of unique clones ( integration sites ) responsible for making up the top 50% of all cells sampled . thus , if a single clone comprises more than 50% of the sample , the uc50 value will be 1 . in contrast , for efficient lentiviral infections of cells in short - term tissue culture , the uc50 values can be in the thousands ( data not shown ) . finally , where available , the vector copy number per cell ( vcn ) , determined separately by qpcr , is added to allow assessment of the efficiency of gene marking in the cell population . the relative abundance of integration sites in or near specific genes is summarized in several ways ( data s2 ) . these include two types of stacked bar graphs ( an example is shown in figure  5b ; both are displayed in data s2 ) . bar graphs either display the number of cells observed ( sonic breaks ) that are associated with integration sites in the samples ( scaled to the most abundant sample ) or the proportional abundance of integration sites in each sample , where every sample is scaled to 100% ( figure  5b ) . each high - abundance integration site is named by the closest gene , which is color coded in the key for each graph . genes are annotated by whether the site is within a transcription unit ( * ) , is within 50 kb of a cancer - related gene (  ) , or is associated with a gene strongly associated with human lymphoma ( ! ) . this last list consists of 38 human genes commonly involved in lymphoid cancers and includes most genes previously implicated in adverse events in human stem cell gene therapy , including lmo2 , mds / evi1 , and ccnd2.29 , 30 , 31 heatmaps provide another type of visualization ( data s2 ) . these have the advantage over stacked bar graphs in that each integration site above a chosen abundance is given equal space , whereas in stacked bar graphs , rarer sites can be shown as thin bands that may be difficult to visualize . a third visualization is provided by line graphs ( data s2 ) , which highlight the behavior of the most abundant clones over time . another set of figures queries integration sites near genes of concern for adverse events in human gene therapy , including lmo2 , ccnd2 , hmga2 , and mecom.29 , 37 in these visualizations , the distance from the transcription start site is shown on the x axis , proportional abundance is shown on the y axis , and the time point is color coded . by this means it can be seen that an integration site near ccnd2 achieves > 10% abundance near the ccnd2 transcription start site ( figure  5c ) . handling integration sites that map to multiple locations in the human genome presents a particular challenge . it could be that an integration site authentically resides in a repeated sequence that is also in the 5 region of a cancer - associated gene and potentially marking an adverse event . to accommodate this , reports include an account of the sonicabundance of cells with integration sites in repeated sequences , allowing tracking of possible sites of concern in multihits . the reports end with word bubbles ( figure  5d ) , which provide a visual summary of genes near integration sites in expanded clones . names of genes that are nearest to the integration sites are used to construct the word bubble . word size is scaled by the sonicabundance measure for each integration site , and each gene name is marked with the same integration flags ( * ,  , and ! ) used in earlier plots . by this means , figure  5 shows excerpts from a reproducible report summarizing monitoring of patient 1 ( p1 ) from the first trial to treat scid - x1 . results are summarized for pbmcs from two time points , 177 and 189.5  months after gene therapy . half a million to a million reads were collected for each sample , allowing investigation of 13,000 cells associated with about 1,000 integration sites . uc50 values for the two time points are 10 and 8 , indicating the presence of expanded cell clones . in early studies based on 454/roche pyrosequencing , the subject was found to have an expanded clone with an integration site near ccnd2 ( 6% of all reads ) , a gene for which a nearby integration event was associated with an adverse event in another scid - x1 gene - corrected patient.25 , 29 thus , it was of interest to monitor the behavior of the clone in this patient over time . analyses using illumina paired - end sequencing are summarized in figure  5 , which shows that the ccnd2 clone has slightly expanded in abundance ( nonparametric comparison of replicate medians yields p  = 0.029 when compared by relative abundances and p  = 0.057 when compared by absolute abundances judged by sonicabundance ) . here , we describe a collection of tools for the analysis and visualization of integration site distributions . integration sites are mapped to the human genome , the abundances of their host cells tabulated using the sonicabundance method , and results stored in a database , allowing flexible downstream analysis . integration site distributions can then be compared with genomic features and sites of epigenetic modification as roc areas . these data are summarized as interactive heatmaps , allowing comparison with random distributions , other integration site datasets , or genomic annotation by simply clicking on a row or column , which outputs comprehensive statistical tests . for applications to human gene therapy , interest often focuses on longitudinal behavior of expanded clones and proximity of integration sites to cancer - associated genes . a standardized report format was developed , allowing interactive comparison among patient datasets and querying multiple aspects of longitudinal behavior . for this , we introduce the uc50 metric , which is generated by ranking progenitors ( integration sites ) from most to fewest daughter cells produced ( linker ligation sites ) and counting the number of progenitors contributing to the top 50% of the distribution . thus , clonal expansion yields low uc50 numbers and highly polyclonal samples high uc50 numbers . these tools were used to query recent clonal behavior in a patient from the first scid - x1 gene therapy trial.29 , 47 , 48 the patient studied has an expanded clone with an integration site near the proto - oncogene ccnd2 . the analysis of integration sites from month 177 to month 189.5 post - treatment revealed stability of this clone , with possible slow expansion . this analysis illustrates how the tools described here can be applied to monitor outcomes in gene therapy . as in cavazzana - calvo et  al . and hacein - bey - abina et  al . , patient 1 ( p1 ) fulfilled the eligibility requirements for first ex  vivo c gene therapy trial ( 19992002 ) at age 11  months . p1 was diagnosed with scid - x1 based on his blood lymphocyte phenotype , revealing a tail - less c receptor expressed at the membrane ( r289 x ) . marrow was harvested and subjected to cd34 cell separation , obtaining 9.8   10 cd34 cells per kilogram of body weight . harvested cells were then exposed to mfg c vector - containing supernatant daily  for 3  days . p1 was then infused with the treated cd34 cells ( 19   10 cells / kg ) without prior chemoablation . as explained in the companion paper , integration sites are identified by sequencing the ltr - host junctions from genomic dna after linker - mediated pcr amplification . genomic dna is randomly sheared by ultrasonication , after which linkers are ligated to the repaired dna for amplification . nested pcr is used to amplify the ltr - host dna junctions by priming from the viral ltr and the linkers , appending the sequences needed for sequencing . samples are sequenced using the illumina paired - end platform , and the output sequencing files are processed by intsitecaller to yield integration site positions on a host draft genome . integration site data and chip - seq data were mapped onto the hg18 genome draft , to match the original draft genome used for analysis of the chip - seq data . as in berry et  al . , receiver operating characteristic ( roc ) areas are used to compare integration sites with random control sites . inspiired is distributed online as a downloadable virtual machine  executable on the windows , mac , and linux operating systems ,  as  well as a github source code repository supported by a conda  software environment ( see https://github.com/bushmanlab/inspiired , which also includes detailed instructions for use and test datasets ) .
analysis of sites of newly integrated dna in cellular genomes is important to several fields , but methods for analyzing and visualizing these datasets are still under development . here , we describe tools for data analysis and visualization that take as input integration site data from our inspiired pipeline . 
 paired - end sequencing allows inference of the numbers of transduced cells as well as the distributions of integration sites in target genomes . 
 we present interactive heatmaps that allow comparison of distributions of integration sites to genomic features and that support numerous user - defined statistical tests . to summarize integration site data from human gene therapy samples , we developed a reproducible report format that catalogs sample population structure , longitudinal dynamics , and integration frequency near cancer - associated genes . 
 we also introduce a novel summary statistic , the uc50 ( unique cell progenitors contributing the most expanded 50% of progeny cell clones ) , which provides a single number summarizing possible clonal expansion . using these tools , 
 we characterize ongoing longitudinal characterization of a patient from the first trial to treat severe combined immunodeficiency - x1 ( scid - x1 ) , showing successful reconstitution for 15 years accompanied by persistence of a cell clone with an  integration site near the cancer - associated gene ccnd2 . 
 software is available at https://github.com/bushmanlab/inspiired .