Repository: bacpop/ska.rust (rust)

Issue:
Add option for repeat mapping
In `ska map`, having an option to automatically mask any repeats in the reference would be useful. I think converting bases to 'N', if mapped (leaving as gapped otherwise)

Pull request:
Mask repeats with Ns during ska map
TODO:
- [x] Basic test (aln)
- [x] Basic test (vcf)
- [x] Test with two chromosomes (aln)
- [x] Test with two chromosomes (vcf)

Closes #32

Patch:
diff --git a/src/cli.rs b/src/cli.rs
index 56cc8d2..aef89cd 100644
--- a/src/cli.rs
+++ b/src/cli.rs
@@ -9,6 +9,8 @@ use super::QualFilter;
 pub const DEFAULT_KMER: usize = 17;
 /// Default single strand (which is equivalent to !rc)
 pub const DEFAULT_STRAND: bool = false;
+/// Default repeat masking behaviour
+pub const DEFAULT_REPEATMASK: bool = false;
 /// Default minimum k-mer count for FASTQ files
 pub const DEFAULT_MINCOUNT: u16 = 10;
 /// Default minimum base quality (PHRED score) for FASTQ files
@@ -189,6 +191,10 @@ pub enum Commands {
         #[arg(short, long, value_enum, default_value_t = FileType::Aln)]
         format: FileType,
 
+        /// Mask any repeats in the alignment with 'N'
+        #[arg(long, default_value_t = DEFAULT_REPEATMASK)]
+        repeat_mask: bool,
+
         /// Number of CPU threads
         #[arg(long, value_parser = valid_cpus, default_value_t = 1)]
         threads: usize,
diff --git a/src/generic_modes.rs b/src/generic_modes.rs
index 9793280..2853a16 100644
--- a/src/generic_modes.rs
+++ b/src/generic_modes.rs
@@ -189,7 +189,7 @@ pub fn weed<IntT: for<'a> UInt<'a>>(
             ska_array.kmer_len(),
             ska_array.rc()
         );
-        let ska_weed = RefSka::new(ska_array.kmer_len(), weed_fasta, ska_array.rc());
+        let ska_weed = RefSka::new(ska_array.kmer_len(), weed_fasta, ska_array.rc(), false);
 
         if !reverse {
             log::info!("Removing weed k-mers");
diff --git a/src/lib.rs b/src/lib.rs
index c48bb1d..179418f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -154,6 +154,7 @@
 //! ```bash
 //! ska map ref.fa seqs.skf -o ref_mapped.aln
 //! ```
+//! Add `--repeat-mask` to mask any repeated split k-mers in the reference with 'N'.
 //!
 //! You can also get a VCF as output, which has rows as variants, and only has the
 //! variable sites (but will include unmapped bases as missing).
@@ -488,6 +489,7 @@ pub fn main() {
             input,
             output,
             format,
+            repeat_mask,
             threads,
         } => {
             check_threads(*threads);
@@ -498,8 +500,12 @@ pub fn main() {
                     ska_array.kmer_len(),
                     ska_array.rc()
                 );
-                let mut ska_ref =
-                    RefSka::<u64>::new(ska_array.kmer_len(), reference, ska_array.rc());
+                let mut ska_ref = RefSka::<u64>::new(
+                    ska_array.kmer_len(),
+                    reference,
+                    ska_array.rc(),
+                    *repeat_mask,
+                );
                 map(&mut ska_array, &mut ska_ref, output, format, *threads);
             } else if let Ok(mut ska_array) = load_array::<u128>(input, *threads) {
                 log::info!(
@@ -507,8 +513,12 @@ pub fn main() {
                     ska_array.kmer_len(),
                     ska_array.rc()
                 );
-                let mut ska_ref =
-                    RefSka::<u128>::new(ska_array.kmer_len(), reference, ska_array.rc());
+                let mut ska_ref = RefSka::<u128>::new(
+                    ska_array.kmer_len(),
+                    reference,
+                    ska_array.rc(),
+                    *repeat_mask,
+                );
                 map(&mut ska_array, &mut ska_ref, output, format, *threads);
             } else {
                 panic!("Could not read input file(s): {input:?}");
diff --git a/src/merge_ska_array.rs b/src/merge_ska_array.rs
index d02e67f..c31feff 100644
--- a/src/merge_ska_array.rs
+++ b/src/merge_ska_array.rs
@@ -63,7 +63,8 @@ use crate::cli::FilterType;
 /// ska_array.delete_samples(&[&"test_1"]);
 ///
 /// // Delete k-mers
-/// let ska_weed = RefSka::new(ska_array.kmer_len(), &"tests/test_files_in/weed.fa", ska_array.rc());
+/// let mask_repeats = false;
+/// let ska_weed = RefSka::new(ska_array.kmer_len(), &"tests/test_files_in/weed.fa", ska_array.rc(), mask_repeats);
 /// let reverse = false;
 /// ska_array.weed(&ska_weed, reverse);
 /// ```
diff --git a/src/ska_ref.rs b/src/ska_ref.rs
index 4d4ac50..b057cbd 100644
--- a/src/ska_ref.rs
+++ b/src/ska_ref.rs
@@ -20,7 +20,8 @@
 //! let ska_dict = load_array::<u64>(&["tests/test_files_in/merge.skf".to_string()], threads).unwrap().to_dict();
 //!
 //! // Index a reference sequence
-//! let mut ref_kmers = RefSka::new(ska_dict.kmer_len(), &"tests/test_files_in/test_ref.fa", ska_dict.rc());
+//! let mask_repeats = false;
+//! let mut ref_kmers = RefSka::new(ska_dict.kmer_len(), &"tests/test_files_in/test_ref.fa", ska_dict.rc(), mask_repeats);
 //!
 //! // Run mapping, output an alignment to stdout
 //! ref_kmers.map(&ska_dict);
@@ -31,6 +32,8 @@
 use std::io::Write;
 use std::str;
 
+use hashbrown::hash_set::Entry::*;
+use hashbrown::HashSet;
 use rayon::prelude::*;
 
 use noodles_vcf::{
@@ -102,6 +105,8 @@ where
     chrom_names: Vec<String>,
     /// Sequence, indexed by chromosome, then position
     seq: Vec<Vec<u8>>,
+    /// Where repeats should be masked with 'N'
+    repeat_coors: Vec<usize>,
 
     /// Mapping information
 
@@ -153,7 +158,7 @@ where
     /// - File doesn't exist or can't be opened.
     /// - File cannot be parsed as FASTA (FASTQ is not supported).
     /// - If there are no valid split k-mers.
-    pub fn new(k: usize, filename: &str, rc: bool) -> Self {
+    pub fn new(k: usize, filename: &str, rc: bool, repeat_mask: bool) -> Self {
         if !(5..=63).contains(&k) || k % 2 == 0 {
             panic!("Invalid k-mer length");
         }
@@ -161,6 +166,8 @@ where
         let mut split_kmer_pos = Vec::new();
         let mut seq = Vec::new();
         let mut chrom_names = Vec::new();
+        let mut singles = HashSet::new();
+        let mut repeats = HashSet::new();
 
         let mut reader =
             parse_fastx_file(filename).unwrap_or_else(|_| panic!("Invalid path/file: {filename}",));
@@ -193,6 +200,9 @@ where
                     chrom,
                     rc,
                 });
+                if repeat_mask {
+                    Self::track_repeats(kmer, &mut singles, &mut repeats);
+                }
                 while let Some((kmer, base, rc)) = kmer_it.get_next_kmer() {
                     pos = kmer_it.get_middle_pos();
                     split_kmer_pos.push(RefKmer {
@@ -202,6 +212,9 @@ where
                         chrom,
                         rc,
                     });
+                    if repeat_mask {
+                        Self::track_repeats(kmer, &mut singles, &mut repeats);
+                    }
                 }
             }
             chrom += 1;
@@ -211,17 +224,68 @@ where
             panic!("{filename} has no valid sequence");
         }
 
+        // Find the repeat ranges, and intersect them
+        let mut repeat_coors = Vec::new();
+        if repeat_mask {
+            let half_split_len = (k - 1) / 2;
+            let mut last_chrom = 0;
+            let mut last_end = 0;
+            let mut chrom_offset = 0;
+            for sk in &split_kmer_pos {
+                if sk.chrom > last_chrom {
+                    chrom_offset += seq[last_chrom].len();
+                    last_chrom = sk.chrom;
+                }
+                if repeats.contains(&sk.kmer) {
+                    let start = sk.pos - half_split_len + chrom_offset;
+                    let end = sk.pos + half_split_len + chrom_offset;
+                    let range = if start > last_end || start == 0 {
+                        std::ops::Range {
+                            start,
+                            end: end + 1,
+                        }
+                    } else {
+                        std::ops::Range {
+                            start: last_end + 1,
+                            end: end + 1,
+                        }
+                    };
+                    for pos in range {
+                        repeat_coors.push(pos);
+                    }
+                    last_chrom = sk.chrom;
+                    last_end = end;
+                }
+            }
+            log::info!("Masking {} unique split k-mer repeats spanning {} bases", repeats.len(), repeat_coors.len());
+        }
+
         Self {
             k,
             seq,
             chrom_names,
             split_kmer_pos,
+            repeat_coors,
             mapped_pos: Vec::new(),
             mapped_variants: Array2::zeros((0, 0)),
             mapped_names: Vec::new(),
         }
     }
 
+    // Keeps track of split k-mers in the ref, any found before are moved
+    // to the repeats set
+    fn track_repeats(kmer: IntT, singles: &mut HashSet<IntT>, repeats: &mut HashSet<IntT>) {
+        if let Vacant(rep_kmer) = repeats.entry(kmer) {
+            match singles.entry(kmer) {
+                Vacant(single_entry) => single_entry.insert(),
+                Occupied(single_entry) => {
+                    single_entry.remove();
+                    rep_kmer.insert();
+                }
+            }
+        }
+    }
+
     /// Map split k-mers in a [`MergeSkaDict`] against the refrence sequence.
     ///
     /// Iterates through the referecne list and searches for split k-mers in
@@ -276,7 +340,8 @@ where
             panic!("No split k-mers mapped to reference");
         }
 
-        let mut seq_writers = vec![AlnWriter::new(&self.seq, self.k); self.mapped_names.len()];
+        let mut seq_writers =
+            vec![AlnWriter::new(&self.seq, self.k, &self.repeat_coors); self.mapped_names.len()];
         rayon::ThreadPoolBuilder::new()
             .num_threads(threads)
             .build_global()
@@ -318,7 +383,7 @@ where
         threads: usize,
     ) -> Result<(), needletail::errors::ParseError> {
         if self.chrom_names.len() > 1 {
-            eprintln!("WARNING: Reference contained multiple contigs, in the output they will be concatenated");
+            log::warn!("Reference contained multiple contigs, in the output they will be concatenated");
         }
 
         let alignments = self.pseudoalignment(threads);
diff --git a/src/ska_ref/aln_writer.rs b/src/ska_ref/aln_writer.rs
index 853425a..fb94efb 100644
--- a/src/ska_ref/aln_writer.rs
+++ b/src/ska_ref/aln_writer.rs
@@ -37,14 +37,16 @@ pub struct AlnWriter<'a> {
     /// Whether the finalise function has been run, filling to the end of the
     /// final contig.
     finalised: bool,
-    /// Buffer for amibguous bases, which are written in at the end.
+    /// Repeats to mask, which happens at finalisation
+    repeat_regions: &'a Vec<usize>,
+    /// Buffer for amibguous bases, which are written at finalisation.
     _ambig_out: Vec<(u8, usize)>,
 }
 
 impl<'a> AlnWriter<'a> {
     /// Create a new [`AlnWriter`] taking the reference sequence mapped against
     /// and the k-mer size used for the mapping
-    pub fn new(ref_seq: &'a Vec<Vec<u8>>, k: usize) -> Self {
+    pub fn new(ref_seq: &'a Vec<Vec<u8>>, k: usize, repeat_regions: &'a Vec<usize>) -> Self {
         let total_size = ref_seq.iter().map(|x| x.len()).sum();
         let half_split_len = (k - 1) / 2;
         Self {
@@ -57,6 +59,7 @@ impl<'a> AlnWriter<'a> {
             seq_out: vec![b'-'; total_size],
             half_split_len,
             finalised: false,
+            repeat_regions,
             _ambig_out: Vec::new(),
         }
     }
@@ -153,6 +156,12 @@ impl<'a> AlnWriter<'a> {
             for (ambig_base, ambig_pos) in &self._ambig_out {
                 self.seq_out[*ambig_pos] = *ambig_base;
             }
+            // Mask repeats
+            for repeat_idx in self.repeat_regions {
+                if self.seq_out[*repeat_idx] != b'-' {
+                    self.seq_out[*repeat_idx] = b'N';
+                }
+            }
             self.finalised = true;
         }
     }


Test patch:
diff --git a/tests/map.rs b/tests/map.rs
index 8cb21fc..ed4ea12 100644
--- a/tests/map.rs
+++ b/tests/map.rs
@@ -224,3 +224,62 @@ fn map_rev_comp() {
         .assert()
         .stdout_eq_path(sandbox.file_string("map_vcf_ss.stdout", TestDir::Correct));
 }
+
+// Tests the --repeat-mask option
+#[test]
+fn repeat_mask() {
+    let sandbox = TestSetup::setup();
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("test_ref.fa", TestDir::Input))
+        .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
+        .arg("-v")
+        .arg("--repeat-mask")
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_aln_k9.masked.stdout", TestDir::Correct));
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("test_ref.fa", TestDir::Input))
+        .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
+        .arg("-v")
+        .arg("--repeat-mask")
+        .args(&["--format", "vcf"])
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_vcf_k9.masked.stdout", TestDir::Correct));
+
+    // Two identical chromosomes. All bases masked
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("test_ref_two_chrom.fa", TestDir::Input))
+        .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
+        .arg("-v")
+        .arg("--repeat-mask")
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_all_repeats.masked.stdout", TestDir::Correct));
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("test_ref_two_chrom_repeats.fa", TestDir::Input))
+        .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
+        .arg("-v")
+        .arg("--repeat-mask")
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_aln_two_chrom.masked.stdout", TestDir::Correct));
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("test_ref_two_chrom_repeats.fa", TestDir::Input))
+        .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
+        .arg("-v")
+        .arg("--repeat-mask")
+        .args(&["--format", "vcf"])
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_vcf_two_chrom.masked.stdout", TestDir::Correct));
+}
diff --git a/tests/test_files_in/test_ref_two_chrom_repeats.fa b/tests/test_files_in/test_ref_two_chrom_repeats.fa
new file mode 100644
index 0000000..1dddbb2
--- /dev/null
+++ b/tests/test_files_in/test_ref_two_chrom_repeats.fa
@@ -0,0 +1,4 @@
+>fake_ref_chrom1
+CTTTTAAGAGTGTTTTGATGGTTTACGTACGTACGT
+>fake_ref_chrom2
+GATTTGAAGGCGTTTATGATATTTAGCTGTACGAGAGTCTTTTAAGAG
diff --git a/tests/test_results_correct/map_all_repeats.masked.stdout b/tests/test_results_correct/map_all_repeats.masked.stdout
new file mode 100644
index 0000000..f55ad4f
--- /dev/null
+++ b/tests/test_results_correct/map_all_repeats.masked.stdout
@@ -0,0 +1,4 @@
+>test_1
+------------NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN------------NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+>test_2
+------------NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN------------NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
diff --git a/tests/test_results_correct/map_aln_k9.masked.stdout b/tests/test_results_correct/map_aln_k9.masked.stdout
new file mode 100644
index 0000000..1d37049
--- /dev/null
+++ b/tests/test_results_correct/map_aln_k9.masked.stdout
@@ -0,0 +1,4 @@
+>test_1
+------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGANNNNNNNNNAANNNNNNNNNGATGGTTN------------
+>test_2
+------------GATTTGAAGGCGTTTATGATATATAGCTGTACGANNNNNNNNNAANNNNNNNNNGATGGTTN------------
diff --git a/tests/test_results_correct/map_aln_two_chrom.masked.stdout b/tests/test_results_correct/map_aln_two_chrom.masked.stdout
new file mode 100644
index 0000000..80ee676
--- /dev/null
+++ b/tests/test_results_correct/map_aln_two_chrom.masked.stdout
@@ -0,0 +1,4 @@
+>test_1
+NNNNNNNNNNNNNNNNGATGGTTN------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGANNNNNNNNNNNNNN
+>test_2
+---NNNNNNNNNNNNNGATGGTTN------------GATTTGAAGGCGTTTATGATATATAGCTGTACGANNNNNNNNNNN---
diff --git a/tests/test_results_correct/map_vcf_k9.masked.stdout b/tests/test_results_correct/map_vcf_k9.masked.stdout
new file mode 100644
index 0000000..b3090ab
--- /dev/null
+++ b/tests/test_results_correct/map_vcf_k9.masked.stdout
@@ -0,0 +1,47 @@
+##fileformat=VCFv4.3
+##contig=<ID=fake_ref>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	test_1	test_2
+fake_ref	1	.	A	.	.	.	.	GT	.	.
+fake_ref	2	.	C	.	.	.	.	GT	.	.
+fake_ref	3	.	G	.	.	.	.	GT	.	.
+fake_ref	4	.	T	.	.	.	.	GT	.	.
+fake_ref	5	.	A	.	.	.	.	GT	.	.
+fake_ref	6	.	C	.	.	.	.	GT	.	.
+fake_ref	7	.	G	.	.	.	.	GT	.	.
+fake_ref	8	.	T	.	.	.	.	GT	.	.
+fake_ref	9	.	A	.	.	.	.	GT	.	.
+fake_ref	10	.	C	.	.	.	.	GT	.	.
+fake_ref	11	.	G	.	.	.	.	GT	.	.
+fake_ref	12	.	T	.	.	.	.	GT	.	.
+fake_ref	35	.	T	A	.	.	.	GT	0	1
+fake_ref	47	.	G	N	.	.	.	GT	1	1
+fake_ref	48	.	A	N	.	.	.	GT	1	1
+fake_ref	49	.	G	N	.	.	.	GT	1	1
+fake_ref	50	.	T	N	.	.	.	GT	1	1
+fake_ref	51	.	C	N	.	.	.	GT	1	1
+fake_ref	52	.	T	N	.	.	.	GT	1	1
+fake_ref	53	.	T	N	.	.	.	GT	1	1
+fake_ref	54	.	T	N	.	.	.	GT	1	1
+fake_ref	55	.	T	N	.	.	.	GT	1	1
+fake_ref	58	.	G	N	.	.	.	GT	1	1
+fake_ref	59	.	A	N	.	.	.	GT	1	1
+fake_ref	60	.	G	N	.	.	.	GT	1	1
+fake_ref	61	.	T	N	.	.	.	GT	1	1
+fake_ref	62	.	G	N	.	.	.	GT	1	1
+fake_ref	63	.	T	N	.	.	.	GT	1	1
+fake_ref	64	.	T	N	.	.	.	GT	1	1
+fake_ref	65	.	T	N	.	.	.	GT	1	1
+fake_ref	66	.	T	N	.	.	.	GT	1	1
+fake_ref	74	.	T	N	.	.	.	GT	1	1
+fake_ref	75	.	A	.	.	.	.	GT	.	.
+fake_ref	76	.	C	.	.	.	.	GT	.	.
+fake_ref	77	.	G	.	.	.	.	GT	.	.
+fake_ref	78	.	T	.	.	.	.	GT	.	.
+fake_ref	79	.	A	.	.	.	.	GT	.	.
+fake_ref	80	.	C	.	.	.	.	GT	.	.
+fake_ref	81	.	G	.	.	.	.	GT	.	.
+fake_ref	82	.	T	.	.	.	.	GT	.	.
+fake_ref	83	.	A	.	.	.	.	GT	.	.
+fake_ref	84	.	C	.	.	.	.	GT	.	.
+fake_ref	85	.	G	.	.	.	.	GT	.	.
+fake_ref	86	.	T	.	.	.	.	GT	.	.
diff --git a/tests/test_results_correct/map_vcf_two_chrom.masked.stdout b/tests/test_results_correct/map_vcf_two_chrom.masked.stdout
new file mode 100644
index 0000000..2240ebb
--- /dev/null
+++ b/tests/test_results_correct/map_vcf_two_chrom.masked.stdout
@@ -0,0 +1,48 @@
+##fileformat=VCFv4.3
+##contig=<ID=fake_ref_chrom1>
+##contig=<ID=fake_ref_chrom2>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	test_1	test_2
+fake_ref_chrom1	1	.	C	N	.	.	.	GT	1	.
+fake_ref_chrom1	2	.	T	N	.	.	.	GT	1	.
+fake_ref_chrom1	3	.	T	N	.	.	.	GT	1	.
+fake_ref_chrom1	4	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	5	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	6	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom1	7	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom1	8	.	G	N	.	.	.	GT	1	1
+fake_ref_chrom1	9	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom1	10	.	G	N	.	.	.	GT	1	1
+fake_ref_chrom1	11	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	12	.	G	N	.	.	.	GT	1	1
+fake_ref_chrom1	13	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	14	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	15	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	16	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	24	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom1	25	.	A	.	.	.	.	GT	.	.
+fake_ref_chrom1	26	.	C	.	.	.	.	GT	.	.
+fake_ref_chrom1	27	.	G	.	.	.	.	GT	.	.
+fake_ref_chrom1	28	.	T	.	.	.	.	GT	.	.
+fake_ref_chrom1	29	.	A	.	.	.	.	GT	.	.
+fake_ref_chrom1	30	.	C	.	.	.	.	GT	.	.
+fake_ref_chrom1	31	.	G	.	.	.	.	GT	.	.
+fake_ref_chrom1	32	.	T	.	.	.	.	GT	.	.
+fake_ref_chrom1	33	.	A	.	.	.	.	GT	.	.
+fake_ref_chrom1	34	.	C	.	.	.	.	GT	.	.
+fake_ref_chrom1	35	.	G	.	.	.	.	GT	.	.
+fake_ref_chrom1	36	.	T	.	.	.	.	GT	.	.
+fake_ref_chrom2	23	.	T	A	.	.	.	GT	0	1
+fake_ref_chrom2	35	.	G	N	.	.	.	GT	1	1
+fake_ref_chrom2	36	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom2	37	.	G	N	.	.	.	GT	1	1
+fake_ref_chrom2	38	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom2	39	.	C	N	.	.	.	GT	1	1
+fake_ref_chrom2	40	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom2	41	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom2	42	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom2	43	.	T	N	.	.	.	GT	1	1
+fake_ref_chrom2	44	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom2	45	.	A	N	.	.	.	GT	1	1
+fake_ref_chrom2	46	.	G	N	.	.	.	GT	1	.
+fake_ref_chrom2	47	.	A	N	.	.	.	GT	1	.
+fake_ref_chrom2	48	.	G	N	.	.	.	GT	1	.