VcfParser.java

package org.pharmgkb.parser.vcf;

import java.io.BufferedReader;
import java.io.Closeable;
import java.io.IOException;
import java.lang.invoke.MethodHandles;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import com.google.common.base.Preconditions;
import org.apache.commons.lang3.StringUtils;
import org.jspecify.annotations.Nullable;
import org.pharmgkb.parser.vcf.model.BaseMetadata;
import org.pharmgkb.parser.vcf.model.ContigMetadata;
import org.pharmgkb.parser.vcf.model.FormatMetadata;
import org.pharmgkb.parser.vcf.model.IdDescriptionMetadata;
import org.pharmgkb.parser.vcf.model.InfoMetadata;
import org.pharmgkb.parser.vcf.model.VcfMetadata;
import org.pharmgkb.parser.vcf.model.VcfPosition;
import org.pharmgkb.parser.vcf.model.VcfSample;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;


/**
 * This class parses a VCF file.
 *
 * @author Mark Woon
 */
public class VcfParser implements Closeable {
  private static final Logger sf_logger = LoggerFactory.getLogger(MethodHandles.lookup().lookupClass());

  private static final char TAB = '\t';
  private static final char COMMA = ',';
  private static final char COLON = ':';
  private static final char SEMICOLON = ';';
  private static final String GLE = "GLE";
  // GLE is a VCF 4.1/4.2 FORMAT String whose own comma-separated genotype:likelihood pairs contain colons. The
  // first group consumes those pairs; the optional second group is the ordinary colon-separated remainder of the
  // sample field for FORMAT keys following GLE. Restrict each pair's genotype to allele indexes so a malformed,
  // arbitrary GLE string (which the parser deliberately preserves) is not mistaken for this special encoding.
  private static final Pattern GLE_VALUE_PATTERN = Pattern.compile(
      "((?:(?:\\d+|\\.)(?:[|/](?:\\d+|\\.))*:[^,:]+)(?:,(?:\\d+|\\.)(?:[|/](?:\\d+|\\.))*:[^,:]+)*)(?::(.*))?");
  // the mandatory fixed fields, in order; in the column-header line CHROM is written as "#CHROM"
  private static final List<String> FIXED_FIELD_NAMES =
      List.of("CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO");

  private final boolean m_rsidsOnly;
  private final BufferedReader m_reader;
  private @Nullable VcfMetadata m_vcfMetadata;
  private final VcfLineParser m_vcfLineParser;

  private int m_lineNumber;
  private boolean m_alreadyFinished;



  private VcfParser(BufferedReader reader, boolean rsidsOnly, VcfLineParser lineParser) {
    m_reader = reader;
    m_rsidsOnly = rsidsOnly;
    m_vcfLineParser = lineParser;
  }


  /**
   * Parses metadata only.
   * This method should be if only the metadata is needed; otherwise, {@link #parse()} is preferred.
   */
  public VcfMetadata parseMetadata() throws IOException {

    if (m_vcfMetadata != null) {
      throw new IllegalStateException("Metadata has already been parsed.");
    }
    VcfMetadata.Builder mdBuilder = new VcfMetadata.Builder();

    // the first physical line must be the single ##fileformat declaration
    String first = m_reader.readLine();
    m_lineNumber++;
    if (first == null || !first.startsWith("##fileformat=")) {
      throw new VcfFormatException("Not a VCF file: the first line must be ##fileformat", m_lineNumber);
    }
    parseMetadataLine(mdBuilder, first);

    String line;
    boolean foundHeader = false;
    while ((line = m_reader.readLine()) != null) {
      m_lineNumber++;
      if (line.startsWith("##")) {
        if (line.startsWith("##fileformat=")) {
          throw new VcfFormatException("Duplicate ##fileformat line", m_lineNumber);
        }
        parseMetadataLine(mdBuilder, line);
      } else if (line.startsWith("#")) {
        try {
          parseColumnInfo(mdBuilder, line);
        } catch (VcfFormatException ex) {
          ex.addMetadata(m_lineNumber, "column (# header)");
          throw ex;
        } catch (RuntimeException e) {
          throw new VcfFormatException(m_lineNumber, "column (# header)", e);
        }
        foundHeader = true;
        break;
      } else {
        throw new VcfFormatException("Unexpected line before the #CHROM header (only ## metadata lines are allowed here)",
            m_lineNumber);
      }
    }
    m_vcfMetadata = mdBuilder.build();
    if (!foundHeader) {
      throw new VcfFormatException("No column header line (starting with #CHROM) was found before the end of the file",
          m_lineNumber);
    }

    // check sample lists
    if (m_vcfMetadata.getNumSamples() == m_vcfMetadata.getSamples().size()) {
      for (int i = 0; i < m_vcfMetadata.getNumSamples(); i++) {
        String sampleName = m_vcfMetadata.getSampleName(i);
        if (!m_vcfMetadata.getSamples().containsKey(sampleName)) {
          sf_logger.warn("Sample {} is missing in the metadata", sampleName);
        }
      }
    } else {
      sf_logger.warn("There are {} samples in the header but {} in the metadata", m_vcfMetadata.getNumSamples(),
          m_vcfMetadata.getSamples().size());
    }

    // deliver the metadata to the line parser once, before any data lines (this method runs at most once)
    m_vcfLineParser.parseMetadata(m_vcfMetadata);

    return m_vcfMetadata;
  }


  /**
   * Gets VCF metadata (if it has already been parsed).
   */
  public @Nullable VcfMetadata getMetadata() {
    return m_vcfMetadata;
  }


  /**
   * Parses the entire VCF file (including the metadata).
   * <p>
   * This is the preferred way to read a VCF file.
   */
  public void parse() throws IOException {
    boolean hasNext = true;
    while (hasNext) {
      hasNext = parseNextLine();
    }
  }

  /**
   * Parses just the next data line available, also reading all the metadata if it has not been read.
   * This is a specialized method; in general calling {@link #parse()} to parse the entire stream is preferred.
   *
   * @return Whether another line may be available to read; false only if and only if this is the last line available
   * @throws IllegalStateException If the stream was already fully parsed
   * @throws VcfFormatException If a VCF formatting error is found
   */
  public boolean parseNextLine() throws IOException, VcfFormatException {

    if (m_alreadyFinished) {
      // prevents user errors from causing infinite loops
      throw new IllegalStateException("Already finished reading the stream");
    }

    if (m_vcfMetadata == null) {
      parseMetadata();
    }

    String line = m_reader.readLine();
    if (line == null) {
      m_alreadyFinished = true;
      return false;
    }
    m_lineNumber++;
    if (line.startsWith("#")) {
      // the VCF spec has no comment syntax: "#" is significant only for "##" metadata lines and the single column
      // header line, both of which are only valid before the first data line (and are already consumed by then)
      throw new VcfFormatException("Unexpected line starting with '#' after the column header; VCF has no comment "
          + "syntax", m_lineNumber);
    }

    try {
      if (StringUtils.stripToNull(line) == null) {
        throw new VcfFormatException("Empty line", m_lineNumber);
      }
      List<String> data = toList(TAB, line);
      if (data.size() != m_vcfMetadata.getNumColumns()) {
        throw new VcfFormatException("Data line does not have expected number of columns (got " + data.size() +
            " vs. " + m_vcfMetadata.getNumColumns() + ")", m_lineNumber);
      }

      // every fixed field is mandatory; an empty field is invalid (the missing value must be ".")
      for (int i = 0; i < FIXED_FIELD_NAMES.size(); i++) {
        if (data.get(i).isEmpty()) {
          throw new VcfFormatException(FIXED_FIELD_NAMES.get(i) + " field is empty; the missing value must be '.'",
              m_lineNumber);
        }
      }

      // CHROM
      String chromosome = data.get(0);

      // POS
      long position;
      try {
        position = Long.parseLong(data.get(1));
      } catch (NumberFormatException e) {
        throw new VcfFormatException("POS '" + data.get(1) + "' is not a number");
      }

      // ID
      List<String> ids = null;
      if (!data.get(2).equals(".")) {
        if (m_rsidsOnly && !VcfUtils.RSID_PATTERN.matcher(data.get(2)).find()) {
          return true;
        }
        ids = toList(SEMICOLON, data.get(2));
      } else if (m_rsidsOnly) {
        return true;
      }

      // REF
      String ref = data.get(3);

      // ALT
      List<String> alt = null;
      if (!data.get(4).equals(".")) {
        alt = toList(COMMA, data.get(4));
      }

      // FILTER
      List<String> filters = null;
      if (!data.get(6).equals("PASS")) {
        filters = toList(SEMICOLON, data.get(6));
      }

      // FORMAT
      List<String> format = null;
      if (data.size() >= 9) {
        format = toList(COLON, data.get(8));
      }

      // QUAL and INFO are parsed lazily by VcfPosition (see setRawQuality/setRawInfo); many consumers never read them.
      VcfPosition pos = new VcfPosition(chromosome, position, ids, ref, alt,
          null, filters, null, format);
      pos.setRawQuality(data.get(5));
      pos.setRawInfo(data.get(7));
      List<VcfSample> samples = new ArrayList<>();
      for (int x = 9; x < data.size(); x++) {
        List<String> values = toSampleValues(format, data.get(x));
        VcfUtils.fillEmptyEntriesWithDot(sf_logger, "sample value", values);
        // per the VCF spec, trailing FORMAT sub-fields may be dropped from a sample; pad any missing ones with the
        // missing value so the sample's value count matches the FORMAT key count
        if (format != null) {
          while (values.size() < format.size()) {
            values.add(".");
          }
        }
        samples.add(new VcfSample(format, values));
      }

      m_vcfLineParser.parseLine(m_vcfMetadata, pos, samples);
      return true;

    } catch (VcfFormatException ex) {
      ex.addMetadata(m_lineNumber, "data");
      throw ex;
    } catch (RuntimeException e) {
      throw new VcfFormatException(m_lineNumber, "data", e);
    }
  }

  @Override
  public void close() {
    try {
      m_reader.close();
    } catch (Exception ex) {
      sf_logger.info("Error closing reader", ex);
    }
  }


  /**
   * Parses a ## metadata line, annotating any failure with the line number.
   */
  private void parseMetadataLine(VcfMetadata.Builder mdBuilder, String line) {
    try {
      parseMetadata(mdBuilder, line);
    } catch (VcfFormatException ex) {
      ex.addMetadata(m_lineNumber, "metadata");
      throw ex;
    } catch (RuntimeException e) {
      throw new VcfFormatException(m_lineNumber, "metadata", e);
    }
  }

  /**
   * Parses a metadata line (starts with ##).
   */
  private void parseMetadata(VcfMetadata.Builder mdBuilder, String line) {

    int idx = line.indexOf("=");
    String propName = line.substring(2, idx).trim();
    String propValue = line.substring(idx + 1).trim();

    sf_logger.debug("{} : {}", propName, propValue);

    switch (propName) {
      case "fileformat":
        mdBuilder.setFileFormat(propValue);
        break;

      case "ALT":
      case "FILTER":
      case "INFO":
      case "FORMAT":
      case "contig":
      case "SAMPLE":
      case "PEDIGREE":
        parseMetadataProperty(mdBuilder, propName, removeAngleBrackets(propValue));
        break;

      case "assembly":
      case "pedigreeDB":
      default:
        mdBuilder.addRawProperty(propName, propValue);
    }
  }

  /**
   * Removes double quotation marks around a string.
   * @throws IllegalArgumentException If angle brackets are not present
   */
  private static String removeAngleBrackets(String string) throws VcfFormatException {
    if (string.startsWith("<") && string.endsWith(">")) {
      return string.substring(1, string.length() - 1);
    }
    throw new VcfFormatException("Angle brackets not present for: '" + string + "'");
  }

  /**
   * Converts metadata name-value pair into object.
   */
  private void parseMetadataProperty(VcfMetadata.Builder mdBuilder,
      String propName, String value) {
    Map<String, String> props = VcfUtils.extractPropertiesFromLine(value);
    switch (propName.toLowerCase()) {
      case "alt":
        mdBuilder.addAlt(new IdDescriptionMetadata(props, true));
        break;
      case "filter":
        mdBuilder.addFilter(new IdDescriptionMetadata(props, true));
        break;
      case "info":
        mdBuilder.addInfo(new InfoMetadata(props));
        break;
      case "format":
        mdBuilder.addFormat(new FormatMetadata(props));
        break;
      case "contig":
        mdBuilder.addContig(new ContigMetadata(props));
        break;
      case "sample":
        mdBuilder.addSample(new IdDescriptionMetadata(props, true));
        break;
      case "pedigree":
        mdBuilder.addPedigree(new BaseMetadata(props));
        break;
    }
  }

  /**
   * Splits {@code string} on the single character {@code delim}, matching
   * {@code Pattern.compile(String.valueOf(delim)).split(string, -1)}: leading, interior, and trailing empty fields
   * are all kept, and a string containing no delimiter yields a single-element list. Avoids the regex engine, which
   * is a meaningful per-line cost when parsing large VCFs.
   * <p>
   * Callers are responsible for handling any empty entries in the result (VCF does not allow zero-length fields); see
   * {@code EMPTY_FIELD_HANDLING.md} for how each field type handles them.
   */
  // package-private for differential testing against Pattern.split
  static List<String> toList(char delim, String string) {
    int idx = string.indexOf(delim);
    if (idx < 0) {
      List<String> single = new ArrayList<>(1);
      single.add(string);
      return single;
    }
    List<String> list = new ArrayList<>();
    int start = 0;
    while (idx >= 0) {
      list.add(string.substring(start, idx));
      start = idx + 1;
      idx = string.indexOf(delim, start);
    }
    list.add(string.substring(start));
    return list;
  }

  /**
   * Splits a sample field according to its FORMAT keys. GLE is the sole VCF 4.1/4.2 exception to ordinary
   * colon-delimited sample values: its value is a String containing comma-separated {@code genotype:likelihood}
   * pairs. See the GLE example in the VCF 4.2 specification.
   */
  private static List<String> toSampleValues(@Nullable List<String> format, String sample) {
    if (format == null) {
      return toList(COLON, sample);
    }
    int gleIndex = format.indexOf(GLE);
    if (gleIndex < 0) {
      return toList(COLON, sample);
    }

    List<String> values = new ArrayList<>();
    int start = 0;
    for (int i = 0; i < gleIndex; i++) {
      int separator = sample.indexOf(COLON, start);
      if (separator < 0) {
        // GLE and every following field were dropped; ordinary splitting lets the caller pad the trailing fields.
        return toList(COLON, sample);
      }
      values.add(sample.substring(start, separator));
      start = separator + 1;
    }

    String gleAndRemainder = sample.substring(start);
    Matcher matcher = GLE_VALUE_PATTERN.matcher(gleAndRemainder);
    if (!matcher.matches()) {
      // A missing GLE value (".") has no internal colon; ordinary splitting handles it and malformed values retain
      // the existing sample-arity validation behavior.
      return toList(COLON, sample);
    }
    values.add(matcher.group(1));
    if (matcher.group(2) != null) {
      values.addAll(toList(COLON, matcher.group(2)));
    }
    return values;
  }

  public int getLineNumber() {
    return m_lineNumber;
  }

  private void parseColumnInfo(VcfMetadata.Builder mdBuilder, String line) {
    List<String> cols = toList(TAB, line);
    if (cols.size() < 8) {
      throw new VcfFormatException("Header line does not have mandatory (tab-delimited) columns", m_lineNumber);
    }
    for (int i = 0; i < FIXED_FIELD_NAMES.size(); i++) {
      String expected = i == 0 ? "#" + FIXED_FIELD_NAMES.get(0) : FIXED_FIELD_NAMES.get(i);
      if (!cols.get(i).equals(expected)) {
        throw new VcfFormatException("Header column " + (i + 1) + " must be '" + expected +
            "' but was '" + cols.get(i) + "'", m_lineNumber);
      }
    }
    if (cols.size() > 8 && !cols.get(8).equals("FORMAT")) {
      throw new VcfFormatException("Header column 9 must be 'FORMAT' when sample columns are present but was '" +
          cols.get(8) + "'", m_lineNumber);
    }
    Set<String> sampleNames = new HashSet<>();
    for (int i = 9; i < cols.size(); i++) {
      if (!sampleNames.add(cols.get(i))) {
        throw new VcfFormatException("Duplicate sample name '" + cols.get(i) + "' in header", m_lineNumber);
      }
    }
    mdBuilder.setColumns(cols);
  }


  public static class Builder {
    private BufferedReader m_reader;
    private Path m_vcfFile;
    private boolean m_rsidsOnly;
    private VcfLineParser m_vcfLineParser;


    /**
     * Provides the {@link Path} to the VCF file to parse.
     * <p>
     * This reads the file as plain, uncompressed text; it does not accept a compressed VCF (e.g. {@code .vcf.gz}) and
     * does not attempt to decompress one. This is deliberate: {@link java.util.zip.GZIPInputStream} cannot reliably read
     * BGZF (block gzip), the format most bioinformatics tools use for compressed VCF, since a BGZF stream is a sequence
     * of independent gzip blocks with trailing extra data that plain {@code GZIPInputStream} does not expect and can
     * mishandle. To read a compressed VCF, decompress it with a library appropriate to how it was compressed (BGZF vs.
     * plain gzip) and pass the result to {@link #fromReader}.
     */
    public Builder fromFile(Path dataFile) {
      Preconditions.checkNotNull(dataFile);
      if (m_reader != null) {
        throw new IllegalStateException("Already loading from reader");
      }
      if (!dataFile.toString().endsWith(".vcf")) {
        throw new IllegalArgumentException("Not a VCF file (doesn't end with .vcf extension");
      }
      m_vcfFile = dataFile;
      return this;
    }

    /**
     * Provides a {@link BufferedReader} to the beginning of the VCF file to parse.
     */
    public Builder fromReader(BufferedReader reader) {
      Preconditions.checkNotNull(reader);
      if (m_vcfFile != null) {
        throw new IllegalStateException("Already loading from file");
      }
      m_reader = reader;
      return this;
    }

    /**
     * Tells parser to ignore data lines that are not associated with an RSID.
     */
    public Builder rsidsOnly() {
      m_rsidsOnly = true;
      return this;
    }

    public Builder parseWith(VcfLineParser lineParser) {
      Preconditions.checkNotNull(lineParser);
      m_vcfLineParser = lineParser;
      return this;
    }


    public VcfParser build() throws IOException {
      if (m_vcfLineParser == null) {
        throw new IllegalStateException("Missing VcfLineParser");
      }
      if (m_vcfFile != null) {
        m_reader = Files.newBufferedReader(m_vcfFile);
      }
      if (m_reader == null) {
        throw new IllegalStateException("Must specify either file or reader to parse");
      }
      return new VcfParser(m_reader, m_rsidsOnly, m_vcfLineParser);
    }
  }
}