BACKGROUND
the advent of next generation sequencing has enabled the interrogation of the human genome and transcriptome with base pair resolution.
this new window into the human genome has expanded our understanding of important somatic genome variants such as single nucleotide polymorphisms , chromosomal rearrangements , and exogenous dna integrations  <cit> .
these variants have primarily been identified using paired-end sequencing, a process that relies on breaking the genomic dna into many small fragments and sequencing each fragment of dna from the ends inward toward the center.
the resulting pairs of sequencing reads share a unique relationship as they originated from the same randomly sheared dna fragment.
by leveraging the unique relationship between the paired-end reads, chromosomal rearrangements and integrations of exogenous dna can be identified and characterized  <cit> .
for example, a rearrangement between two chromosomes would result in at least one chimeric chromosome .
by mapping the sequencing reads to a reference genome, a subset of the paired-end reads will support the chromosomal rearrangement by spanning the break point between the two chromosomes with one of the paired-end reads mapping to chromosome a, while the respective other paired-end read maps to chromosome b .
in addition to identifying rearrangements, this technique has also been applied to identify integrations into the human genome of: transposable elements , viral genomes  <cit> , and bacterial dna  <cit> .
if the integration site can be rebuilt in silico, the structure and sequence of the integration site can be used to determine a potential mechanism of integration  <cit> , elucidate potential functional implications of novel integrations, and lay the foundation for future experiments.fig.
<dig> a brief schematic illustrates using next generation sequencing to identify structural variation.
a recombination event is illustrated between two chromosomes .
in order to identify structural variation, paired-end reads are identified that map to different chromosomes indicating a fusion of the two different chromosomes has occurred
previously, we identified paired-end reads supporting the integration of pseudomonas-like rrna genes into the 5′-utr of genes in human stomach adenocarcinoma  genomes using rna-seq data from the cancer genome atlas   <cit> .
these paired-end reads have one read mapping exclusively to the 16s or 23s rrna genes of a pseudomonas-like bacteria  <cit> , while the paired read maps uniquely to the 5′-utr of the human ceacam <dig>  ceacam <dig>  cd <dig>  or tmsb <dig> genes.
as such, these read pairs support bacterial dna integrations by spanning the junction of the bacterial and human dna.
to identify these regions in the human genome, >4× human coverage in a single sequencing run was required.
despite this level of coverage, single reads that traverse the integration site were not identified, likely owing to the length of the sequencing reads and limits of alignment algorithms.
such reads would enable the assembly of the integration with the base pair resolution needed to determine the potential mechanism of integration  <cit> .
instead, here, we have used the paired-end reads to model the most likely structure of the integration of pseudomonas-like rrna gene fragments into ceacam <dig>  ceacam <dig>  cd <dig>  and tmsb <dig>
RESULTS
establishing the boundaries of the bacterial rrna integrations
in all cases, reads were only recovered spanning one side of the junction between bacterial and human dna, likely because the integration is close to the transcriptional start site.
despite numerous attempts, the integrations could not be assembled using the original reads supporting an integration event or with additional reads located near the integrations that were identified using either alignment- or blast-based algorithms.
neither split reads, those that are part bacterial and part human on the breakpoint of the integrations, nor soft-clipped reads could be identified in the original alignments or in alignments to a custom reference with both the human genome and the pseudomonas 16s & 23s rrna genes.
further blast-based examination of the unmapped read in read pairs that had only one read aligned to the human genome in the region flanking the integration or the pseudomonas rrna gene reference near the integration also failed to identify split reads.
to examine this further, a bacterial-human dna integration was constructed with the bacterial dna directly abutting the human dna.
a mock dataset was created of all  <dig> possible combinations of 100-bp paired end reads spanning the integration breakpoint in this artificial sequence mimicking a bacterial-human dna integration.
the first read generated was entirely bacterial and ended at the integration breakpoint.
each of the  <dig> subsequent reads in the mock dataset shifted by 1-bp, such that the dataset included a mock read for every position across the integration beginning with an entirely bacterial read and ending with an entirely human read.
the second read in the pair was held constant and corresponded to a sequence 225-bp downstream of the break point.
lgtseek identified only  <dig>  reads that cover the breakpoint, none of which were soft-clipped as the differences with respect to mapping were similar to those arising from sequencing errors.
therefore, we conclude that lgtseek, and more specifically the version of bwa used in lgtseek, is unable to identify reads that span the junction between bacterial and human dna in this data set.
given that the bacterial dna integrations could not be assembled, the focus shifted to estimating the location of the bacterial rrna gene fragment integrations into the human genome by examining the structure of the human transcript and the reads supporting the bacterial dna integrations.
the integration breakpoint must be downstream of the transcriptional start site  of each human gene for three reasons.
first, the integrations were identified in an rna-seq data set derived from transcripts so they must be within the transcript boundaries.
second, examination of the expression of these genes across all participants from the stad and breast cancer  data sets from tcga data available for download from the sra between september 18th–20th,  <dig>  <cit>  are consistent with an accurately annotated tss .
third, the introns and other noncoding regions of the human genome do not have low-level sequencing coverage that would suggest the presence of contaminating genomic dna.
therefore, the integration must be downstream of the annotated tss and the tss is the left-most boundary for the possible location of bacterial dna integration, relative to the direction of transcription.
the right-most boundary for the bacterial dna integration break point can be delineated from the position of human reads supporting the bacterial dna integrations.
more specifically, the bacterial dna integration must be upstream of the human reads supporting the bacterial dna integration, relative to the direction of transcription.
therefore, the left most boundary for the site of integration in the human genes, in all cases, should be considered the tss, while the right most boundary position for the integration is the left most position of the consensus sequence for the human reads supporting the bacterial dna integration.
using the library insert-size to refine the location of bacterial dna integration
the library insert-size and its distribution can be used to refine the location of the bacterial rrna gene fragment integration.
in order to model the integrations with the greatest resolution and accuracy, two calculations were used to determine the number of bases needed between the bacterial and human consensus sequences to mirror the library insert-size distribution.
suppose that the distance between the bacterial and human fragments is x, where  <dig> ≤ x ≤  <dig> , and that each integration has n total reads supporting it.
the insert-size, ii , of the ith read pair spanning the junction is calculated based on the known positions of each read within the bacterial and human fragment consensus sequences and the assumed distance, x, between these fragments.
the mean absolute value, ad, of the differences between the median sequencing library insert-size, msl, and the read pairs supporting the integration can be calculated as adx=∑i=1niix−msln, where n is the total number of reads supporting the integration.
the location of the junction between bacterial rrna and the human gene of interest was estimated to be within the range, , where l is the position of the left end of the human fragment consensus sequence within the human gene and xad is the value of x which has the minimum value of ad .fig.
<dig> the location of the bacterial dna integration was refined using the library insert-size.
in this hypothetical case, the paired-end reads support the integration of a fragment of the bacterial 16s rrna gene into a human gene .
using the consensus sequence of the two fragments, the number of bases, xx, needed between the two fragments was titrated so that the insert-sizes of the reads supporting the integration most closely resemble those of the sequencing library.
two calculations were used to determine the optimum distance between the two fragments.
the first calculation identified the number of bases  yielding the minimum average difference  between insert-size of the reads supporting the integration, i
i, and the median sequencing library insert-size, m
sl.
the second calculation determined the number of bases  resulting in the minimum jensen-shannon distance  for the insert-size distributions of the sequencing library, d
sl, and the reads supporting the integration, d
i
the second approach of estimating the distance between the bacterial and human consensus sequences uses a dissimilarity measure called the jensen-shannon distance, jsd, between the insert-size distributions of the sequencing library, dsl, and reads supporting the integration, di .
the jensen-shannon distance between two distributions takes values between  <dig> and  <dig>  with  <dig> implying the identity of the distributions and  <dig> indicating major discordance between the distributions.
as in the average distance case, the location of the junction between bacterial 16s rrna fragment and the human gene of interest was estimated to lie within the range, , where l is the position of the left end of the human fragment consensus sequence within the human gene and xjs is the distance value for the minimum value of jsd .
the majority of the time the ad and jsd calculations were in agreement.
when the two calculations differed, the jsd was preferred since it accounts for the distribution of insert size values for the two populations.
validating the method
to validate this method of modeling the bacterial dna integrations, a bacterial dna integration in a cancer cell line genome was identified, modeled, and experimentally validated.
using the cancer cell line encyclopedia  <cit>  and a similar cancer cell line sequencing project that is publicly available through genentech inc  <cit> , paired-end rna-seq reads were identified that support the likely in vitro integration of a kanamycin resistance gene  near the  <dig>  kbp position in chromosome  <dig> of the kpl- <dig> cell line.
the paired-end reads supporting this integration spanned the junction of bacterial and human dna at the 3' - side of the kanamycin gene, with respect to the direction of transcription.
given that this integration was likely in vitro, a viral promoter is expected to drive the kanamycin gene expression, which would prevent the bacterial sequence focused pipeline from detecting the 5' - junction between the integrated bacterial dna and human chromosome.
in both independently sequenced datasets, the jsd and the ad predicted that there should been 0 bp between the consensus bacterial and human dna .
subsequently, the genomic dna of the kpl- <dig> cell line was acquired from leibniz institute dsmz and the junction between the integrated kanamycin resistance gene and chromosome  <dig> was pcr amplified, cloned, and sequenced.
the sequencing revealed that the jsd method correctly predicted the sequence and approximate position of the integration site, but that the exact position was off by 2 bp, with the actual position of integration being 2 bp from that predicted by the model .fig.
<dig> a bacterial dna integration into the kpl- <dig> cell line genome was modeled using the jsd and ad calculations.
based on the data from the cancer cell line encyclopedia  <cit>  and a similar data set from genentech  <cit> , the jsd and ad methods predicted that the bacterial & human consensus sequences are 0 bp apart .
however, pcr amplification and sequencing revealed an additional 2 bp of sequence between the two consensus sequences as illustrated .
consistent with this, a plot of the sequence coverage of the reads across the junction is relatively constant .
no further insertions or deletions are observed in the underlying reads
kpl <dig> breakpoint
while investigating the kpl <dig> integration, soft-clipped reads at the kpl <dig> bacteria-human dna junction were identified manually.
to better understand why these reads were not identified by the lgtseek pipeline, a mock dataset was created of all  <dig> possible combinations of 100-bp paired end reads spanning the cloned and sequenced integration breakpoint.
the first read generated was entirely bacterial and ended at the integration breakpoint.
each of the  <dig> subsequent reads in the mock dataset shifted by 1-bp, such that the dataset included a mock read for every position across the integration beginning with an entirely bacterial read and ending with an entirely human read.
the second read in the pair was held constant and corresponded to a sequence 190-bp downstream of the break point, which was selected based on the insert-size distribution of the ccle data.
no split reads were identified from running the entire kpl <dig> dataset spiked with this mock dataset through the lgtseek pipeline due to a 24-bp inversion in the bacterial dna near the kpl <dig> breakpoint .
bacterial rrna gene integrations into ceacam <dig> & ceacam6
from the tcga stad data set that was previously analyzed  <cit> , two participants have paired-end reads supporting pseudomonas dna integrating into exon  <dig> of ceacam <dig> .
in participant a, the integration of the pseudomonas 16s rrna gene  fragment is supported by  <dig> paired-end reads .
using the jsd, the calculated distance between the bacterial and human fragment consensus sequences is 26 bp, with a minimum jsd value of  <dig>  , while the ad placed it between 11–17 bp .
a visual inspection of the read insert-size distribution over a variety of values of x  further supports the jsd prediction.
therefore, given that the human fragment consensus sequence starts at position  <dig> in the ceacam <dig> gene, the bacterial dna integration is estimated to occur between positions 48– <dig> in the ceacam <dig> gene .fig.
<dig> the model of the pseudomonas 16s rrna gene fragment integration into ceacam <dig> of participant a is presented.
the structure of the pseudomonas 16s rrna gene fragment, the first exon ceacam <dig>  and the paired reads supporting the integration of the rrna gene fragment into ceacam <dig> are illustrated .
calculations using the jsd  and ad  support that there are 26 bp between the pseudomonas 16s rrna & ceacam <dig> fragments and that the integration is in 48–73 bp of ceacam <dig> .
the ceacam <dig> sequence upstream of the bacterial 16s rrna fragment is for illustrative purposes only .
the insert-size for each paired-end read is color-coded with lighter colors being closer to the median insert-size, red designating those insert-sizes larger than the median, and green for insert-sizes less than the median as further clarified in additional file 7: figure s <dig>  the optimum model for the structure of the integration  is compared to alternative distances  between the two fragments in , where an additional  <dig> ,  <dig> , − <dig> , & - <dig>  median absolute deviations  are placed between the two fragments represented by the gray region.
the actual dna sequence of the gray region is unknown
participant b has  <dig> paired-end reads supporting the pseudomonas 16s rrna gene integration into ceacam <dig> .
both the jsd and the ad calculations are in agreement that the distance between the two fragments is 0 bp , and that the 16s rrna gene fragment has integrated at the 18 bp in ceacam <dig>  however, manually inspecting the insert-size of the reads suggests that the human read farthest upstream is an outlier that does not fit the insert-size distribution of the other reads .
after removing this potential outlier, a model of the integration was calculated .
the jsd and ad support that the distance between the consensus sequences is 7 bp .
in addition, both the jsd and ad decreased , suggesting that this may be a more accurate representation of the integration, particularly given that decreasing the number of reads in the analysis typically increases the jsd as described below.
therefore, the integration of the 16s rrna gene fragment is likely between positions 47–53 bp of ceacam <dig>  for the subsequent comparisons in this manuscript, only this latter model is discussed.fig.
<dig> a model of the pseudomonas 16s rrna gene fragment integration at the 18 bp of ceacam <dig> from participant b is illustrated .
based on the jsd  and ad  calculations, the model has 0 bp between the bacterial and human fragment consensus sequences.
after removing the potential outlier paired-end read, the optimum structure for the integration is illustrated in d. by removing the potential outlier read, both the jsd and ad decreased , supporting that there are 7 bp between the two fragments.
this indicates the pseudomonas 16s rrna likely integrated between positions 45–53 bp of ceacam5
names indicate the participant with letters a through f, their respective rrna gene integration , and the human gene it has integrated into
if the potential outlier paired-end read is excluded, the modeling suggests that the 16s rrna integrations into the ceacam <dig> gene of participants a & b are nearly identical , with the models only differing by a few bp .
in addition, the human reads from participants a & b map to similar positions in ceacam <dig> .
in addition to the 16s rrna gene integration described above, participant a also has  <dig> paired-end reads supporting an integration of 1626–1699 bp from the pseudomonas 23s rrna gene  into the ceacam <dig> gene.
the jsd and ad calculations support that the integration is between positions 49– <dig> of ceacam <dig>  in addition, participant b has  <dig> reads supporting integration of a nearly identical fragment of the 23s rrna gene  integrating into an almost identical region in ceacam <dig> .
therefore, the data support two different pseudomonas dna integrations into approximately the same position in ceacam <dig> for both participants a & b.
similar to ceacam <dig>  participants a & c have similar pseudomonas 16s rrna gene fragment integrations into the related ceacam <dig> .
both participants have nearly identical fragments of the 16s rrna gene  integrated into ceacam <dig> at the same predicted position  .
in addition, these ceacam <dig> integrations are similar to the  position of the 16s & 23s rrna integrations and  the sequence of the 16s rrna integrations into ceacam <dig>
cd <dig> integrations
both participants d and e have fragments of the pseudomonas 16s and 23s rrna integrated in the cd <dig> gene .
both participants have almost identical fragments of the 16s rrna gene  integrated into overlapping positions in the cd <dig> gene.
for participant d, the calculated model supports that the 16s rrna fragment integrated between positions 180–187 bp of cd <dig>  while the fragment of 16s rrna integrated into participant e’s cd <dig> gene between positions 163–205 bp .
the data supporting these two integrations are consistent and suggests that the two participants have similar integrations.
the data also support a second fragment further downstream in the 16s rrna gene  integrating into positions 136–147 bp of cd <dig> in participant d.
there are  <dig> different fragments of the 23s rrna gene fragments with support for integration into the cd <dig> gene of participants d & e. all three integrations have similar locations in the cd <dig> gene  .
of the three integrations, two integrations, one each from participant d & e, are predicted to have integrated into the same region of cd <dig> .
the remaining 23s rrna gene fragment with support for integration into cd <dig> is predicted to have integrated into the 201–225 bp of cd <dig> of participant e .
tmsb <dig> integrations
unlike the bacterial dna integrations into ceacam <dig>  ceacam <dig>  and cd <dig>  the tmsb <dig> gene  has only one participant with the integration of a pseudomonas 16s rrna fragment .
the model estimates that the 16s rrna fragment integrated into the first 13 bp of tmsb <dig> .
calibrating the jsd
as described previously, a calculated jsd value of  <dig> supports that the two populations are identical, while a value of  <dig> indicates major discordance between the distributions.
the average jsd for the models of bacterial dna integration presented here was  <dig> .
therefore, the jsd calculation was investigated further to determine why it was not closer to zero.
to test the effect that the number of reads had on the jsd, the jsd was calculated for  <dig> samples for specific values of k, such that  <dig> ≤ k ≤  <dig> , where k represents the number of randomly sampled read pairs that mapped to the ceacam <dig> transcript from the sequencing library of participant b .
samples with 3– <dig> randomly selected reads had a jsd of ~ <dig> , consistent with the results for modeling the integrations.
as the number of randomly sampled reads increases, the jsd decreases.
the average jsd continues to decrease as the number of reads increases until around  <dig> reads when the jsd starts to plateau, eventually reaching a minimum of  <dig>  at  <dig>  reads, the maximum number tested.fig.
<dig> boxplots illustrate the jsd  calculated for subsets of rna-seq data mapped to the human transcript.
the insert-size distribution of the sequencing library of participant b was compared to the distribution of insert-sizes of specific values of randomly selected reads that mapped to the ceacam <dig> transcript from the sequencing library of participant b .
panel a has a discontinuous x-axis to better illustrate the jsd over the various number of reads tested ; panel b illustrates the same data with a continuous x-axis
guanine rich sequences
while modeling the integrations, it was observed that both the bacterial and human fragment consensus sequences seemed to be particularly rich in guanine.
consistent with this observation, the median %gc for the bacterial and human fragment consensus sequences was  <dig>  % and  <dig>  %, respectively.
the %gc and %ag were compared for the consensus sequences supporting the bacterial dna integrations using the wilcoxon rank-sum test  .
the bacterial and human consensus sequences were significantly different than each other .
both sets of consensus sequences were also significantly different than the complete refseq pseudomonas genomes and the human genomic reference .
the human fragment consensus sequences were also significantly different that the pseudomonas 16s & 23s rrna genes available in the silva database .
considering that the predicted integration sites are typically into the first 150 bp of each human transcript, the %gc of the fragment consensus sequences supporting the integrations was compared with the first 150 bp of each transcript from the human transcriptome reference and found that the %gc for the human fragment consensus sequences was significantly different .
given that the expressed transcripts may have a different gc profile  <cit> , the %gc for the fragment consensus sequences was compared to the expressed transcripts  from the participants with the bacterial dna integrations.
the %gc of the human consensus sequences was significantly different than the expressed transcripts .
next, by stratifying the expressed genes into quartiles based on their rpkm, the %gc for the quartiles of expressed transcripts was compared to consensus sequences.
the human fragment consensus sequences were significantly different than the two most highly expressed quartiles  .
in addition, with a less stringent but acceptable significance level of  <dig> , the %gc for the lowest two quartiles  are likely significantly different than the human fragment consensus sequences .table  <dig> wilcoxon rank-sum test comparison of gc & ag content
* = p <  <dig> ; na not applicable
in addition to gc-content, the enrichment of guanine suggests that these are purine-rich sequences, which have been shown to have a role in regulating transcription .
therefore, the percentage of bases that are purines  for each fragment was calculated using the transcribed strand of each gene and compared to various references using the wrt .
the median %ag for the human and bacterial fragment consensus sequences was  <dig>  &  <dig> , respectively.
the %ag of the bacterial fragment consensus sequences was significantly different than the pseudomonas 16s & 23s rrna gene database, and both of the complete and first 150 bp of the expressed human transcripts in the participants with the bacterial dna integrations .
on the other hand, the %ag of the human fragment consensus sequences were statistically different than the expressed transcripts from the participants with the bacterial dna integrations .
due to the observed higher guanine content of these sequences, the human genes with bacterial dna integrations were investigated for guanine associated motifs.
a database search of published g-quartets in the human genome  <cit>  did not reveal evidence for g-quartets near the integrations.
examination of the first exon and  <dig> kbp upstream of the tss for ceacam <dig>  ceacam <dig>  and cd <dig> did not reveal cpg islands .
on the other hand, all three algorithms support cpg islands in the  <dig> kbp region upstream of the tmsb <dig> tss.
however, only  <dig> of the  <dig> algorithms suggested the cpg island may extend downstream of the tss near the predicted integration into tmsb <dig> .
cpg islands are typically defined as ≥200 bp regions that have more observed than expected cpg dinucleotides  and high gc content .
next, the window size for identifying cpg islands was decreased to refine the location of the cpg island signal in the first exon of tmsb <dig>  the cpg signal was determined to originate from the 3' end of the first exon of tmsb <dig>  while the side of the first exon near the predicted integration site had no cpg signal .
therefore, it is unlikely the bacterial integration into tmsb <dig> disrupted cpg sites.
rrna structure
based on the position of the bacterial dna integrations into 5′-utr of the human genes, it is possible that the bacterial dna integrations altered the expression of the human genes by disrupting a regulatory element in the 5′-utr.
one reason to expect that this might be the case is that the integrations arise from the 16s & 23s rrna genes that contain numerous stem-loop structures that have the potential to alter secondary structure and thus transcription.
therefore, the consensus sequence of the integrated bacterial dna was examined for the presence of stem-loop structures from the bacterial rrna gene by mapping the bacterial fragments onto the known secondary structure .
it was concluded that 11/ <dig> integrations have large stem-loop structures in the middle of the bacterial fragments; while the remaining three fragments of integrated bacterial dna have stem-loop structures near the ends of each fragment.
subsequently, the predicted rna structure of the human transcripts was examined to determine if the bacterial integrations may have disrupted the secondary structure of the human transcripts.
by mapping the predicted position of each integration onto the ceacam <dig>  ceacam <dig>  cd <dig>  and tmsb <dig> predicted mrna transcript structures, it was determined that the integrations are either in, or near, complex secondary structures .
relative position of integration
to explore if the bacterial dna integrations may be occurring in the same relative position within each human gene, the relative position of the predicted site for integration within the first exon was calculated .
while the ceacam <dig> and ceacam <dig> integrations are enriched around the middle of the first exon, the two genes also share homology and have similar structures.
the ceacam genes, cd <dig>  and tmsb <dig> all have different relative position of integrations .table  <dig> relative position of the bacterial dna integrations
discussion
obstacles in assembling the integrations
here, the jensen-shannon distance between the insert-size distributions of the paired-end reads supporting the pseudomonas-like dna integrations relative to the respective sequencing library was used to model bacterial dna integration in ceacam <dig>  ceacam <dig>  cd <dig>  & tmsb <dig>  the jsd values supporting these models are within the expected range of the jsd values given the number of reads supporting the bacterial dna integrations.
ideally, the sequencing reads would be assembled to determine the exact integration site and sequence.
there are likely three compounding factors that prevent the assembly of these bacterial dna integrations.
first, the reads are only 51 bp.
in order to find a read that is on the integration breakpoint, the read on the breakpoint needs to be able to be split into uniquely human and bacterial portions.
as a result of the short read length, a read would have to be exactly half human and half bacterial in order to give a sufficient sequence length of ~25 bp to identify each portion.
this is highly unlikely to occur.
second, there is limited coverage of the integrations.
the limited coverage is likely due to the proximity of the integrations to the beginnings of transcripts where rna-seq coverage is poor.
in addition, it is possible that the low rna-seq coverage is the result of the bacterial dna integrations inducing decreased expression.
lastly, the data suggest a heterogeneous population of transcript variants in a tumor that makes de novo assembly difficult.
in a mock-experiment with an artificial bacterial-human dna integration, lgtseek identified only  <dig>  reads that cover the breakpoint.
these reads align to the bacterial dna, with only 1–3 bp at the 3'-end of the read covering the breakpoint.
lgtseek uses a version of bwa that may not identify these reads.
while a more recent version of bwa identifies soft-clipped reads better, it also works on the assumption that all reads should match the reference and as such has been found to erroneously align bacterial reads to the human reference.
in kpl <dig>  the cloned and sequenced bacteria-human dna junction was used to generate a similar mock dataset that demonstrates that split and soft-clipped reads could not be identified due to a 24-bp inversion.
given that when the bacterial and human dna are flush, only 3 % of read pairs are identified, and that there are only  <dig> read pairs supporting stad bacterial dna integrations on average, we do not expect that soft-clipped reads would be identified even if the bacterial dna is flush to the human dna.
furthermore, such soft clipped reads that do exist would be difficult to distinguish from sequencing error or sequencing polymorphism between the bacterial sequence and the reference genome.
additionally, the average distance between the bacterial and human fragments of dna in the models described here was found to be ~24 bp, which suggests that the two may not be flush, further explaining the difficulty in detecting such reads.
heterogeneous integrations
participants a through e all have data supporting multiple integrations of different fragments of the 16s and 23s rrna pseudomonas gene into the same human gene within each participant.
in addition, the data support multiple participants with similar, if not identical, bacterial dna integrations.
for example, participants a & b have 16s and 23s rrna integrations into the 50–70 bp of ceacam <dig>  this latter result suggests that the tumors sequenced have a heterogeneous population of cancer cells with different integrations.
characterizing the integrations
these bacterial dna integrations were investigated to determine if they are enriched in absolute position as well as relative position within the utr and first exon of each gene.
however, no pattern for the position of integration into these genes was identified, suggesting that each gene may have a different feature that is affected by, or permissive to, the bacterial dna integrations.
the bacterial and human fragments supporting the integrations are high in guanine, demonstrated by the combination of high %gc and high %ag.
despite the high guanine content, it does not appear that the bacterial dna integrations disrupted g quartets or cpg islands, but the disruption of other g-rich motifs cannot be ruled out.
functional consequences of the bacterial dna integrations
the integrations of the pseudomonas 16s and 23s rrna genes may be important and have a biological effect by altering the transcriptional regulation of the human genes.
one possibility is that the integration of bacterial dna introduced new structures at the breakpoint between bacterial and human dna that alter the transcriptional regulation of these human genes.
given that all of the bacterial integrations have stem-loop structures, secondary structure may play an important role.
alternatively, the integrated bacterial sequence could alter the regulation of these genes by disrupting either the transcript stability, or the availability of the transcriptional and/or translation machinery to bind to these genes.
considering that the integrations are into the 5'-utr of these genes, before the translational start site, it is unlikely that the integrations introduced a frameshift or premature stop codon.
it is also unlikely that these integrations have interfered with splicing since the integrations are not near a splice site recognition sequence.
however, it cannot ruled out that the integrations have created novel transcriptional start sites, translational start sites, or splice sites that would have functional consequences.
CONCLUSIONS
based on the available data, models were generated for the bacterial 16s and 23s rrna fragment integrations into the human ceacam <dig>  ceacam <dig>  cd <dig>  and tmsb <dig>  these mutations are especially intriguing because of their positions in the 5′-utr near the transcriptional start site.
the models presented here lay the groundwork for further in vitro experiments reconstructing these sequences to test if the bacterial rrna integrations may alter the expression of the human genes.

 the advent of next generation sequencing has enabled the interrogation of the human genome and transcriptome with base pair resolution. this new window into the human genome has expanded our understanding of important somatic genome variants such as single nucleotide polymorphisms, chromosomal rearrangements, and exogenous dna integrations. these variants have primarily been identified using paired-end sequencing, a process that relies on breaking the genomic dna into many small fragments and sequencing each fragment of dna from the ends inward toward the center. the resulting pairs of sequencing reads share a unique relationship as they originated from the same randomly sheared dna fragment. by leveraging the unique relationship between the paired-end reads, chromosomal rearrangements and integrations of exogenous dna can be identified and characterized. for example, a rearrangement between two chromosomes would result in at least one chimeric chromosome. by mapping the sequencing reads to a reference genome, a subset of the paired-end reads will support the chromosomal rearrangement by spanning the break point between the two chromosomes with one of the paired-end reads mapping to chromosome a, while the respective other paired-end read maps to chromosome b. in addition to identifying rearrangements, this technique has also been applied to identify integrations into the human genome of: transposable elements, viral genomes, and bacterial dna. if the integration site can be rebuilt in silico, the structure and sequence of the integration site can be used to determine a potential mechanism of integration, elucidate potential functional implications of novel integrations, and lay the foundation for future experiments.fig. a brief schematic illustrates using next generation sequencing to identify structural variation. a recombination event is illustrated between two chromosomes. in order to identify structural variation, paired-end reads are identified that map to different chromosomes indicating a fusion of the two different chromosomes has occurred previously, we identified paired-end reads supporting the integration of pseudomonas-like rrna genes into the 5′-utr of genes in human stomach adenocarcinoma genomes using rna-seq data from the cancer genome atlas. these paired-end reads have one read mapping exclusively to the 16s or 23s rrna genes of a pseudomonas-like bacteria, while the paired read maps uniquely to the 5′-utr of the human ceacam ceacam cd or tmsb genes. as such, these read pairs support bacterial dna integrations by spanning the junction of the bacterial and human dna. to identify these regions in the human genome, >4× human coverage in a single sequencing run was required. despite this level of coverage, single reads that traverse the integration site were not identified, likely owing to the length of the sequencing reads and limits of alignment algorithms. such reads would enable the assembly of the integration with the base pair resolution needed to determine the potential mechanism of integration. instead, here, we have used the paired-end reads to model the most likely structure of the integration of pseudomonas-like rrna gene fragments into ceacam ceacam cd and tmsb  establishing the boundaries of the bacterial rrna integrations in all cases, reads were only recovered spanning one side of the junction between bacterial and human dna, likely because the integration is close to the transcriptional start site. despite numerous attempts, the integrations could not be assembled using the original reads supporting an integration event or with additional reads located near the integrations that were identified using either alignment- or blast-based algorithms. neither split reads, those that are part bacterial and part human on the breakpoint of the integrations, nor soft-clipped reads could be identified in the original alignments or in alignments to a custom reference with both the human genome and the pseudomonas 16s & 23s rrna genes. further blast-based examination of the unmapped read in read pairs that had only one read aligned to the human genome in the region flanking the integration or the pseudomonas rrna gene reference near the integration also failed to identify split reads. to examine this further, a bacterial-human dna integration was constructed with the bacterial dna directly abutting the human dna. a mock dataset was created of all possible combinations of 100-bp paired end reads spanning the integration breakpoint in this artificial sequence mimicking a bacterial-human dna integration. the first read generated was entirely bacterial and ended at the integration breakpoint. each of the subsequent reads in the mock dataset shifted by 1-bp, such that the dataset included a mock read for every position across the integration beginning with an entirely bacterial read and ending with an entirely human read. the second read in the pair was held constant and corresponded to a sequence 225-bp downstream of the break point. lgtseek identified only reads that cover the breakpoint, none of which were soft-clipped as the differences with respect to mapping were similar to those arising from sequencing errors. therefore, we conclude that lgtseek, and more specifically the version of bwa used in lgtseek, is unable to identify reads that span the junction between bacterial and human dna in this data set. given that the bacterial dna integrations could not be assembled, the focus shifted to estimating the location of the bacterial rrna gene fragment integrations into the human genome by examining the structure of the human transcript and the reads supporting the bacterial dna integrations. the integration breakpoint must be downstream of the transcriptional start site of each human gene for three reasons. first, the integrations were identified in an rna-seq data set derived from transcripts so they must be within the transcript boundaries. second, examination of the expression of these genes across all participants from the stad and breast cancer data sets from tcga data available for download from the sra between september 18th–20th, are consistent with an accurately annotated tss. third, the introns and other noncoding regions of the human genome do not have low-level sequencing coverage that would suggest the presence of contaminating genomic dna. therefore, the integration must be downstream of the annotated tss and the tss is the left-most boundary for the possible location of bacterial dna integration, relative to the direction of transcription. the right-most boundary for the bacterial dna integration break point can be delineated from the position of human reads supporting the bacterial dna integrations. more specifically, the bacterial dna integration must be upstream of the human reads supporting the bacterial dna integration, relative to the direction of transcription. therefore, the left most boundary for the site of integration in the human genes, in all cases, should be considered the tss, while the right most boundary position for the integration is the left most position of the consensus sequence for the human reads supporting the bacterial dna integration. using the library insert-size to refine the location of bacterial dna integration the library insert-size and its distribution can be used to refine the location of the bacterial rrna gene fragment integration. in order to model the integrations with the greatest resolution and accuracy, two calculations were used to determine the number of bases needed between the bacterial and human consensus sequences to mirror the library insert-size distribution. suppose that the distance between the bacterial and human fragments is x, where ≤ x ≤, and that each integration has n total reads supporting it. the insert-size, ii, of the ith read pair spanning the junction is calculated based on the known positions of each read within the bacterial and human fragment consensus sequences and the assumed distance, x, between these fragments. the mean absolute value, ad, of the differences between the median sequencing library insert-size, msl, and the read pairs supporting the integration can be calculated as adx=∑i=1niix−msln, where n is the total number of reads supporting the integration. the location of the junction between bacterial rrna and the human gene of interest was estimated to be within the range,, where l is the position of the left end of the human fragment consensus sequence within the human gene and xad is the value of x which has the minimum value of ad.fig. the location of the bacterial dna integration was refined using the library insert-size. in this hypothetical case, the paired-end reads support the integration of a fragment of the bacterial 16s rrna gene into a human gene. using the consensus sequence of the two fragments, the number of bases, xx, needed between the two fragments was titrated so that the insert-sizes of the reads supporting the integration most closely resemble those of the sequencing library. two calculations were used to determine the optimum distance between the two fragments. the first calculation identified the number of bases yielding the minimum average difference between insert-size of the reads supporting the integration, i i, and the median sequencing library insert-size, m sl. the second calculation determined the number of bases resulting in the minimum jensen-shannon distance for the insert-size distributions of the sequencing library, d sl, and the reads supporting the integration, d i the second approach of estimating the distance between the bacterial and human consensus sequences uses a dissimilarity measure called the jensen-shannon distance, jsd, between the insert-size distributions of the sequencing library, dsl, and reads supporting the integration, di. the jensen-shannon distance between two distributions takes values between and with implying the identity of the distributions and indicating major discordance between the distributions. as in the average distance case, the location of the junction between bacterial 16s rrna fragment and the human gene of interest was estimated to lie within the range,, where l is the position of the left end of the human fragment consensus sequence within the human gene and xjs is the distance value for the minimum value of jsd. the majority of the time the ad and jsd calculations were in agreement. when the two calculations differed, the jsd was preferred since it accounts for the distribution of insert size values for the two populations. validating the method to validate this method of modeling the bacterial dna integrations, a bacterial dna integration in a cancer cell line genome was identified, modeled, and experimentally validated. using the cancer cell line encyclopedia and a similar cancer cell line sequencing project that is publicly available through genentech inc, paired-end rna-seq reads were identified that support the likely in vitro integration of a kanamycin resistance gene near the kbp position in chromosome of the kpl- cell line. the paired-end reads supporting this integration spanned the junction of bacterial and human dna at the 3' - side of the kanamycin gene, with respect to the direction of transcription. given that this integration was likely in vitro, a viral promoter is expected to drive the kanamycin gene expression, which would prevent the bacterial sequence focused pipeline from detecting the 5' - junction between the integrated bacterial dna and human chromosome. in both independently sequenced datasets, the jsd and the ad predicted that there should been 0 bp between the consensus bacterial and human dna. subsequently, the genomic dna of the kpl- cell line was acquired from leibniz institute dsmz and the junction between the integrated kanamycin resistance gene and chromosome was pcr amplified, cloned, and sequenced. the sequencing revealed that the jsd method correctly predicted the sequence and approximate position of the integration site, but that the exact position was off by 2 bp, with the actual position of integration being 2 bp from that predicted by the model.fig. a bacterial dna integration into the kpl- cell line genome was modeled using the jsd and ad calculations. based on the data from the cancer cell line encyclopedia and a similar data set from genentech, the jsd and ad methods predicted that the bacterial & human consensus sequences are 0 bp apart. however, pcr amplification and sequencing revealed an additional 2 bp of sequence between the two consensus sequences as illustrated. consistent with this, a plot of the sequence coverage of the reads across the junction is relatively constant. no further insertions or deletions are observed in the underlying reads kpl breakpoint while investigating the kpl integration, soft-clipped reads at the kpl bacteria-human dna junction were identified manually. to better understand why these reads were not identified by the lgtseek pipeline, a mock dataset was created of all possible combinations of 100-bp paired end reads spanning the cloned and sequenced integration breakpoint. the first read generated was entirely bacterial and ended at the integration breakpoint. each of the subsequent reads in the mock dataset shifted by 1-bp, such that the dataset included a mock read for every position across the integration beginning with an entirely bacterial read and ending with an entirely human read. the second read in the pair was held constant and corresponded to a sequence 190-bp downstream of the break point, which was selected based on the insert-size distribution of the ccle data. no split reads were identified from running the entire kpl dataset spiked with this mock dataset through the lgtseek pipeline due to a 24-bp inversion in the bacterial dna near the kpl breakpoint. bacterial rrna gene integrations into ceacam & ceacam6 from the tcga stad data set that was previously analyzed, two participants have paired-end reads supporting pseudomonas dna integrating into exon of ceacam. in participant a, the integration of the pseudomonas 16s rrna gene fragment is supported by paired-end reads. using the jsd, the calculated distance between the bacterial and human fragment consensus sequences is 26 bp, with a minimum jsd value of, while the ad placed it between 11–17 bp. a visual inspection of the read insert-size distribution over a variety of values of x further supports the jsd prediction. therefore, given that the human fragment consensus sequence starts at position in the ceacam gene, the bacterial dna integration is estimated to occur between positions 48– in the ceacam gene.fig. the model of the pseudomonas 16s rrna gene fragment integration into ceacam of participant a is presented. the structure of the pseudomonas 16s rrna gene fragment, the first exon ceacam and the paired reads supporting the integration of the rrna gene fragment into ceacam are illustrated. calculations using the jsd and ad support that there are 26 bp between the pseudomonas 16s rrna & ceacam fragments and that the integration is in 48–73 bp of ceacam. the ceacam sequence upstream of the bacterial 16s rrna fragment is for illustrative purposes only. the insert-size for each paired-end read is color-coded with lighter colors being closer to the median insert-size, red designating those insert-sizes larger than the median, and green for insert-sizes less than the median as further clarified in additional file 7: figure s the optimum model for the structure of the integration is compared to alternative distances between the two fragments in, where an additional,, −, & - median absolute deviations are placed between the two fragments represented by the gray region. the actual dna sequence of the gray region is unknown participant b has paired-end reads supporting the pseudomonas 16s rrna gene integration into ceacam. both the jsd and the ad calculations are in agreement that the distance between the two fragments is 0 bp, and that the 16s rrna gene fragment has integrated at the 18 bp in ceacam however, manually inspecting the insert-size of the reads suggests that the human read farthest upstream is an outlier that does not fit the insert-size distribution of the other reads. after removing this potential outlier, a model of the integration was calculated. the jsd and ad support that the distance between the consensus sequences is 7 bp. in addition, both the jsd and ad decreased, suggesting that this may be a more accurate representation of the integration, particularly given that decreasing the number of reads in the analysis typically increases the jsd as described below. therefore, the integration of the 16s rrna gene fragment is likely between positions 47–53 bp of ceacam for the subsequent comparisons in this manuscript, only this latter model is discussed.fig. a model of the pseudomonas 16s rrna gene fragment integration at the 18 bp of ceacam from participant b is illustrated. based on the jsd and ad calculations, the model has 0 bp between the bacterial and human fragment consensus sequences. after removing the potential outlier paired-end read, the optimum structure for the integration is illustrated in d. by removing the potential outlier read, both the jsd and ad decreased, supporting that there are 7 bp between the two fragments. this indicates the pseudomonas 16s rrna likely integrated between positions 45–53 bp of ceacam5 names indicate the participant with letters a through f, their respective rrna gene integration, and the human gene it has integrated into if the potential outlier paired-end read is excluded, the modeling suggests that the 16s rrna integrations into the ceacam gene of participants a & b are nearly identical, with the models only differing by a few bp. in addition, the human reads from participants a & b map to similar positions in ceacam. in addition to the 16s rrna gene integration described above, participant a also has paired-end reads supporting an integration of 1626–1699 bp from the pseudomonas 23s rrna gene into the ceacam gene. the jsd and ad calculations support that the integration is between positions 49– of ceacam in addition, participant b has reads supporting integration of a nearly identical fragment of the 23s rrna gene integrating into an almost identical region in ceacam. therefore, the data support two different pseudomonas dna integrations into approximately the same position in ceacam for both participants a & b. similar to ceacam participants a & c have similar pseudomonas 16s rrna gene fragment integrations into the related ceacam. both participants have nearly identical fragments of the 16s rrna gene integrated into ceacam at the same predicted position. in addition, these ceacam integrations are similar to the position of the 16s & 23s rrna integrations and the sequence of the 16s rrna integrations into ceacam cd integrations both participants d and e have fragments of the pseudomonas 16s and 23s rrna integrated in the cd gene. both participants have almost identical fragments of the 16s rrna gene integrated into overlapping positions in the cd gene. for participant d, the calculated model supports that the 16s rrna fragment integrated between positions 180–187 bp of cd while the fragment of 16s rrna integrated into participant e’s cd gene between positions 163–205 bp. the data supporting these two integrations are consistent and suggests that the two participants have similar integrations. the data also support a second fragment further downstream in the 16s rrna gene integrating into positions 136–147 bp of cd in participant d. there are different fragments of the 23s rrna gene fragments with support for integration into the cd gene of participants d & e. all three integrations have similar locations in the cd gene. of the three integrations, two integrations, one each from participant d & e, are predicted to have integrated into the same region of cd. the remaining 23s rrna gene fragment with support for integration into cd is predicted to have integrated into the 201–225 bp of cd of participant e. tmsb integrations unlike the bacterial dna integrations into ceacam ceacam and cd the tmsb gene has only one participant with the integration of a pseudomonas 16s rrna fragment. the model estimates that the 16s rrna fragment integrated into the first 13 bp of tmsb. calibrating the jsd as described previously, a calculated jsd value of supports that the two populations are identical, while a value of indicates major discordance between the distributions. the average jsd for the models of bacterial dna integration presented here was. therefore, the jsd calculation was investigated further to determine why it was not closer to zero. to test the effect that the number of reads had on the jsd, the jsd was calculated for samples for specific values of k, such that ≤ k ≤, where k represents the number of randomly sampled read pairs that mapped to the ceacam transcript from the sequencing library of participant b. samples with 3– randomly selected reads had a jsd of ~, consistent with the results for modeling the integrations. as the number of randomly sampled reads increases, the jsd decreases. the average jsd continues to decrease as the number of reads increases until around reads when the jsd starts to plateau, eventually reaching a minimum of at reads, the maximum number tested.fig. boxplots illustrate the jsd calculated for subsets of rna-seq data mapped to the human transcript. the insert-size distribution of the sequencing library of participant b was compared to the distribution of insert-sizes of specific values of randomly selected reads that mapped to the ceacam transcript from the sequencing library of participant b. panel a has a discontinuous x-axis to better illustrate the jsd over the various number of reads tested; panel b illustrates the same data with a continuous x-axis guanine rich sequences while modeling the integrations, it was observed that both the bacterial and human fragment consensus sequences seemed to be particularly rich in guanine. consistent with this observation, the median %gc for the bacterial and human fragment consensus sequences was % and  %, respectively. the %gc and %ag were compared for the consensus sequences supporting the bacterial dna integrations using the wilcoxon rank-sum test. the bacterial and human consensus sequences were significantly different than each other. both sets of consensus sequences were also significantly different than the complete refseq pseudomonas genomes and the human genomic reference. the human fragment consensus sequences were also significantly different that the pseudomonas 16s & 23s rrna genes available in the silva database. considering that the predicted integration sites are typically into the first 150 bp of each human transcript, the %gc of the fragment consensus sequences supporting the integrations was compared with the first 150 bp of each transcript from the human transcriptome reference and found that the %gc for the human fragment consensus sequences was significantly different. given that the expressed transcripts may have a different gc profile, the %gc for the fragment consensus sequences was compared to the expressed transcripts from the participants with the bacterial dna integrations. the %gc of the human consensus sequences was significantly different than the expressed transcripts. next, by stratifying the expressed genes into quartiles based on their rpkm, the %gc for the quartiles of expressed transcripts was compared to consensus sequences. the human fragment consensus sequences were significantly different than the two most highly expressed quartiles. in addition, with a less stringent but acceptable significance level of, the %gc for the lowest two quartiles are likely significantly different than the human fragment consensus sequences.table wilcoxon rank-sum test comparison of gc & ag content * = p < ; na not applicable in addition to gc-content, the enrichment of guanine suggests that these are purine-rich sequences, which have been shown to have a role in regulating transcription. therefore, the percentage of bases that are purines for each fragment was calculated using the transcribed strand of each gene and compared to various references using the wrt. the median %ag for the human and bacterial fragment consensus sequences was &, respectively. the %ag of the bacterial fragment consensus sequences was significantly different than the pseudomonas 16s & 23s rrna gene database, and both of the complete and first 150 bp of the expressed human transcripts in the participants with the bacterial dna integrations. on the other hand, the %ag of the human fragment consensus sequences were statistically different than the expressed transcripts from the participants with the bacterial dna integrations. due to the observed higher guanine content of these sequences, the human genes with bacterial dna integrations were investigated for guanine associated motifs. a database search of published g-quartets in the human genome did not reveal evidence for g-quartets near the integrations. examination of the first exon and kbp upstream of the tss for ceacam ceacam and cd did not reveal cpg islands. on the other hand, all three algorithms support cpg islands in the kbp region upstream of the tmsb tss. however, only of the algorithms suggested the cpg island may extend downstream of the tss near the predicted integration into tmsb. cpg islands are typically defined as ≥200 bp regions that have more observed than expected cpg dinucleotides and high gc content. next, the window size for identifying cpg islands was decreased to refine the location of the cpg island signal in the first exon of tmsb the cpg signal was determined to originate from the 3' end of the first exon of tmsb while the side of the first exon near the predicted integration site had no cpg signal. therefore, it is unlikely the bacterial integration into tmsb disrupted cpg sites. rrna structure based on the position of the bacterial dna integrations into 5′-utr of the human genes, it is possible that the bacterial dna integrations altered the expression of the human genes by disrupting a regulatory element in the 5′-utr. one reason to expect that this might be the case is that the integrations arise from the 16s & 23s rrna genes that contain numerous stem-loop structures that have the potential to alter secondary structure and thus transcription. therefore, the consensus sequence of the integrated bacterial dna was examined for the presence of stem-loop structures from the bacterial rrna gene by mapping the bacterial fragments onto the known secondary structure. it was concluded that 11/ integrations have large stem-loop structures in the middle of the bacterial fragments; while the remaining three fragments of integrated bacterial dna have stem-loop structures near the ends of each fragment. subsequently, the predicted rna structure of the human transcripts was examined to determine if the bacterial integrations may have disrupted the secondary structure of the human transcripts. by mapping the predicted position of each integration onto the ceacam ceacam cd and tmsb predicted mrna transcript structures, it was determined that the integrations are either in, or near, complex secondary structures. relative position of integration to explore if the bacterial dna integrations may be occurring in the same relative position within each human gene, the relative position of the predicted site for integration within the first exon was calculated. while the ceacam and ceacam integrations are enriched around the middle of the first exon, the two genes also share homology and have similar structures. the ceacam genes, cd and tmsb all have different relative position of integrations.table relative position of the bacterial dna integrations discussion obstacles in assembling the integrations here, the jensen-shannon distance between the insert-size distributions of the paired-end reads supporting the pseudomonas-like dna integrations relative to the respective sequencing library was used to model bacterial dna integration in ceacam ceacam cd & tmsb the jsd values supporting these models are within the expected range of the jsd values given the number of reads supporting the bacterial dna integrations. ideally, the sequencing reads would be assembled to determine the exact integration site and sequence. there are likely three compounding factors that prevent the assembly of these bacterial dna integrations. first, the reads are only 51 bp. in order to find a read that is on the integration breakpoint, the read on the breakpoint needs to be able to be split into uniquely human and bacterial portions. as a result of the short read length, a read would have to be exactly half human and half bacterial in order to give a sufficient sequence length of ~25 bp to identify each portion. this is highly unlikely to occur. second, there is limited coverage of the integrations. the limited coverage is likely due to the proximity of the integrations to the beginnings of transcripts where rna-seq coverage is poor. in addition, it is possible that the low rna-seq coverage is the result of the bacterial dna integrations inducing decreased expression. lastly, the data suggest a heterogeneous population of transcript variants in a tumor that makes de novo assembly difficult. in a mock-experiment with an artificial bacterial-human dna integration, lgtseek identified only reads that cover the breakpoint. these reads align to the bacterial dna, with only 1–3 bp at the 3'-end of the read covering the breakpoint. lgtseek uses a version of bwa that may not identify these reads. while a more recent version of bwa identifies soft-clipped reads better, it also works on the assumption that all reads should match the reference and as such has been found to erroneously align bacterial reads to the human reference. in kpl the cloned and sequenced bacteria-human dna junction was used to generate a similar mock dataset that demonstrates that split and soft-clipped reads could not be identified due to a 24-bp inversion. given that when the bacterial and human dna are flush, only 3 % of read pairs are identified, and that there are only read pairs supporting stad bacterial dna integrations on average, we do not expect that soft-clipped reads would be identified even if the bacterial dna is flush to the human dna. furthermore, such soft clipped reads that do exist would be difficult to distinguish from sequencing error or sequencing polymorphism between the bacterial sequence and the reference genome. additionally, the average distance between the bacterial and human fragments of dna in the models described here was found to be ~24 bp, which suggests that the two may not be flush, further explaining the difficulty in detecting such reads. heterogeneous integrations participants a through e all have data supporting multiple integrations of different fragments of the 16s and 23s rrna pseudomonas gene into the same human gene within each participant. in addition, the data support multiple participants with similar, if not identical, bacterial dna integrations. for example, participants a & b have 16s and 23s rrna integrations into the 50–70 bp of ceacam this latter result suggests that the tumors sequenced have a heterogeneous population of cancer cells with different integrations. characterizing the integrations these bacterial dna integrations were investigated to determine if they are enriched in absolute position as well as relative position within the utr and first exon of each gene. however, no pattern for the position of integration into these genes was identified, suggesting that each gene may have a different feature that is affected by, or permissive to, the bacterial dna integrations. the bacterial and human fragments supporting the integrations are high in guanine, demonstrated by the combination of high %gc and high %ag. despite the high guanine content, it does not appear that the bacterial dna integrations disrupted g quartets or cpg islands, but the disruption of other g-rich motifs cannot be ruled out. functional consequences of the bacterial dna integrations the integrations of the pseudomonas 16s and 23s rrna genes may be important and have a biological effect by altering the transcriptional regulation of the human genes. one possibility is that the integration of bacterial dna introduced new structures at the breakpoint between bacterial and human dna that alter the transcriptional regulation of these human genes. given that all of the bacterial integrations have stem-loop structures, secondary structure may play an important role. alternatively, the integrated bacterial sequence could alter the regulation of these genes by disrupting either the transcript stability, or the availability of the transcriptional and/or translation machinery to bind to these genes. considering that the integrations are into the 5'-utr of these genes, before the translational start site, it is unlikely that the integrations introduced a frameshift or premature stop codon. it is also unlikely that these integrations have interfered with splicing since the integrations are not near a splice site recognition sequence. however, it cannot ruled out that the integrations have created novel transcriptional start sites, translational start sites, or splice sites that would have functional consequences.  based on the available data, models were generated for the bacterial 16s and 23s rrna fragment integrations into the human ceacam ceacam cd and tmsb these mutations are especially intriguing because of their positions in the 5′-utr near the transcriptional start site. the models presented here lay the groundwork for further in vitro experiments reconstructing these sequences to test if the bacterial rrna integrations may alter the expression of the human genes. 
modeling of these integrations suggests that the previously identified sequences cover most of the sequence flanking the junction between the bacterial and human dna. the functional and biological consequences of these bacterial dna integrations remain unknown. further examination of these reads reveals that these integrations are rich in guanine nucleotides and the integrated bacterial dna may have complex transcript secondary structures. the models presented here lay the foundation for future experiments to test if bacterial dna integrations alter the transcription of the human genes. cancer is a disease driven by the accumulation of genomic alterations, including the integration of exogenous dna into the human somatic genome.
 cancer is a disease driven by the accumulation of genomic alterations, including the integration of exogenous dna into the human somatic genome. we previously identified in silico evidence of dna fragments from a pseudomonas-like bacteria integrating into the 5′-utr of four proto-oncogenes in stomach cancer sequencing data. the functional and biological consequences of these bacterial dna integrations remain unknown.  modeling of these integrations suggests that the previously identified sequences cover most of the sequence flanking the junction between the bacterial and human dna. further examination of these reads reveals that these integrations are rich in guanine nucleotides and the integrated bacterial dna may have complex transcript secondary structures.  the models presented here lay the foundation for future experiments to test if bacterial dna integrations alter the transcription of the human genes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords dna integrationsomatic genomegenome variationcancer genomicshost-bacteria interactionshttp://dx.doi.org/ /100000052nih office of the director 1-dp2-od007372dunning hotopp julie c. http://dx.doi.org/http://dx.doi.org/ /100000062national institute of diabetes and digestive and kidney diseases gi-t dk067872-10http://dx.doi.org/ /100000054national cancer institute1-r01-ca206188dunning hotopp julie c. issue-copyright-statement© the author 2016 