MAF processing¶

Open this tutorial as a notebook

MAF alignment files can be summarised into per-base coverage and identity against the reference, and exported as BED or wiggle tracks for a genome browser.

How BioKotlin summarises coverage and identity from a MAF alignment

%use biokotlin

import biokotlin.genome.*
import org.jetbrains.kotlinx.dataframe.api.print
import org.jetbrains.kotlinx.dataframe.io.writeCSV
import java.io.File

// A small single-chromosome MAF from the test suite, and a scratch directory
// for the files produced below.
val mafFile = "../src/test/kotlin/biokotlin/testData/LineA.maf"
val mafDir = File(mafFile).parent
val outputDir = File(System.getProperty("java.io.tmpdir"), "biokotlin-maf-tutorial")
outputDir.mkdirs()
println("writing output to $outputDir")

writing output to /tmp/biokotlin-maf-tutorial

Coverage and identity across a region¶

getCoverageAndIdentityFromMAFs scans every .maf file in a directory in parallel and returns two arrays: coverage and identity at each reference position in the requested interval. With one alignment file the maximum value of either is 1.

val contig = "1"
val start = 13
val stop = 50

val coverageAndIdentity = GetCovIDFromMAFMultiThread()
    .getCoverageAndIdentityFromMAFs(contig, start, stop, mafDir)

getCoverageAndIdentityFromMAFs-MF: outside the loop, before combineResults: 0.026737939 seconds
getCoverageAndIdentityFromMAFs-MF: processing to closing result channel took: 0.283589679 seconds
combineResults-MF: finished  combineResults in : 0.253754871 seconds

Over a short interval it is easy enough to print the arrays directly:

println("position  coverage  identity")
for (idx in 0..(stop - start)) {
    println("${start + idx}\t${coverageAndIdentity.first[idx]}\t${coverageAndIdentity.second[idx]}")
}

position  coverage  identity
13  1   1
14  1   1
15  1   1
16  1   1
17  1   1
18  1   1
19  1   1
20  1   0
21  1   0
22  1   1
23  1   1
24  1   1
25  1   1
26  1   1
27  1   1
28  1   1
29  1   1
30  1   1
31  1   1
32  1   1
33  1   1
34  1   1
35  1   1
36  1   1
37  1   1
38  1   1
39  1   1
40  1   1
41  1   1
42  1   1
43  1   1
44  1   1
45  1   1
46  1   1
47  1   1
48  1   1
49  1   1
50  1   1

Exporting a BED file¶

createBedFileFromCoverageIdentity writes the intervals that meet a minimum coverage and identity. Set both thresholds to the number of aligned assemblies to require full support.

val minCov = 1
val minId = 1
val bedFile = File(outputDir, "coverage_${minCov}_and_${minId}.bed").path

createBedFileFromCoverageIdentity(
    coverageAndIdentity.first, coverageAndIdentity.second,
    contig, start, minCov, minId, bedFile
)
println(File(bedFile).readText())

createBedFileFromCoverageIdentity: File written to /tmp/biokotlin-maf-tutorial/coverage_1_and_1.bed
#gffTags
1   12  19  Name=1_12_19
1   21  50  Name=1_21_50

Exporting wiggle tracks¶

createWiggleFilesFromCoverageIdentity writes one wiggle file for coverage and one for identity, ready to load into IGV alongside a set of genic regions.

createWiggleFilesFromCoverageIdentity(
    coverageAndIdentity.first, coverageAndIdentity.second, contig, outputDir.path + "/"
)
outputDir.listFiles()?.sortedBy { it.name }?.forEach { println(it.name) }

createWiggleFilesFromCoverageIdentity: Identity written to /tmp/biokotlin-maf-tutorial//identity_1.wig
createWiggleFilesFromCoverageIdentity: Coverage written to /tmp/biokotlin-maf-tutorial//coverage_1.wig
coverage_1.wig
coverage_1_and_1.bed
identity_1.wig

Converting wiggle to BigWig¶

Wiggle files get large quickly. To load them into IGV comfortably, convert them to BigWig with the wigToBigWig utility, which needs a tab-delimited file of chromosome sizes with the chromosome name in the first column and its length in the second:

wigToBigWig input.wig b73NAM_chrom.sizes output.bw

Percent coverage and identity per chromosome¶

getCoverageIdentityPercentForMAF summarises a single MAF file, reporting how much of each reference contig the alignment covered and how much of that matched. The result is a DataFrame.

val covIdDF = getCoverageIdentityPercentForMAF(mafFile)
covIdDF!!.print()

getCoverageIdentityPercentForMAF: begin reading MAF file into blocks ...
getCoverageIdentifyPercentForMAF: time to read MAF file to blocks:  0.001825507 seconds
Processing the MAF blocks per-chrom
finished chrom 1
   contig numRegionBPs percentCov percentId
 0      1       130000  83.760769 80.618462

It writes out as CSV like any other DataFrame:

val csv = File(outputDir, "chrom_stats.csv").path
covIdDF!!.writeCSV(csv)
println(File(csv).readText())

contig,numRegionBPs,percentCov,percentId
1,130000,83.76076923076923,80.61846153846155

Pass a region in contig:start-end form to restrict the summary to part of a chromosome.

val regionDF = getCoverageIdentityPercentForMAF(mafFile, "1:9000-9300")
regionDF!!.print()

getCoverageIdentityPercentForMAF: begin reading MAF file into blocks ...
getCoverageIdentifyPercentForMAF: time to read MAF file to blocks:  0.002294465 seconds
Processing the MAF blocks per-chrom
finished chrom 1
   contig numRegionBPs percentCov percentId
 0      1          301  86.378738 63.787375