tool_intersect.cpp.bk

// ******************************************************
// vcfCTools (c) 2011 Alistair Ward
// Marth Lab, Department of Biology, Boston College
// All rights reserved.
// ------------------------------------------------------
// Last modified: 18 February 2011
// ------------------------------------------------------
// Calculate the intersection of two vcf files or a vcf
// file and a bed file.
// ******************************************************

#include "tool_intersect.h"

using namespace std;
using namespace vcfCTools;

// intersectTool imlementation.
intersectTool::intersectTool(void)
  : AbstractTool()
{
  recordsInMemory = 100;
  passFilters = false;
  findCommon = false;
  findUnion  = false;
  findUnique = false;
  currentReferenceSequence = "";
}

// Destructor.
intersectTool::~intersectTool(void) {}

// Help
int intersectTool::Help(void) {
  cout << "Intersect help" << endl;
  cout << "Usage: ./vcfCTools intersect [options]." << endl;
  cout << endl;
  cout << "Options:" << endl;
  cout << "  -h, --help" << endl;
  cout << "	display intersect help." << endl;
  cout << "  -i, --in" << endl;
  cout << "	input vcf files (two, or one if intersecting with bed file)." << endl;
  cout << "  -b, --bed" << endl;
  cout << "	input bed file." << endl;
  cout << "  -o, --output" << endl;
  cout << "	output vcf file." << endl;
  cout << "  -c, --common" << endl;
  cout << "	output variants present in both files." << endl;
  cout << "  -u, --union" << endl;
  cout << "	output variants present in either file." << endl;
  cout << "  -q, --unique" << endl;
  cout << "	output variants unique to one of the files." << endl;
  cout << "  -p, --pass-filters" << endl;
  cout << "	Only variants that pass filters are considered." << endl;
  cout << "  -1, --snps" << endl;
  cout << "     analyse SNPs." << endl;
  cout << "  -2, --mnps" << endl;
  cout << "     analyse MNPs." << endl;
  cout << "  -3, --indels" << endl;
  cout << "     analyse indels." << endl;
  cout << endl;
  cout << "Additional information:" << endl;
  cout << "  The -c, -u and -q options require either 'a', 'b' or 'q' (not valid for -q) as an argument." << endl;
  cout << endl;
  cout << "  a: Write out records from the first file." << endl;
  cout << "  b: Write out records from the second file." << endl;
  cout << "  q: Write out records with the highest variant quality." << endl;
  cout << endl;
  exit(0);

  return 0;
}

// Parse the command line and get all required and optional arguments.
int intersectTool::parseCommandLine(int argc, char* argv[]) {
  commandLine = argv[0];
  for (int i = 2; i < argc; i++) {
    commandLine += " ";
    commandLine += argv[i];
  }

  int argument; // Counter for getopt.

  // Define the long options.
  static struct option long_options[] = {
    {"help", no_argument, 0, 'h'},
    {"bed", required_argument, 0, 'b'},
    {"in", required_argument, 0, 'i'},
    {"out", required_argument, 0, 'o'},
    {"pass-filters", no_argument, 0, 'p'},
    {"common", required_argument, 0, 'c'},
    {"union", required_argument, 0, 'u'},
    {"unique", required_argument, 0, 'q'},
    {"snps", no_argument, 0, '1'},
    {"mnps", no_argument, 0, '2'},
    {"indels", no_argument, 0, '3'},

    {0, 0, 0, 0}
  };

  while (true) {
    int option_index = 0;
    argument = getopt_long(argc, argv, "hb:i:o:pc:u:q:123", long_options, &option_index);

    if (argument == -1) {break;}
    switch (argument) {
      // Input bed file.
      case 'b':
        bedFile = optarg;
        break;

      // Input vcf file - required input.
      case 'i':
        vcfFiles.push_back(optarg);
        break;

      // Help.
      case 'h':
        return Help();

      // Output file.
      case 'o':
        outputFile = optarg;
        break;

      // Common variants.
      case 'c':
        findCommon = true;
        writeFrom = optarg;
        break;

      // Find the union.
      case 'u':
        findUnion = true;
        writeFrom = optarg;
        break;

      // Unique variants.
      case 'q':
        findUnique = true;
        writeFrom = optarg;
        break;

      // Only consider variants if they pass filters.
      case 'p':
        passFilters = true;
        break;

      // Analyse SNPs.
      case '1':
        processSnps = true;
        break;

      // Analyse MNPs.
      case '2':
        processMnps = true;
        break;

      // Analyse indels.
      case '3':
        processIndels = true;
        break;
      
      //
      case '?':
        cout << "Unknown option: " << argv[optind - 1] << endl;
        exit(1);
 
      // default
      default:
        abort ();
    }
  }

// Remaining arguments are unknown, so terminate with an error.
  if (optind < argc - 1) {
    cerr << "Unknown options." << endl;
    exit(1);
  }

// Check that either two vcf files or one vcf and one bed file is specified.
  if (vcfFiles.size() == 0 || (vcfFiles.size() == 1 and bedFile == "") || (bedFile != "" && vcfFiles.size() != 1) || vcfFiles.size() > 2) {
    cerr << "Two vcf files or a vcf and a bed file must be specified (--in, -i, --bed, -b)." << endl;
    exit(1);
  }

// Check whether finding common, union or unique variants.  If more than one of
// these options are selected, terminate the program.
  if ( (findCommon + findUnion + findUnique) > 1) {
    cerr << "Only one operation (-c [--common], -u [--union] or -q [--unique]) can be performed at once." << endl;
    exit(1);
  } else if (!findCommon && !findUnion && !findUnique) {
    cerr << "One operation (-c [--common], -u [--union] or -q [--unique]) must be selected." << endl;
    exit(1);
  } else {
    if (writeFrom == "a") {cerr << "Writing out records from file: " << vcfFiles[0] << endl;}
    else if (writeFrom == "b") {cerr << "Writing out records from file: " << vcfFiles[1] << endl;}
    else if ( (findCommon || findUnion) && writeFrom == "q") {cerr << "Writing out records with the highest quality value." << endl;}
    else {
      cerr << "The file from which the records are to be written needs to be selected." << endl;
      cerr << "    a - write records from the first inputted file." << endl;
      cerr << "    b - write records from the second inputted file." << endl;
      if (findCommon) {cerr << "    q - write records with the highest variant quality." << endl;}
      exit(1);
    }
  }
    
  return 0;
}

// Intersect two vcf files.  Intersect by variant position only.
void intersectTool::intersectVcf(vcf& v1, vcf& v2) {
  v1.success = v1.getRecord(currentReferenceSequence);
  v2.success = v2.getRecord(currentReferenceSequence);
  v1.update = true;
  v2.update = true;

// Define writeWhileParse.  If parsing through one of the vcf files to find the next
// common coordinate, this Boolean will instruct the parseVcf routine on whether or
// not to write the records to the output file.
  bool write;

  currentReferenceSequence = v1.variantRecord.referenceSequence;
  write = ( (findUnique && writeFrom == "a") || findUnion) ? true : false;
  v1.success = v1.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
  write = ( (findUnique && writeFrom == "b") || findUnion) ? true : false;
  v2.success = v2.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);

// Set the iterators to the start of each structure.
  v1.variantsIter = v1.variants.begin();
  v2.variantsIter = v2.variants.begin();

// Compare the variant structures for the two vcf files.  As variants are deleted,
// new ones are added until there are no variants left.  Terminate the comparisons
// as soon as the end of a file is reached.
  while (v1.variants.size() != 0 && v2.variants.size() != 0) {
    if (v1.variantsIter->first == v2.variantsIter->first) {
      if (!findUnique) {v1.writeRecord(output);}

      // Erase the variants from the structure.
      v1.variants.erase(v1.variantsIter);
      v2.variants.erase(v2.variantsIter);

      // Add the next variants to the structure if any more exist for this reference
      // sequence.
      if (v1.update && v1.success) {
        v1.addVariantToStructure();
        v1.success = v1.getRecord(currentReferenceSequence);
        if (v1.variantRecord.referenceSequence != currentReferenceSequence || !v1.success) {v1.update = false;}
      }
      if (v2.update && v2.success) {
        v2.addVariantToStructure();
        v2.success = v2.getRecord(currentReferenceSequence);
        if (v2.variantRecord.referenceSequence != currentReferenceSequence || !v2.success) {v2.update = false;}
      }

    // If the variant in v1 is after that in v2, remove the variant from v2 and update
    // it's structure.
    } else if (v1.variantsIter->first > v2.variantsIter->first) {
      if ( (findUnique && writeFrom == "b") || findUnion) {v2.writeRecord(output);}
      v2.variants.erase(v2.variantsIter);
      if (v2.update) {
        v2.addVariantToStructure();
        v2.success = v2.getRecord(currentReferenceSequence);
        if (v2.variantRecord.referenceSequence != currentReferenceSequence || !v2.success) {v2.update = false;}
      }

    // If the variant in v1 is after that in v2, remove the variant from v2 and update
    // it's structure.
    } else if (v2.variantsIter->first > v1.variantsIter->first) {
      if ( (findUnique && writeFrom == "a") || findUnion) {v1.writeRecord(output);}
      v1.variants.erase(v1.variantsIter);
      if (v1.update) {
        v1.addVariantToStructure();
        v1.success = v1.getRecord(currentReferenceSequence);
        if (v1.variantRecord.referenceSequence != currentReferenceSequence || !v1.success) {v1.update = false;}
      }
    }

    // If either of the variant structures is empty, but the end of the file hasn't been
    // reached, process the variants remaining in the other vcf files variant structure,
    // then rebuild the structures for the next reference sequence.
    //
    // If v1 is empty or both v1 and v2 are empty, the first conditional in the if
    // statement is true.  parseRemainingRef will then clear the contents of v2 if it 
    // wasn't empty.  If v1 isn't empy, but v2 is, the else condition is triggered.
    // Whichever way, if either v1 or v2 are empty, both of them will be empty and then
    // rebuilt after the if statement.
    if (v1.variants.size() == 0 && v1.success) {
      string fileId = "b";
      parseRemainingRef(v2, currentReferenceSequence, fileId);
      currentReferenceSequence = v1.variantRecord.referenceSequence;
      v1.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
      v2.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
    } else if (v2.variants.size() == 0 && v2.success) {
      string fileId = "a";
      parseRemainingRef(v1, currentReferenceSequence, fileId);
      currentReferenceSequence = v1.variantRecord.referenceSequence;
      v1.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
      v2.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
    }

    // Reset the iterators to the first elements in the structures.
    v1.variantsIter = v1.variants.begin();
    v2.variantsIter = v2.variants.begin();
  }

  // Check that both variant structures are empty.  At most, one may still have variants in it.
  if (v1.variants.size() != 0 && v2.variants.size() != 0) {
    cerr << "The intersection loop has terminated and yet both variant" << endl;
    cerr << "structures still contain entries.  A bug must exist in the" << endl;
    cerr << "intersection algorithm." << endl;
    cerr << "PROGRAM TERMINATED UNEXPECTEDLY." << endl;
    exit(1);
  }

  // If v1 still contains variants, process them.
  if (v1.variants.size() != 0) {

    // Empty the variant structure.
    string fileId = "a";
    parseRemainingRef(v1, currentReferenceSequence, fileId);

    // Parse any remaining variants in the vcf file.  These are all unique, since the second
    // vcf file has given up all it's variants
    while (v1.success) {
      currentReferenceSequence = v1.variantRecord.referenceSequence;
      v1.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
      parseRemainingRef(v1, currentReferenceSequence, fileId);
    }

  // If v2 still contains variants, process them.
  } else if (v2.variants.size() != 0) {

    // Empty the variant structure.
    string fileId = "a";
    parseRemainingRef(v1, currentReferenceSequence, fileId);

    // Parse any remaining variants in the vcf file.  These are all unique, since the second
    // vcf file has given up all it's variants
    while (v1.success) {
      currentReferenceSequence = v1.variantRecord.referenceSequence;
      v1.buildVariantStructure(recordsInMemory, currentReferenceSequence, write, output);
      parseRemainingRef(v1, currentReferenceSequence, fileId);
    }
  }
}

// Parse through the rest of the variants in the vcf file for a particular reference sequence,
// then process all variants remaining in the variant structure.
void intersectTool::parseRemainingRef(vcf& v, string& currentReferenceSequence, string& fileId) {
  while (v.success && v.variantRecord.referenceSequence == currentReferenceSequence) {
    v.addVariantToStructure();
    v.variantsIter = v.variants.begin();
    if ( (findUnique && writeFrom == fileId) || findUnion) {v.writeRecord(output);}
    v.variants.erase(v.variantsIter);
    v.success = v.getRecord(currentReferenceSequence);
  }

  // Process the variants remaining in the structure.
  for (v.variantsIter = v.variants.begin(); v.variantsIter != v.variants.end(); v.variantsIter++) {
    if ( (findUnique && writeFrom == fileId) || findUnion) {v.writeRecord(output);}
    v.variants.erase(v.variantsIter);
  }
}

// Intersect a vcf file and a bed file.  It is assumed that the 
// two files are sorted by genomic coordinates and the reference
// sequences are in the same order.  Do not group together variants
// in common reference sequence.
void intersectTool::intersectVcfBed(vcf& v, bed& b) {
  bool successBed = b.getRecord();
  bool successVcf = v.getRecord(currentReferenceSequence);
  string currentReferenceSequence = v.referenceSequence;
  v.update = true;

// As soon as the end of the first file is reached, there are no
// more intersections and the program can terminate.
  while (successVcf && successBed) {
    if (v.referenceSequence == b.referenceSequence) {
      if (v.position < b.start) {successVcf = v.parseVcf(b.referenceSequence, b.start, false, output, passFilters);}
      else if (v.position > b.end) {successBed = b.parseBed(v.referenceSequence, v.position);}
      else {
        *output << v.record << endl;
        successVcf = v.getRecord(currentReferenceSequence);
      }
    } else {
      if (v.referenceSequence == currentReferenceSequence) {successVcf = v.parseVcf(b.referenceSequence, b.start, false, output, passFilters);}
      if (b.referenceSequence == currentReferenceSequence) {successBed = b.parseBed(v.referenceSequence, v.position);}
      currentReferenceSequence = v.referenceSequence;
    }
  }
}

// Run the tool.
int intersectTool::Run(int argc, char* argv[]) {
  int getOptions = intersectTool::parseCommandLine(argc, argv);
  output = openOutputFile(outputFile);

// If intersection is between a vcf file and a bed file, create a vcf and a bed object
// and intersect.
  if (bedFile != "") {
    cerr << "Not updated bed intersection tool." << endl;
    exit(1);
    //vcf v; // Create a vcf object.
    //bed b; // Create a bed object.

    //v.openVcf(vcfFiles[0]);
    //b.openBed(bedFile);
    //v.parseHeader();

// Write the header to the output file.
    //string taskDescription = "##vcfCtools=intersect " + vcfFiles[0] + ", " + bedFile;
    //writeHeader(output, v, false, taskDescription);

// Intersect the files.
    //if (groupVariants) {
      //intersectVariantGroupsBed(v, b, output);
      //intersectVariantGroupsBed(v, b);
    //} else {
      //intersectVcfBed(v, b, output);
      //intersectVcfBed(v, b);
    //}

// Check that the input files had the same list of reference sequences.
// If not, it is possible that there were some problems.
    //checkReferenceSequences(v.referenceSequenceVector, b.referenceSequenceVector); // tools.cpp

// Close the vcf file and return.
    //v.closeVcf();
    //b.closeBed();
  } else {
    vcf v1; // Create a vcf object.
    vcf v2; // Create a vcf object.
    
// Open the vcf files.
    v1.openVcf(vcfFiles[0]);
    v2.openVcf(vcfFiles[1]);

// Read in the header information.
    v1.parseHeader();
    v2.parseHeader();
    checkDataSets(v1, v2); // tools.cpp

// Check that the header for the two files contain the same samples.
    if (v1.samples != v2.samples) {
      cerr << "vcf files contain different samples (or sample order)." << endl;
      exit(1);
    } else {
      string taskDescription = "##vcfCTools=intersect " + vcfFiles[0] + ", " + vcfFiles[1];
      writeHeader(output, v1, false, taskDescription); // tools.cpp
    }

// Intersect the two vcf files.
    intersectVcf(v1, v2);

// Check that the input files had the same list of reference sequences.
// If not, it is possible that there were some problems.
    checkReferenceSequences(v1.referenceSequenceVector, v2.referenceSequenceVector); // tools.cpp

// Close the vcf files.
    v1.closeVcf();
    v2.closeVcf();
  }

  return 0;
}