<jupyter_start><jupyter_text>IntroductionThe BioGraph format is a compressed index of sequencing data that holds a rapidly querable representation of reads. The BioGraph has two main components. The Seqset and the Readmap. Briefly, the Seqset is the compressed index of every kmer from the foward and reverse complement of all reads in the sequencing experiment. Each of these kmers is called a SeqsetEntry or simply an entry. The Readmap is a lookup of what entries correspond to original input reads and keep information such as the read's length or if the read has a mate-pair. This tutorial is an introduction to basic queries and functionality of the BioGraphSDK. SetupTo follow along with this tutorial, the BioGraphSDK should have been installed and the `bgtools install_tests` have been run. For information on how to do this, see **Spiral Genetics’ BioGraph Toolkit and SDK User Guide** in the _Installation_ section About the DataA handful of SNPs, indels, and SVs were randomly simulated in the "J02459.1 Enterobacteriaphage lambda, complete genome". These variants were then given genotypes mocking a familytrio. The sample genomes for the proband, mother, and father can be found in`references/samples/` as well a vcf in `variants/family.vcf`.For each sample, 150bp paired-end reads were simulated using [mason](http://www.seqan.de/apps/mason/).For the examples, we use the lambdaToyData. Point to the benchmark data or the results from a run of `bgtools install_tests` to follow along.  Opening a BioGraphTo start, let's look at the single sample proband_lambda.bg by opening it as a `BioGraph` object and looking at some of its metadata.Update the `LAMBDA_PATH` variable to point to the root directory where converted lambdaToyData exists.<jupyter_code>LAMBDA_PATH = "/scratch/lambda/lambdaToyData/benchmark"
import os
import biograph as bgsdk
#Load the BioGraph
my_bg = bgsdk.BioGraph(os.path.join(LAMBDA_PATH, "proband_lambda.bg"))
ref = bgsdk.Reference(os.path.join(LAMBDA_PATH, "ref_lambda"))

print("BioGraph Version", my_bg.metadata.version)
print("AccessionId:", my_bg.metadata.accession_id)
print("BioGraphId:", my_bg.metadata.biograph_id)
print("Sample Info", my_bg.metadata.samples)<jupyter_output>BioGraph Version 4.0.1-dev
AccessionId: proband
BioGraphId: 52f41046-0e66-9b4a-b5a9-00a1a8029aa5
Sample Info {'proband': '47b43dbb0306fb5da00bca9c33b2fb3de9db7bf4'}<jupyter_text>AccessionId is a user-specified sample identifier for the data.BioGraphId is an internal sha1 that's used to identify the sample for merging.Sample Info is a dictionary holding Key/Value of AccesionId to BioGraphId.We also looked at a couple propeties of the reads contained within the BioGraph.The AccessionIds are most useful for BioGraphs comprising multiple, merged samples. We'll cover details about merged BioGraphs in a later section.  Seqsets and ReadmapsBioGraphs contain two main data structures: Seqsets and Readmaps.  A readmap lists information about all reads that were seen.  A Seqset contains a list of all sequences and subsequences of bases present in any readmap.The Seqset data structure is designed for efficient searching of arbitrary sequences, and effecient traversal of overlapping sequences.The Readmap data structure tracks which Seqset entries have reads that support them, and associated data about those reads such as pairing data.  Using a SeqsetYou can use the "find" method to look up a seqset entry by sequence, and the "sequence" method to look up the sequence associated with a seqset entry.  However, try to avoid using these methods in performance critical code paths; it is much more efficient to use the seqset traversal methods such as "push_front" and "pop_front".  In fact, "find" uses "push_front" internally, and "sequence" uses "pop_front" internally.<jupyter_code># Look up a seqset entry by sequence.
seqset = my_bg.seqset
print("This seqset has %d entries" % (seqset.size()))
query_seq = bgsdk.Sequence("GTAATCTTTTAAT" + "T" + "TTAAATAAGTTA")
missing_query_seq = bgsdk.Sequence("GTAATCTTTTAAT" + "C" + "TTAAATAAGTTA")
for seq in [query_seq, missing_query_seq]:
    entry = seqset.find(seq)
    if entry:
        print("Sequence %s is present in the seqset" % (entry.sequence()))
    else:
        print("Sequence %s is not present in the seqset" % (seq))<jupyter_output>This seqset has 97340 entries
Sequence GTAATCTTTTAATTTTAAATAAGTTA is present in the seqset
Sequence GTAATCTTTTAATCTTAAATAAGTTA is not present in the seqset<jupyter_text>However, it is much more efficient to use the seqset traversal primitives to find neighboring entries if you're not looking up a whole sequence of bases from scratch.<jupyter_code>seq = bgsdk.Sequence("AATCTTTTAATTTTAAATAAGTTA")
entry = seqset.find(seq)
for base in "TGATC":
    seq = base + seq
    entry = entry.push_front(base)
    if entry:
        print("%s exists in the seqset" % (seq))
    else:
        print("%s does not exist in the seqset" % (seq))
        break<jupyter_output>TAATCTTTTAATTTTAAATAAGTTA exists in the seqset
GTAATCTTTTAATTTTAAATAAGTTA exists in the seqset
AGTAATCTTTTAATTTTAAATAAGTTA exists in the seqset
TAGTAATCTTTTAATTTTAAATAAGTTA does not exist in the seqset<jupyter_text>Using a Readmap If your BioGraph only contains a single sample, you will not have to specify the accession ID when opening the readmap.  However, BioGraphs can contain multiple samples; see the section on Merged BioGraphs for details.<jupyter_code>my_rm = my_bg.open_readmap()

# Print general statistics about the readmap
print("Readmap contains %d reads and %d bases" % (my_rm.get_read_count(), my_rm.get_num_bases()))
p = my_rm.get_pair_stats()
print("\tPaired\tUnpaired")
print("Reads\t{p.paired_reads}\t{p.unpaired_reads}".format(p=p))
print("Bases\t{p.paired_bases}\t{p.unpaired_bases}".format(p=p))
print("Including both forward and reverse direction reads, readmap contains %d entries" % (my_rm.size()))<jupyter_output>Readmap contains 48956 reads and 7235142 bases
	Paired	Unpaired
Reads	39974	8982
Bases	5908084	1327058
Including both forward and reverse direction reads, readmap contains 97912 entries<jupyter_text>Looking up readsIf you've found a seqset entry and you'd like to see what reads relate to it, there are two methods available.One method lets you find all reads that are entirely contained within a seqset entry, and the other lets you find all reads that start with a given seqset entry and are entirely contained within it.Once you have a specific read, you can use "get_seqset_entry" to translate it back to its seqset entry.<jupyter_code># Find all reads that start with the given seqset entry and are entirely contained within it
query_seq = bgsdk.Sequence("AAAGAAGATTTCCAATAATCAGAACAAGTCGGCTCCTGTTTAGTTACGAGCGACATTGCTCCGTGTATTCACTCGTTGGAATGAATACACAGTGCAGTGTTTATTCTGTTATTTATGCCAAAAATAAAGGCCACTATCAGGCAGCTTTGT")
print("Original sequence:")
print(query_seq)
print("")
entry = seqset.find(query_seq)
for read in my_rm.get_prefix_reads(entry):
    if read.is_original_orientation():
        print("Found a read of length %d.  This is how it appeared in the input data:" % (len(read)))
        print(read.get_seqset_entry().sequence())
    else:
        print("Found a read of length %d.  It appeared in the input data as its reverse complement:" % (len(read)))
        print(read.get_seqset_entry().sequence().rev_comp())
    if read.has_mate():
        print("  The preceding read has a mate: %s" % (read.get_mate().get_seqset_entry().sequence()))
    else:
        print("  The preceding read was not part of a pair, however its reverse complement is: %s" % (read.get_rev_comp().get_seqset_entry().sequence()))
    print("")
# Find all reads that contain a given seqset entry anywhere in the read.
query_seq = bgsdk.Sequence("CTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAAC")
display_offset = 11
print("%s%s" % (" "*display_offset, query_seq))
print("")
entry = seqset.find(query_seq)
for (offset, read) in my_rm.get_reads_containing(entry):
    print("%s%s" % (" "*(display_offset - offset), read.get_seqset_entry().sequence()))<jupyter_output>CTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAAC

           CTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAACCAGGTTCGTG
         GACTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAACCAGGTTCG
        TGACTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAACCAGGTTC
       ATGACTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAACCAGGTT
      CATGACTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAGGCCTGTTATTAGCTCAGTAATGTAGATGGTCATCTTTTAACTCCATATACCGCCAATACCCGTTTCATCGCGGCACTCTGGCGACACTCCTTAAAAACCAG
      CATGACTTCCCTCTCCCCCAAATAAAAAGGCCTGCGATTACCAGCAG[...]<jupyter_text>Merged BioGraphsMerged BioGraphs can have multiple samples present in the BioGraph.  In this case, the seqset will contain the union of all sequences contained in any input sample, and there will be a separate readmap for each sample.For these examples, we'll use the `family_lambda.bg` and take a look at all the properties.<jupyter_code>my_bg = bgsdk.BioGraph(os.path.join(LAMBDA_PATH, "family_lambda.bg"))
print("AccessionId:", my_bg.metadata.accession_id)
print("BioGraphId:", my_bg.metadata.biograph_id)
print("Number of Samples", len(my_bg.metadata.samples))
for sample in my_bg.metadata.samples:
    print("-- Opening Sample", sample, "--")
    readmap = my_bg.open_readmap(sample)
    print("   Number of Reads:", readmap.get_read_count())
    print("   Number of Bases:", readmap.get_num_bases())<jupyter_output>AccessionId: family
BioGraphId: 3f64b88d-e1b5-f29d-1569-256e95595e6d
Number of Samples 3
-- Opening Sample mother --
   Number of Reads: 48964
   Number of Bases: 7235912
-- Opening Sample father --
   Number of Reads: 48930
   Number of Bases: 7228335
-- Opening Sample proband --
   Number of Reads: 48956
   Number of Bases: 7235142<jupyter_text>Once you open a specific readmap, you can use all the same seqset and readmap calls in a merged BioGraph as you can use in a BioGraph with only one sample.  Using a ReferenceSince BioGraph files are reference agnostic, there is no information about a reference inside of it. Instead, information between the BioGraph and a Reference are pieced together on the fly. For example, let's use the lambda reference to find all of the reads that contain a portion of its chromosome.<jupyter_code>my_ref = bgsdk.Reference(os.path.join(LAMBDA_PATH, "ref_lambda"))
ref_range = my_ref.make_range("lambda", 19575, 19600)
print("Looking for reads containing %s:%d-%d (%s)" % (ref_range.chromosome, 
                                                      ref_range.start, 
                                                      ref_range.end, 
                                                      ref_range.sequence))
entry = my_bg.seqset.find(ref_range.sequence)
for sample in my_bg.metadata.samples:
    readmap = my_bg.open_readmap(sample)
    reads = readmap.get_reads_containing(entry)
    print("%d matching reads found in %s" % (len(list(reads)), sample))<jupyter_output>Looking for reads containing lambda:19575-19600 (TAAATTCTGATTAGCCAGGTAACAC)
119 matching reads found in mother
107 matching reads found in father
114 matching reads found in proband<jupyter_text>Additionally, we can query the Reference for locations where a sequence exists. Note that this is more of a mapping than a proper alignment. This means that we report all of the places in the reference a sequence exactly matches to the direct strand. There currently is no method to do a 'fuzzy' match like what's needed for traditional alignment operations.<jupyter_code>lookup = my_ref.find("CGTGCTGTC")
print(lookup.matches, "matches to reference found.")
print("chrom start end")
for i in range(lookup.matches):
    mat = lookup.get_match(i)
    print(mat.chromosome, mat.start, mat.end)<jupyter_output>2 matches to reference found.
chrom start end
lambda 12589 12598
lambda 537 546
{'fraction_non_alphanumeric': 0.06905184945889194, 'fraction_numerical': 0.01663705378775642, 'mean_word_length': 5.1915, 'pattern_counts': {'":': 1, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '489931', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4119, 'n_tokens_neox': 3811, 'n_words': 1306, 'path': 'python/jupyter/biograph/BioGraph_SDK_Tutorial.ipynb', 'repo_name': 'spiralgenetics/biograph', 'size': '12382'}