Genomic ranges¶

Open this tutorial as a notebook

A SeqPosition is a single site, optionally tied to a SeqRecord holding the sequence it came from. An SRange is an interval between two of them, and an SRangeSet is a sorted collection of intervals. Together they cover much of what you would otherwise reach for bedtools to do.

All coordinates here are 1-based and inclusive at both ends, matching GFF and BED conventions rather than Kotlin's 0-based array indexing.

%use biokotlin

import biokotlin.genome.*
import biokotlin.seq.*
import org.jetbrains.kotlinx.dataframe.api.print

Sequences and sequence records¶

Seq and NucSeq build a nucleotide sequence; ProteinSeq builds a peptide. A NucSeqRecord pairs a sequence with an identifier, which is what ranges refer back to.

val seq = Seq("GCAGAT")

val rec1 = NucSeqRecord(NucSeq("ATAACACAGAGATATATC"), "1")
val rec1a = NucSeqRecord(seq, "1")
println(rec1)
println(rec1a)

ID: 1
Sequence: ATAACACAGAGATATATC

ID: 1
Sequence: GCAGAT

Subsetting a sequence¶

Index a record with an IntRange to pull out part of its sequence. The range is inclusive at both ends, so 1..6 returns six bases.

rec1[1..6]

TAACAC

Positions¶

A SeqPosition may be created with or without a record. With one, the site must fall inside the stored sequence.

val seqPos1 = SeqPosition(rec1, 8)
println(seqPos1)

val seqPosNoRecord = SeqPosition(null, 8)
println(seqPosNoRecord)

1:8
null:8

Ranges¶

Build an SRange from a record and a coordinate range, or from two SeqPosition objects directly.

val sRange = rec1.range(8..12)
println(sRange)

// The `..` operator also works directly on positions.
val sRange2 = seqPos1..seqPos1.plus(4)
println(sRange2)

1:8..1:12
1:8..1:12

Flanking¶

Flanking works like bedtools flank. flankBoth returns two new intervals, one on each side of the original; flankLeft and flankRight return one. When the range carries a record, flanking is clipped to the length of that sequence, so no interval runs off the end of the chromosome.

val sRangeFlanked = sRange.flankBoth(5)
println("sRange:        $sRange")
println("flanked both:  $sRangeFlanked")

sRange:        1:8..1:12
flanked both:  [1:3..1:7, 1:13..1:17]

// The sequence in rec1 is only 18 bases long, so flanking right by 10 from
// position 12 stops at 18 rather than reaching 22.
val sRangeFlankRight = sRange.flankRight(10)
println("flanked right: $sRangeFlankRight")

val sRangeFlankLeft = sRange.flankLeft(4)
println("flanked left:  $sRangeFlankLeft")

flanked right: 1:13..1:18
flanked left:  1:4..1:7

Sets of ranges¶

An SRangeSet keeps its ranges sorted by a comparator you supply. SeqRangeSort provides the usual building blocks: numberThenAlphaSort orders sequence identifiers the way a person would, and leftEdge orders ranges on the same sequence by their start.

val dnaString = "ACGTGGTGAATATATATGCGCGCGTGCGTGGATCAGTCAGTCATGCATGCATGTGTGTACACACATGTGATCGTAGCTAGCTAGCTGACTGACTAGCTGACCGTACGTACGTATCAGTCAGCTGACACGTGGTGAATATATATGCGCGCGTGCGTGGATCAGTCAGTCATGCATGCATGTGTGTACACA"
val dnaString2 = "ACGTGGTGAATATATATGCGCGCGTGCGTGGACGTACGTACGTACGTATCAGTCAGCTGAC"
val dnaString3 = "TCAGTGATGATGATGCACACACACACACGTAGCTAGCTGCTAGCTAGTGATACGTAGCAAAAAATTTTTT"

val record1 = NucSeqRecord(NucSeq(dnaString), "Seq1")
val record2 = NucSeqRecord(NucSeq(dnaString2), "Seq2")
val record3 = NucSeqRecord(NucSeq(dnaString3), "Seq3")

val sort = SeqRangeSort.by(SeqRangeSort.numberThenAlphaSort, SeqRangeSort.leftEdge)

val set1 = nonCoalescingSetOf(
    sort,
    record1.range(27..40), record1.range(44..58), record1.range(1..15),
    record3.range(18..33), record2.range(3..13), record2.range(25..35)
)
println("SRangeSet 1:")
set1.toDataFrame().print()

val set2 = nonCoalescingSetOf(
    sort,
    record1.range(30..35), record1.range(40..50), record1.range(18..22),
    record3.range(1..10), record2.range(10..13), record2.range(45..55)
)
println("SRangeSet 2:")
set2.toDataFrame().print()

SRangeSet 1:
     ID start end  range
 0 Seq1     1  15  1..15
 1 Seq1    27  40 27..40
 2 Seq1    44  58 44..58
 3 Seq2     3  13  3..13
 4 Seq2    25  35 25..35
 5 Seq3    18  33 18..33

SRangeSet 2:
     ID start end  range
 0 Seq1    18  22 18..22
 1 Seq1    30  35 30..35
 2 Seq1    40  50 40..50
 3 Seq2    10  13 10..13
 4 Seq2    45  55 45..55
 5 Seq3     1  10  1..10

Intersections¶

intersect returns the overlapping portions of two sets, the same way bedtools intersect does.

val intersections = set1.intersect(set2)
println("intersection size: ${intersections.size}")
intersections.toDataFrame().print()

intersection size: 4
     ID start end  range
 0 Seq1    30  35 30..35
 1 Seq1    40  40 40..40
 2 Seq1    44  50 44..50
 3 Seq2    10  13 10..13

Coalescing and non-coalescing sets¶

A coalescing set merges overlapping ranges as they are added; a non-coalescing set keeps them separate. Below, the two ranges on Sequence 1 overlap, so only the coalescing set merges them.

val recordA = NucSeqRecord(
    NucSeq("ACGTGGTGAATATATATGCGCGCGTGCGTGGATCAGTCAGTCATGCATGCATGTGTGTACACACATGTGATCGTAGCTAGCTAGCTGACTGACTAGCTGAC"),
    "Sequence 1", description = "The first sequence", annotations = mapOf("key1" to "value1")
)
val recordB = NucSeqRecord(
    NucSeq("ACGTGGTGAATATATATGCGCGCGTGCGTGGACGTACGTACGTACGTATCAGTCAGCTGAC"),
    "Sequence 2", description = "The second sequence", annotations = mapOf("key1" to "value1")
)

val srangeList = listOf(
    SeqPositionRanges.of(recordA, 8..28),
    recordB.range(25..40),
    SeqPositionRanges.of(SeqPosition(recordA, 27), SeqPosition(recordA, 40)),
    SeqPositionRanges.of(recordB, 3..19)
)

println("Ranges in the list:")
srangeList.forEach { println("  $it") }

Ranges in the list:
  Sequence 1:8..Sequence 1:28
  Sequence 2:25..Sequence 2:40
  Sequence 1:27..Sequence 1:40
  Sequence 2:3..Sequence 2:19

val comparator: Comparator<SRange> = SeqRangeSort.by(SeqRangeSort.numberThenAlphaSort, SeqRangeSort.leftEdge)

println("Non-coalescing set:")
nonCoalescingSetOf(comparator, srangeList).forEach { println("  $it") }

println()
println("Coalescing set:")
coalescingSetOf(comparator, srangeList).forEach { println("  $it") }

Non-coalescing set:
  Sequence 1:8..Sequence 1:28
  Sequence 1:27..Sequence 1:40
  Sequence 2:3..Sequence 2:19
  Sequence 2:25..Sequence 2:40

Coalescing set:
  Sequence 1:8..Sequence 1:40
  Sequence 2:3..Sequence 2:19
  Sequence 2:25..Sequence 2:40

Reading BED files¶

Given a FASTA file and a BED file with coordinates relative to it, bedfileToSRangeSet builds an SRangeSet. From there the whole range API applies, and toDataFrame renders the result as a Kotlin DataFrame.

val fasta = "../src/test/kotlin/biokotlin/genome/chr9chr10short.fa"
val bedFile = "../src/test/kotlin/biokotlin/genome/chr9chr10_SHORTwithOverlaps.bed"

val srangeSet = bedfileToSRangeSet(bedFile, fasta)
println("Size of srangeSet: ${srangeSet.size}")
srangeSet.toDataFrame().print()

bedfileToSRangeSet: calling fastaToNucSeq
fastaToNucSeq: finished chrom chr9
fastaToNucSeq: finished chrom chr10
Size of srangeSet: 10
      ID start end    range
 0  chr9    21  24   21..24
 1  chr9    51  75   51..75
 2  chr9   151 225 151..225
 3  chr9   201 300 201..300
 4  chr9   341 375 341..375
 5  chr9   501 600 501..600
 6 chr10    71 109  71..109
 7 chr10   131 178 131..178
 8 chr10   176 300 176..300
 9 chr10   503 678 503..678

Intersecting a single range against a set¶

intersectingRanges finds every range in a set that overlaps one range of interest. A single-position range is a convenient way to ask "which features cover this SNP?".

// Reuse a record that is already in the set so the coordinates line up.
val seqRecT = srangeSet.elementAt(0).start.seqRecord

val snp = SeqPosition(seqRecT, 220)..SeqPosition(seqRecT, 220)
val snpHits = snp.intersectingRanges(srangeSet)
println("ranges covering position 220: $snpHits")
snpHits.toDataFrame().print()

ranges covering position 220: [chr9:151..chr9:225, chr9:201..chr9:300]
     ID start end    range
 0 chr9   151 225 151..225
 1 chr9   201 300 201..300

The same call works for a wider interval:

val window = SeqPosition(seqRecT, 70)..SeqPosition(seqRecT, 200)
val windowHits = window.intersectingRanges(srangeSet)
println("window: $window")
println("ranges overlapping it: $windowHits")
windowHits.toDataFrame().print()

window: chr9:70..chr9:200
ranges overlapping it: [chr9:51..chr9:75, chr9:151..chr9:225]
     ID start end    range
 0 chr9    51  75   51..75
 1 chr9   151 225 151..225