path: root/Master/texmf-dist/tex/lualatex/pgfmolbio/pgfmolbio.chromatogram.lua

--
-- This is file `pgfmolbio.chromatogram.lua',
-- generated with the docstrip utility.
--
-- The original source files were:
--
-- pgfmolbio.dtx  (with options: `pmb-chr-lua')
--
-- Copyright (C) 2013 by Wolfgang Skala
--
-- This work may be distributed and/or modified under the
-- conditions of the LaTeX Project Public License, either version 1.3
-- of this license or (at your option) any later version.
-- The latest version of this license is in
--   http://www.latex-project.org/lppl.txt
-- and version 1.3 or later is part of all distributions of LaTeX
-- version 2005/12/01 or later.
--
module("pgfmolbio.chromatogram", package.seeall)


if luatexbase then
  luatexbase.provides_module{
    name          = "pgfmolbio.chromatogram",
    version       = 0.2,
    date          = "2012/10/01",
    description   = "DNA sequencing chromatograms",
    author        = "Wolfgang Skala",
    copyright     = "Wolfgang Skala",
    license       = "LPPL",
  }
end

local ALL_BASES = {"A", "C", "G", "T"}
local PGFKEYS_PATH = "/pgfmolbio/chromatogram/"

local stringToDim = pgfmolbio.stringToDim
local dimToString = pgfmolbio.dimToString
local packageError = pgfmolbio.packageError
local packageWarning = pgfmolbio.packageWarning
local getRange = pgfmolbio.getRange

local function stdProbStyle(prob)
  local color = ""
  if prob >= 0 and prob < 10 then
    color = "black"
  elseif prob >= 10 and prob < 20 then
    color = "pmbTraceRed"
  elseif prob >= 20 and prob < 30 then
    color = "pmbTraceYellow"
  else
    color = "pmbTraceGreen"
  end
  return "ultra thick, " .. color
end

local function findBasesInStr(target)
  if not target then return end
  local result = {}
  for _, v in ipairs(ALL_BASES) do
    if target:upper():find(v) then
      table.insert(result, v)
    end
  end
  return result
end

local function readInt(file, n, offset)
  if offset then file:seek("set", offset) end
  local result = 0
  for i = 1, n do
    result = result * 0x100 + file:read(1):byte()
  end
  return result
end

Chromatogram = {}

function Chromatogram:new()
  newChromatogram = {
    sampleMin = 1,
    sampleMax = 500,
    sampleStep = 1,
    peakMin = -1,
    peakMax = -1,
    xUnit = stringToDim("0.2mm"),
    yUnit = stringToDim("0.01mm"),
    samplesPerLine = 500,
    baselineSkip = stringToDim("3cm"),
    canvasHeight = stringToDim("2cm"),
    traceStyle = {
      A = PGFKEYS_PATH .. "trace A style",
      C = PGFKEYS_PATH .. "trace C style",
      G = PGFKEYS_PATH .. "trace G style",
      T = PGFKEYS_PATH .. "trace T style"
    },
    tickStyle = {
      A = PGFKEYS_PATH .. "tick A style",
      C = PGFKEYS_PATH .. "tick C style",
      G = PGFKEYS_PATH .. "tick G style",
      T = PGFKEYS_PATH .. "tick T style"
    },
    tickLength = stringToDim("1mm"),
    baseLabelText = {
      A = "\\pgfkeysvalueof{" .. PGFKEYS_PATH .. "base label A text}",
      C = "\\pgfkeysvalueof{" .. PGFKEYS_PATH .. "base label C text}",
      G = "\\pgfkeysvalueof{" .. PGFKEYS_PATH .. "base label G text}",
      T = "\\pgfkeysvalueof{" .. PGFKEYS_PATH .. "base label T text}"
    },
    baseLabelStyle = {
      A = PGFKEYS_PATH .. "base label A style",
      C = PGFKEYS_PATH .. "base label C style",
      G = PGFKEYS_PATH .. "base label G style",
      T = PGFKEYS_PATH .. "base label T style"
    },
    showBaseNumbers = true,
    baseNumberMin = -1,
    baseNumberMax = -1,
    baseNumberStep = 10,
    probDistance = stringToDim("0.8cm"),
    probStyle = stdProbStyle,
    tracesDrawn = ALL_BASES,
    ticksDrawn = "ACGT",
    baseLabelsDrawn = "ACGT",
    probabilitiesDrawn = "ACGT",
  }
  setmetatable(newChromatogram, self)
  self.__index = self
  return newChromatogram
end

function Chromatogram:getMinMaxProbability()
  local minProb = 0
  local maxProb = 0
  for _, currPeak in ipairs(self.selectedPeaks) do
    for __, currProb in pairs(currPeak.prob) do
      if currProb > maxProb then maxProb = currProb end
      if currProb < minProb then minProb = currProb end
    end
  end
  return minProb, maxProb
end

function Chromatogram:getSampleAndPeakIndex(baseIndex, isLowerLimit)
  local sampleId, peakId

  sampleId = tonumber(baseIndex)
  if sampleId then
    for i, v in ipairs(self.peaks) do
      if isLowerLimit then
        if v.offset >= sampleId then
          peakId = i
          break
        end
      else
        if v.offset == sampleId then
          peakId = i
          break
        elseif v.offset > sampleId then
          peakId = i - 1
          break
        end
      end
    end
  else
    peakId = tonumber(baseIndex:match("base%s*(%d+)"))
    if peakId then
      sampleId = self.peaks[peakId].offset
    end
  end
  return sampleId, peakId
end

function Chromatogram:readScfFile(filename)
  if filename ~= self.lastScfFile then
    self.lastScfFile = filename
    local scfFile, errorMsg = io.open(filename, "rb")
    if not scfFile then packageError(errorMsg) end

    self.samples = {A = {}, C = {}, G = {}, T = {}}
    self.peaks = {}
    self.header = {
      magicNumber = readInt(scfFile, 4, 0),
      samplesNumber = readInt(scfFile, 4),
      samplesOffset = readInt(scfFile, 4),
      basesNumber = readInt(scfFile, 4),
      leftClip = readInt(scfFile, 4),
      rightClip = readInt(scfFile, 4),
      basesOffset = readInt(scfFile, 4),
      comments = readInt(scfFile, 4),
      commentsOffset = readInt(scfFile, 4),
      version = readInt(scfFile, 4),
      sampleSize = readInt(scfFile, 4),
      codeSet = readInt(scfFile, 4),
      privateSize = readInt(scfFile, 4),
      privateOffset = readInt(scfFile, 4)
    }
    if self.header.magicNumber ~= 0x2E736366 then
      packageError(
        "Magic number in scf scfFile '" ..
        self.lastScfFile ..
        "' corrupt!"
      )
    end
    if self.header.version ~= 0x332E3030 then
      packageError(
        "Scf scfFile '" ..
        self.lastScfFile ..
        "' is not version 3.00!"
      )
    end
    scfFile:seek("set", self.header.samplesOffset)
    for baseIndex, baseName in ipairs(ALL_BASES) do
      for i = 1, self.header.samplesNumber do
        self.samples[baseName][i] =
          readInt(scfFile, self.header.sampleSize)
      end

      for _ = 1, 2 do
        local preValue = 0
        for i = 1, self.header.samplesNumber do
          self.samples[baseName][i] = self.samples[baseName][i] + preValue
          if self.samples[baseName][i] > 0xFFFF then
            self.samples[baseName][i] = self.samples[baseName][i] - 0x10000
          end
          preValue = self.samples[baseName][i]
        end
      end
    end
    for i = 1, self.header.basesNumber do
      self.peaks[i] = {
        offset = readInt(scfFile, 4),
        prob = {A, C, G, T},
        base
      }
    end

    for i = 1, self.header.basesNumber do
      self.peaks[i].prob.A = readInt(scfFile, 1)
    end

    for i = 1, self.header.basesNumber do
      self.peaks[i].prob.C = readInt(scfFile, 1)
    end

    for i = 1, self.header.basesNumber do
      self.peaks[i].prob.G = readInt(scfFile, 1)
    end

    for i = 1, self.header.basesNumber do
      self.peaks[i].prob.T = readInt(scfFile, 1)
    end

    for i = 1, self.header.basesNumber do
      self.peaks[i].base = string.char(readInt(scfFile, 1))
    end

    scfFile:close()
  end
end

function Chromatogram:setParameters(newParms)
  local keyHash = {
    sampleRange = function(v)
      local sampleRangeMin, sampleRangeMax, sampleRangeStep =
        getRange(
          v:trim(),
          "^([base]*%s*%d+)%s*%-",
          "%-%s*([base]*%s*%d+)",
          "step%s*(%d+)$"
        )
      self.sampleMin, self.peakMin =
        self:getSampleAndPeakIndex(sampleRangeMin, true)
      self.sampleMax, self.peakMax =
        self:getSampleAndPeakIndex(sampleRangeMax, false)
      if self.sampleMin >= self.sampleMax then
        packageError("Sample range is smaller than 1.")
      end
      self.sampleStep = sampleRangeStep or self.sampleStep
    end,
    xUnit = stringToDim,
    yUnit = stringToDim,
    samplesPerLine = tonumber,
    baselineSkip = stringToDim,
    canvasHeight = stringToDim,
    tickLength = stringToDim,
    showBaseNumbers = function(v)
      if v == "true" then return true else return false end
    end,
    baseNumberRange = function(v)
      local baseNumberRangeMin, baseNumberRangeMax, baseNumberRangeStep =
        getRange(
          v:trim(),
          "^([auto%d]*)%s+%-",
          "%-%s+([auto%d]*$)"
        )
      if tonumber(baseNumberRangeMin) then
        self.baseNumberMin = tonumber(baseNumberRangeMin)
      else
        self.baseNumberMin = self.peakMin
      end
      if tonumber(baseNumberRangeMax) then
        self.baseNumberMax = tonumber(baseNumberRangeMax)
      else
        self.baseNumberMax = self.peakMax
      end
      if self.baseNumberMin >= self.baseNumberMax then
        packageError("Base number range is smaller than 1.")
      end
      if self.baseNumberMin < self.peakMin then
        self.baseNumberMin = self.peakMin
        packageWarning("Lower base number range is smaller than lower sample range. It was adjusted to " .. self.baseNumberMin .. ".")
      end
      if self.baseNumberMax > self.peakMax then
        self.baseNumberMax = self.peakMax
        packageWarning("Upper base number range exceeds upper sample range. It was adjusted to " .. self.baseNumberMax .. ".")
      end
      self.baseNumberStep = tonumber(baseNumberRangeStep)
        or self.baseNumberStep
    end,
    probDistance = stringToDim,
    probStyle = function(v) return v end,
    tracesDrawn = findBasesInStr,
    ticksDrawn = function(v) return v end,
    baseLabelsDrawn = function(v) return v end,
    probabilitiesDrawn = function(v) return v end,
    probStyle = function(v) return v end
  }
  for key, value in pairs(newParms) do
    if keyHash[key] then
      self[key] = keyHash[key](value)
    end
  end
end

function Chromatogram:printTikzChromatogram()
  if pgfmolbio.errorCatched then return end
  self.selectedPeaks = {}
  local tIndex = 1
  for rPeakIndex, currPeak in ipairs(self.peaks) do
    if currPeak.offset >= self.sampleMin
        and currPeak.offset <= self.sampleMax then
      self.selectedPeaks[tIndex] = {
        offset = currPeak.offset + 1 - self.sampleMin,
        base = currPeak.base,
        prob = currPeak.prob,
        baseIndex = rPeakIndex,
        probXRight = self.sampleMax + 1 - self.sampleMin
      }
      if tIndex > 1 then
        self.selectedPeaks[tIndex-1].probXRight =
          (self.selectedPeaks[tIndex-1].offset
          + self.selectedPeaks[tIndex].offset) / 2
      end
      tIndex = tIndex + 1
    end
  end

  if tIndex > 1 then
    if self.baseNumberMin == -1 then
      self.baseNumberMin = self.selectedPeaks[1].baseIndex
    end
    if self.baseNumberMax == -1 then
      self.baseNumberMax = self.selectedPeaks[tIndex-1].baseIndex
    end
  end

  local samplesLeft = self.sampleMax - self.sampleMin + 1
  local currLine = 0
  while samplesLeft > 0 do
    local yLower = -currLine * self.baselineSkip
    local yUpper = -currLine * self.baselineSkip + self.canvasHeight
    local xRight =
      (math.min(self.samplesPerLine, samplesLeft) - 1) * self.xUnit
    tex.sprint(
      "\n\t\\draw [" .. PGFKEYS_PATH .. "canvas style] (" ..
      dimToString(0) ..
      ", " ..
      dimToString(yLower) ..
      ") rectangle (" ..
      dimToString(xRight) ..
      ", " ..
      dimToString(yUpper) ..
      ");"
    )
    samplesLeft = samplesLeft - self.samplesPerLine
    currLine = currLine + 1
  end

  for _, baseName in ipairs(self.tracesDrawn) do
    tex.sprint("\n\t\\draw [" .. self.traceStyle[baseName] .. "] ")
    local currSampleIndex = self.sampleMin
    local sampleX = 1
    local x = 0
    local y = 0
    local currLine = 0
    local firstPointInLine = true

    while currSampleIndex <= self.sampleMax do
      x = ((sampleX - 1) % self.samplesPerLine) * self.xUnit
      y = self.samples[baseName][currSampleIndex] * self.yUnit
        - currLine * self.baselineSkip
      if sampleX % self.sampleStep == 0 then
        if not firstPointInLine then
          tex.sprint(" -- ")
        else
          firstPointInLine = false
        end
        tex.sprint(
          "(" ..
          dimToString(x) ..
          ", " ..
          dimToString(y) ..
          ")"
        )
      end
      if sampleX ~= self.sampleMax + 1 - self.sampleMin then
        if sampleX >= (currLine + 1) * self.samplesPerLine then
          currLine = currLine + 1
          tex.sprint(";\n\t\\draw [" .. self.traceStyle[baseName] .. "] ")
          firstPointInLine = true
        end
      else
        tex.sprint(";")
      end
    sampleX = sampleX + 1
    currSampleIndex = currSampleIndex + 1
    end
  end

  local currLine = 0
  local lastProbX = 1
  local probRemainder = false

  for _, currPeak in ipairs(self.selectedPeaks) do
    while currPeak.offset > (currLine + 1) * self.samplesPerLine do
      currLine = currLine + 1
    end

    local x = ((currPeak.offset - 1) % self.samplesPerLine) * self.xUnit
    local yUpper = -currLine * self.baselineSkip
    local yLower = -currLine * self.baselineSkip - self.tickLength
    local tickOperation = ""
    if self.ticksDrawn:upper():find(currPeak.base) then
      tickOperation = "--"
    end

    tex.sprint(
      "\n\t\\draw [" ..
      self.tickStyle[currPeak.base] ..
      "] (" ..
      dimToString(x) ..
      ", " ..
      dimToString(yUpper) ..
      ") " ..
      tickOperation ..
      " (" ..
      dimToString(x) ..
      ", " ..
      dimToString(yLower) ..
      ")"
    )
    if self.baseLabelsDrawn:upper():find(currPeak.base) then
      tex.sprint(
        " node [" ..
        self.baseLabelStyle[currPeak.base] ..
        "] {" ..
        self.baseLabelText[currPeak.base] ..
        "}"
      )
    end

    if self.showBaseNumbers
        and currPeak.baseIndex >= self.baseNumberMin
        and currPeak.baseIndex <= self.baseNumberMax
        and (currPeak.baseIndex - self.baseNumberMin)
          % self.baseNumberStep == 0 then
      tex.sprint(
        " node [" ..
        PGFKEYS_PATH ..
        "base number style] {\\strut " ..
        currPeak.baseIndex ..
        "}"
      )
    end
    tex.sprint(";")

    if probRemainder then
      tex.sprint(probRemainder)
      probRemainder = false
    end
    local drawCurrProb =
      self.probabilitiesDrawn:upper():find(currPeak.base)
    local xLeft = lastProbX - 1 - currLine * self.samplesPerLine
    if xLeft < 0 then
      local xLeftPrev = (self.samplesPerLine + xLeft) * self.xUnit
      local xRightPrev = (self.samplesPerLine - 1) * self.xUnit
      local yPrev = -(currLine-1) * self.baselineSkip - self.probDistance
      if drawCurrProb then
        tex.sprint(
          "\n\t\\draw [" ..
          self.probStyle(currPeak.prob[currPeak.base]) ..
          "] (" ..
          dimToString(xLeftPrev) ..
          ", " ..
          dimToString(yPrev) ..
          ") -- (" ..
          dimToString(xRightPrev) ..
          ", " ..
          dimToString(yPrev) ..
          ");"
        )
      end
      xLeft = 0
    else
      xLeft = xLeft * self.xUnit
    end

    local xRight = currPeak.probXRight - 1 - currLine * self.samplesPerLine
    if xRight >= self.samplesPerLine then
      if drawCurrProb then
        local xRightNext = (xRight - self.samplesPerLine) * self.xUnit
        local yNext = -(currLine+1) * self.baselineSkip - self.probDistance
        probRemainder =
          "\n\t\\draw [" ..
          self.probStyle(currPeak.prob[currPeak.base]) ..
          "] (" ..
          dimToString(0) ..
          ", " ..
          dimToString(yNext) ..
          ") -- (" ..
          dimToString(xRightNext) ..
          ", " ..
          dimToString(yNext) ..
          ");"
      end
      xRight = (self.samplesPerLine - 1) * self.xUnit
    else
      xRight = xRight * self.xUnit
    end

    local y = -currLine * self.baselineSkip - self.probDistance
    if drawCurrProb then
      tex.sprint(
        "\n\t\\draw [" ..
        self.probStyle(currPeak.prob[currPeak.base]) ..
        "] (" ..
        dimToString(xLeft) ..
        ", " ..
        dimToString(y) ..
        ") -- (" ..
        dimToString(xRight) ..
        ", " ..
        dimToString(y) ..
        ");"
      )
    end
    lastProbX = currPeak.probXRight
  end
end
--
-- End of file `pgfmolbio.chromatogram.lua'.