Merge pull request #1791 from jltsiren/master

Better haplotype generation

Author: jltsiren

Diff for: deps/gbwt

@@ -45,6 +45,8 @@ void help_index(char** argv) {
          << "    -H, --write-haps FILE  store the threads as sequences in FILE" << endl
          << "    -F, --thread-db FILE   write thread database to FILE" << endl
          << "    -P, --force-phasing    replace unphased genotypes with randomly phased ones" << endl
+         << "    -o, --discard-overlaps skip overlapping alternate alleles if the overlap cannot be resolved" << endl
+         << "    -O, --check-overlaps   print information on overlapping variants to stderr" << endl
          << "    -B, --batch-size N     number of samples per batch (default 200)" << endl
          << "    -R, --range X..Y       process samples X to Y (inclusive)" << endl
          << "    -r, --rename V=P       rename contig V in the VCFs to path P in the graph (may repeat)" << endl
@@ -137,12 +139,13 @@ int main_index(int argc, char** argv) {
     // GBWT
     bool index_haplotypes = false, index_paths = false, index_gam = false;
     vector<string> gam_file_names;
-    bool force_phasing = false;
+    bool force_phasing = false, discard_overlaps = false, check_overlaps = false;
     size_t samples_in_batch = 200; // Samples per batch.
     std::pair<size_t, size_t> sample_range(0, ~(size_t)0); // The semiopen range of samples to process.
     map<string, string> path_to_vcf; // Path name conversion from --rename.
     map<string, pair<size_t, size_t>> regions; // Region restrictions for contigs, in VCF name space, as 0-based exclusive-end ranges.
     unordered_set<string> excluded_samples; // Excluded sample names from --exclude.
+    std::set<std::pair<gbwt::size_type, gbwt::size_type>> overlaps; // Unresolved overlaps in the haplotypes.
 
     // GCSA
     gcsa::size_type kmer_size = gcsa::Key::MAX_LENGTH;
@@ -180,6 +183,8 @@ int main_index(int argc, char** argv) {
             {"gbwt-name", required_argument, 0, 'G'},
             {"write-haps", required_argument, 0, 'H'},
             {"force-phasing", no_argument, 0, 'P'},
+            {"discard-overlaps", no_argument, 0, 'o'},
+            {"check-overlaps", no_argument, 0, 'O'},
             {"batch-size", required_argument, 0, 'B'},
             {"range", required_argument, 0, 'R'},
             {"rename", required_argument, 0, 'r'},
@@ -207,7 +212,7 @@ int main_index(int argc, char** argv) {
         };
 
         int option_index = 0;
-        c = getopt_long (argc, argv, "b:t:px:F:v:TG:H:PB:R:r:I:E:g:i:f:k:X:Z:Vd:maANDP:CM:h",
+        c = getopt_long (argc, argv, "b:t:px:F:v:TM:G:H:PoOB:R:r:I:E:g:i:f:k:X:Z:Vd:maANDCh",
                 long_options, &option_index);
 
         // Detect the end of the options.
@@ -262,6 +267,12 @@ int main_index(int argc, char** argv) {
         case 'P':
             force_phasing = true;
             break;
+        case 'o':
+            discard_overlaps = true;
+            break;
+        case 'O':
+            check_overlaps = true;
+            break;
         case 'B':
             samples_in_batch = std::max(parse<size_t>(optarg), 1ul);
             break;
@@ -616,6 +627,9 @@ int main_index(int argc, char** argv) {
 
                     // Determine the reference nodes for the current variant and create a variant site.
                     // If the variant is not an insertion, there should be a path for the ref allele.
+                    // Otherwise the reference position can be determined from the predecessors of the
+                    // alternate alleles.
+                    // TODO: What if the reference visits the same node several times?
                     var.position--; // Use a 0-based position to get the correct var_name.
                     std::string var_name = make_variant_id(var);
                     std::string ref_path_name = "_alt_" + var_name + "_0";
@@ -634,19 +648,24 @@ int main_index(int argc, char** argv) {
                         bool found = false;
                         for (size_t alt_index = 1; alt_index < var.alleles.size(); alt_index++) {
                             std::string alt_path_name = "_alt_" + var_name + "_" + to_string(alt_index);
+                            size_t candidate_pos = 0;
+                            bool candidate_found = false;
                             auto alt_path_iter = alt_paths.find(alt_path_name);
                             if (alt_path_iter != alt_paths.end()) {
                                 gbwt::vector_type pred_nodes = predecessors(*xg_index, alt_path_iter->second);
                                 for (auto node : pred_nodes) {
                                     size_t pred_pos = variants.firstOccurrence(node);
                                     if (pred_pos != variants.invalid_position()) {
-                                        ref_pos = pred_pos + 1;
+                                        candidate_pos = std::max(candidate_pos, pred_pos + 1);
+                                        candidate_found = true;
                                         found = true;
-                                        break;
                                     }
                                 }
-                                if (found) {
-                                    break;
+                                // For each alternate allele, find the rightmost reference node among
+                                // its predecessors. If multiple alleles have candidates for the
+                                // reference position, choose the leftmost one.
+                                if (candidate_found) {
+                                    ref_pos = std::min(ref_pos, candidate_pos);
                                 }
                             }
                         }
@@ -695,6 +714,9 @@ int main_index(int argc, char** argv) {
                     }
                     cerr << "- Phasing information: " << gbwt::inMegabytes(phasing_bytes) << " MB" << endl;
                 }
+                if (check_overlaps) {
+                    gbwt::checkOverlaps(variants, cerr, true);
+                }
 
                 // Save memory:
                 // - Delete the alt paths if we no longer need them.
@@ -724,12 +746,25 @@ int main_index(int argc, char** argv) {
                                 << "_" << haplotype.phase
                                 << "_" << haplotype.count;
                             store_thread(haplotype.path, sn.str());
+                        },
+                        [&](gbwt::size_type site, gbwt::size_type allele) -> bool {
+                            if (check_overlaps) {
+                                overlaps.insert(std::make_pair(site, allele));
+                            }
+                            return discard_overlaps;
                         });
                     if (show_progress) {
                         cerr << "- Processed samples " << phasings[batch].offset() << " to " << (phasings[batch].offset() + phasings[batch].size() - 1) << endl;
                     }
                 }
-            } // End of contig.
+                if (check_overlaps && !overlaps.empty()) {
+                    cerr << overlaps.size() << " unresolved overlaps:" << endl;
+                    for (auto overlap : overlaps) {
+                        cerr << "- site " << overlap.first << ", allele " << overlap.second << endl;
+                    }
+                    overlaps.clear();
+                }
+            } // End of contigs.
         } // End of haplotypes.
 
         // Store the thread database. Write it to disk if a filename is given,

Diff for: src/subcommand/index_main.cpp

@@ -113,6 +113,8 @@ void VGset::to_xg(xg::XG& index, bool store_threads, const regex& paths_to_take,
 
                 // Sort out all the mappings from the paths we pulled out
                 for(Path& path : paths_taken) {
+                    mappings[path.name()] = map<int64_t, Mapping>(); // We want to include empty paths as well.
+
                     for(size_t i = 0; i < path.mapping_size(); i++) {
                         // For each mapping, file it under its rank if a rank is
                         // specified, or at the last rank otherwise.

Diff for: src/vg_set.cpp

@@ -189,7 +189,7 @@ rm -f xy.vg xy.xg
 
 vg construct -r small/x.fa -v small/x.vcf.gz -a >x.vg
 vg index -x x.xg -v small/x.vcf.gz -H haps.bin x.vg
-is $(du -b haps.bin | cut -f 1) 345 "threads may be exported to binary for use in GBWT construction"
+is $(du -b haps.bin | cut -f 1) 329 "threads may be exported to binary for use in GBWT construction"
 
 rm -f x.vg x.xg part.vg x.gcsa haps.bin x.gbwt
 

Diff for: test/t/06_vg_index.t

@@ -26,11 +26,11 @@ is $(vg stats -N y.vg) 44 "pruning with path restoring leaves the correct number
 is $(vg stats -E y.vg) 48 "pruning with path restoring leaves the correct number of edges"
 rm -f y.vg
 
-# Unfold threads: 1 component, 57 nodes, 68 edges
+# Unfold threads: 1 component, 60 nodes, 72 edges
 vg prune -u -g x.gbwt -e 1 x.vg > y.vg
 is $(vg stats -s y.vg | wc -l) 1 "pruning with path unfolding produces the correct number of components"
-is $(vg stats -N y.vg) 57 "pruning with path unfolding produces the correct number of nodes"
-is $(vg stats -E y.vg) 68 "pruning with path unfolding produces the correct number of edges"
+is $(vg stats -N y.vg) 60 "pruning with path unfolding produces the correct number of nodes"
+is $(vg stats -E y.vg) 72 "pruning with path unfolding produces the correct number of edges"
 rm -f y.vg
 
 rm -f x.vg x.gbwt
@@ -46,8 +46,8 @@ vg index -G y.gbwt -v small/xy2.vcf.gz y.vg
 # Check the range in the node mapping at each point
 is $(head -c 16 xy.mapping | hexdump -v -e '4/4 "%u_"') 99_0_99_0_ "empty mapping starts from the right node id"
 vg prune -u -g x.gbwt -e 1 -a -m xy.mapping x.vg > /dev/null
-is $(head -c 16 xy.mapping | hexdump -v -e '4/4 "%u_"') 99_0_125_0_ "the first unfolded graph adds the correct number of nodes to the mapping"
+is $(head -c 16 xy.mapping | hexdump -v -e '4/4 "%u_"') 99_0_128_0_ "the first unfolded graph adds the correct number of nodes to the mapping"
 vg prune -u -g y.gbwt -e 1 -a -m xy.mapping y.vg > /dev/null
-is $(head -c 16 xy.mapping | hexdump -v -e '4/4 "%u_"') 99_0_150_0_ "the second unfolded graph adds the correct number of nodes to the mapping"
+is $(head -c 16 xy.mapping | hexdump -v -e '4/4 "%u_"') 99_0_156_0_ "the second unfolded graph adds the correct number of nodes to the mapping"
 
 rm -f x.vg y.vg x.gbwt y.gbwt xy.mapping