[wip] rust-htslib integration

.gitignore

@@ -2,4 +2,6 @@
 /Cargo.lock
 test.mmi
 */target
-/*.mmi
+/*.mmi
+test_data/*.bam
+test_data/*.mmi

Cargo.toml

@@ -9,25 +9,36 @@ repository = "https://github.com/jguhlin/minimap2-rs"
 categories = ["science"]
 keywords = ["bioinformatics", "fasta", "alignment", "fastq"]
 exclude = [
-	"**/*.fasta",
-	"libsfasta/test_data/",
-	"*.profdata",
-	"*.mmi",
-	"**/*.mmi",
-	"minimap2-sys/"
+  "**/*.fasta",
+  "libsfasta/test_data/",
+  "*.profdata",
+  "*.mmi",
+  "**/*.mmi",
+  "minimap2-sys/",
 ]
 
 
 [dependencies]
 libc = "0.2.134"
-bytelines = "2.4.0"
-simdutf8 = "0.1.4"
-flate2 = { version = "1.0.17", features = ["zlib-ng"], default-features = false }
+
+simdutf8 = { version = "0.1.4", optional = true }
+
+flate2 = { version = "1.0.17", features = [
+  "zlib"
+], default-features = false, optional = true }
 
 # Dep for development
 #minimap2-sys = { path = "./minimap2-sys" }
 minimap2-sys = "0.1.7"
-fffx = "0.1.1"
+
+fffx = {version = "0.1.2", optional = true }
+#fffx = { path = "../fffx", optional=true }
+rust-htslib = { version = "0.40.2", optional = true }
+
+[features]
+default = ['map-file']
+htslib = ['rust-htslib']
+map-file = ['flate2', 'simdutf8', 'fffx']
 
 # [profile.release]
 # opt-level = 3
@@ -39,3 +50,4 @@ debug = true
 
 [profile.dev.package."*"]
 opt-level = 3
+

src/htslib.rs

@@ -0,0 +1,169 @@
+use crate::{Aligner, Mapping, Strand};
+use core::ffi;
+use minimap2_sys::mm_idx_t;
+use rust_htslib::bam::header::HeaderRecord;
+use rust_htslib::bam::record::{Cigar, CigarString};
+use rust_htslib::bam::{Header, Record};
+
+pub fn mapping_to_record(
+    mapping: Option<&Mapping>,
+    seq: &[u8],
+    header: Header,
+    qual: Option<&[u8]>,
+    query_name: Option<&[u8]>,
+) -> Record {
+    let mut rec = Record::new();
+    let qname = query_name.unwrap_or(b"query");
+    // FIXFIX: there's probably a better way of setting a default value
+    // for the quality string
+    let qual = match qual {
+        Some(q) => Vec::from(q),
+        None => {
+            let q = vec![255; seq.len()];
+            q
+        }
+    };
+
+    let cigar: Option<CigarString> = mapping
+        .and_then(|m| m.alignment.clone()) // FIXFIX: we probably don't need a clone here
+        .and_then(|a| a.cigar)
+        .map(|c| cigar_to_cigarstr(&c));
+
+    rec.set(qname, cigar.as_ref(), seq, &qual[..]);
+    match mapping {
+        Some(m) => {
+            println!("Strand {m:?}");
+            if m.strand == Strand::Reverse {
+                println!("here");
+                rec.set_reverse();
+            }
+            // TODO: set secondary/supplementary flags
+            rec.set_pos(m.target_start as i64);
+            rec.set_mapq(m.mapq as u8);
+            rec.set_mpos(-1);
+            // TODO: set tid from sequences listed in header
+            rec.set_mtid(-1);
+            rec.set_insert_size(0);
+        }
+        None => {
+            rec.set_unmapped();
+            rec.set_tid(-1);
+            rec.set_pos(-1);
+            rec.set_mapq(255);
+            rec.set_mpos(-1);
+            rec.set_mtid(-1);
+            rec.set_insert_size(-1);
+        }
+    };
+    // TODO: set AUX flags for cs/md if available
+    rec
+}
+
+fn cigar_to_cigarstr(cigar: &Vec<(u32, u8)>) -> CigarString {
+    let op_vec: Vec<Cigar> = cigar
+        .to_owned()
+        .iter()
+        .map(|(len, op)| match op {
+            0 => Cigar::Match(*len),
+            1 => Cigar::Ins(*len),
+            2 => Cigar::Del(*len),
+            3 => Cigar::RefSkip(*len),
+            4 => Cigar::SoftClip(*len),
+            5 => Cigar::HardClip(*len),
+            6 => Cigar::Pad(*len),
+            7 => Cigar::Equal(*len),
+            8 => Cigar::Diff(*len),
+            _ => panic!("Unexpected cigar operation"),
+        })
+        .collect();
+    CigarString(op_vec)
+}
+
+#[derive(Debug, PartialEq, Eq)]
+pub struct SeqMetaData {
+    pub name: String,
+    pub length: u32,
+    pub is_alt: bool,
+}
+
+#[derive(Debug)]
+pub struct MMIndex {
+    pub inner: mm_idx_t,
+}
+
+impl MMIndex {
+    pub fn n_seq(&self) -> u32 {
+        self.inner.n_seq
+    }
+
+    pub fn seqs(&self) -> Vec<SeqMetaData> {
+        let mut seqs: Vec<SeqMetaData> = Vec::with_capacity(self.n_seq() as usize);
+        for i in 0..self.n_seq() {
+            let _seq = unsafe { *(self.inner.seq).offset(i as isize) };
+            let c_str = unsafe { ffi::CStr::from_ptr(_seq.name) };
+            let rust_str = c_str.to_str().unwrap().to_string();
+            seqs.push(SeqMetaData {
+                name: rust_str,
+                length: _seq.len,
+                is_alt: _seq.is_alt != 0,
+            });
+        }
+        seqs
+    }
+
+    pub fn get_header(&self) -> Header {
+        let mut header = Header::new();
+        for seq in self.seqs() {
+            header.push_record(
+                HeaderRecord::new(b"SQ")
+                    .push_tag(b"SN", &seq.name)
+                    .push_tag(b"LN", &seq.length),
+            );
+        }
+        header
+    }
+}
+
+impl From<&Aligner> for MMIndex {
+    fn from(aligner: &Aligner) -> Self {
+        MMIndex {
+            inner: aligner.idx.unwrap(),
+        }
+    }
+}
+
+#[cfg(test)]
+#[cfg(feature = "htslib")]
+mod tests {
+    use super::*;
+
+    fn get_aligner() -> Aligner {
+        let aligner = Aligner::builder()
+            .with_threads(1)
+            .with_index("test_data/MT-human.fa", Some("test_data/MT-human.mmi"))
+            .unwrap();
+        aligner
+    }
+
+    #[test]
+    fn test_index() {
+        let aligner = get_aligner();
+        let idx = MMIndex::from(&aligner);
+        let seqs = idx.seqs();
+        assert_eq!(
+            seqs,
+            vec![SeqMetaData {
+                name: "MT_human".to_string(),
+                length: 16569u32,
+                is_alt: false
+            }]
+        );
+
+        let header = idx.get_header();
+
+        let records = header.to_hashmap();
+        let observed = records.get("SQ").unwrap().first().unwrap();
+        assert_eq!(observed.get("SN").unwrap(), "MT_human");
+        assert_eq!(observed.get("LN").unwrap(), "16569");
+    }
+}

src/lib.rs

@@ -1,20 +1,27 @@
 use std::cell::RefCell;
-use std::io::BufRead;
+
 use std::io::Read;
 use std::mem::MaybeUninit;
 use std::num::NonZeroI32;
 use std::path::Path;
 
-use flate2::read::GzDecoder;
 use minimap2_sys::*;
+
+#[cfg(feature = "map-file")]
+use flate2::read::GzDecoder;
+#[cfg(feature = "map-file")]
 use simdutf8::basic::from_utf8;
 
+#[cfg(feature = "htslib")]
+pub mod htslib;
+
 /// Alias for mm_mapop_t
 pub type MapOpt = mm_mapopt_t;
 
 /// Alias for mm_idxopt_t
 pub type IdxOpt = mm_idxopt_t;
 
+#[cfg(feature = "map-file")]
 pub use fffx::{Fasta, Fastq, Sequence};
 
 // TODO: Probably a better way to handle this...
@@ -456,7 +463,7 @@ impl Aligner {
 
         let mut idx: MaybeUninit<*mut mm_idx_t> = MaybeUninit::uninit();
 
-        let mut idx_reader = unsafe { idx_reader.assume_init() };
+        let idx_reader = unsafe { idx_reader.assume_init() };
 
         unsafe {
             // Just a test read? Just following: https://github.com/lh3/minimap2/blob/master/python/mappy.pyx#L147
@@ -480,8 +487,8 @@ impl Aligner {
 
     /// Map a single sequence to an index
     /// not implemented yet!
-    pub fn with_seq(mut self, seq: &[u8]) -> Result<Self, &'static str> {
-        let seq = match std::ffi::CString::new(seq) {
+    pub fn with_seq(self, seq: &[u8]) -> Result<Self, &'static str> {
+        let _seq = match std::ffi::CString::new(seq) {
             Ok(seq) => seq,
             Err(_) => return Err("Invalid sequence"),
         };
@@ -625,7 +632,7 @@ impl Aligner {
                             let mut m_cs_string: libc::c_int = 0i32;
 
                             let cs_str = if cs {
-                                let cs_len = mm_gen_cs(
+                                let _cs_len = mm_gen_cs(
                                     km,
                                     &mut cs_string,
                                     &mut m_cs_string,
@@ -644,7 +651,7 @@ impl Aligner {
                             };
 
                             let md_str = if md {
-                                let md_len = mm_gen_MD(
+                                let _md_len = mm_gen_MD(
                                     km,
                                     &mut cs_string,
                                     &mut m_cs_string,
@@ -675,7 +682,6 @@ impl Aligner {
                     } else {
                         None
                     };
-
                     mappings.push(Mapping {
                         target_name: Some(
                             std::ffi::CStr::from_ptr(contig)
@@ -716,6 +722,7 @@ impl Aligner {
     ///
     /// TODO: Remove cs and md and make them options on the struct
     ///
+    #[cfg(feature = "map-file")]
     pub fn map_file(&self, file: &str, cs: bool, md: bool) -> Result<Vec<Mapping>, &'static str> {
         // Make sure index is set
         if self.idx.is_none() {
@@ -762,7 +769,7 @@ impl Aligner {
         }
 
         // If gzipped, open it with a reader...
-        let mut reader: Box<dyn Read> = if compression_type == CompressionType::GZIP {
+        let reader: Box<dyn Read> = if compression_type == CompressionType::GZIP {
             Box::new(GzDecoder::new(std::fs::File::open(file).unwrap()))
         } else {
             Box::new(std::fs::File::open(file).unwrap())
@@ -828,6 +835,7 @@ pub enum FileFormat {
 }
 
 #[allow(dead_code)]
+#[cfg(feature = "map-file")]
 pub fn detect_file_format(buffer: &[u8]) -> Result<FileFormat, &'static str> {
     let buffer = from_utf8(&buffer).expect("Unable to parse file as UTF-8");
     if buffer.starts_with(">") {
@@ -909,7 +917,7 @@ mod tests {
 
     #[test]
     fn test_builder() {
-        let mut aligner = Aligner::builder().preset(Preset::MapOnt);
+        let _aligner = Aligner::builder().preset(Preset::MapOnt);
     }
 
     #[test]
@@ -930,7 +938,7 @@ mod tests {
         let mappings = aligner.map("ACGGTAGAGAGGAAGAAGAAGGAATAGCGGACTTGTGTATTTTATCGTCATTCGTGGTTATCATATAGTTTATTGATTTGAAGACTACGTAAGTAATTTGAGGACTGATTAAAATTTTCTTTTTTAGCTTAGAGTCAATTAAAGAGGGCAAAATTTTCTCAAAAGACCATGGTGCATATGACGATAGCTTTAGTAGTATGGATTGGGCTCTTCTTTCATGGATGTTATTCAGAAGGAGTGATATATCGAGGTGTTTGAAACACCAGCGACACCAGAAGGCTGTGGATGTTAAATCGTAGAACCTATAGACGAGTTCTAAAATATACTTTGGGGTTTTCAGCGATGCAAAA".as_bytes(), false, false, None, None).unwrap();
         println!("{:#?}", mappings);
 
-        let mut aligner = aligner.with_cigar();
+        let aligner = aligner.with_cigar();
 
         aligner
             .map(
@@ -1045,4 +1053,5 @@ b"GTTTATGTAGCTTATTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACACGCCCCATAAACAAATAGGT
             assert_eq!(align.cs.is_some(), *cs);
         }
     }
+
 }

test_data/cDNA_reads.fa

@@ -0,0 +1,4 @@
+>cdna.fwd	el:Z:chr1:0:1720:1
+ATGCCTAGAAGTGTGTGATCGCATTGCTGCCAAGTATTCGATGCATCTGTTACCCAGAGGTGCTCCTCACTACAGCCAGGTCATGGACTTCTTCTCAGGAAAGTCTGGTCGCGGTCGTGGTTGTCCGAGAAACGCATCACCCACAGATAAAATCAGTTATTACAGTTGGACCTTCATGTCAAACCAGAGACCCGTATTTCAAATCGCACATACTGCGTCGTGCAATGCCGGGCGCTAACGGCTCAATATCACGCTGCGTCACTATGGCTACCCCAAAGCGGGGGGGGCATCGACGGGCTGCATAGTGGCAGTATCGAGACGACGACCGTTAAAGAATTTCGGAGGCTGGGTGCTGTACTGTAATCGCCTATTGCAGCACTCAGATTTTAACCTGGTGCCG
+>cdna.rev	el:Z:chr1:360:1720:0
+CGGCACCAGGTTAAAATCTGAGTGCTGCAATAGGCGATTACAGTACAGCACCCAGCCTCCGAAATTCTTTAACGGTCGTCGTCTCGATACTGCCACTATGGTTGGTAGGTCCCAATCTTCGAGCACAGCGGCTCTAACGTGGCCGAAAGCAGTGCCACGATCCGCATGATCTTATGAGCCATACCACATGTAGTTTAGAGCAGCCCGTCGATGCCCCCCCCGCTTTGGGGTAGCCATAGTGACGCAGCGTGATATTGAGCCGTTAGCGCCCGGCATTGCACGACGCAGTATGTGCGATTTTACTCCTGAAGCAGCCGTGAGTACTGTACTCTAATGCACACGTACTTGTTTATGCGAGAAGGTTACCGTCTAATTTTAATGGCCTAGGCGATGGCTATAA

test_data/cDNA_vs_genome.sam

@@ -0,0 +1,6 @@
+@HD	VN:1.6	SO:queryname
+@SQ	SN:chr1	LN:1720
+@PG	ID:minimap2	PN:minimap2	VN:2.24-r1122	CL:minimap2 -ay -x splice --MD --cs test_data/genome.fa test_data/cDNA_reads.fa
+@PG	ID:samtools	PN:samtools	PP:minimap2	VN:1.16.1	CL:samtools sort -n -o test_data/cDNA_vs_genome.sam
+cdna.fwd	0	chr1	1	60	100M620N100M80N100M620N100M	*	0	0	ATGCCTAGAAGTGTGTGATCGCATTGCTGCCAAGTATTCGATGCATCTGTTACCCAGAGGTGCTCCTCACTACAGCCAGGTCATGGACTTCTTCTCAGGAAAGTCTGGTCGCGGTCGTGGTTGTCCGAGAAACGCATCACCCACAGATAAAATCAGTTATTACAGTTGGACCTTCATGTCAAACCAGAGACCCGTATTTCAAATCGCACATACTGCGTCGTGCAATGCCGGGCGCTAACGGCTCAATATCACGCTGCGTCACTATGGCTACCCCAAAGCGGGGGGGGCATCGACGGGCTGCATAGTGGCAGTATCGAGACGACGACCGTTAAAGAATTTCGGAGGCTGGGTGCTGTACTGTAATCGCCTATTGCAGCACTCAGATTTTAACCTGGTGCCG	*	NM:i:0	ms:i:400	AS:i:304	nn:i:0	ts:A:+	tp:A:P	cm:i:124	s1:i:372	s2:i:87	de:f:0	MD:Z:400	rl:i:0	el:Z:chr1:0:1720:1
+cdna.rev	16	chr1	361	60	100M440N100M260N100M260N100M	*	0	0	TTATAGCCATCGCCTAGGCCATTAAAATTAGACGGTAACCTTCTCGCATAAACAAGTACGTGTGCATTAGAGTACAGTACTCACGGCTGCTTCAGGAGTAAAATCGCACATACTGCGTCGTGCAATGCCGGGCGCTAACGGCTCAATATCACGCTGCGTCACTATGGCTACCCCAAAGCGGGGGGGGCATCGACGGGCTGCTCTAAACTACATGTGGTATGGCTCATAAGATCATGCGGATCGTGGCACTGCTTTCGGCCACGTTAGAGCCGCTGTGCTCGAAGATTGGGACCTACCAACCATAGTGGCAGTATCGAGACGACGACCGTTAAAGAATTTCGGAGGCTGGGTGCTGTACTGTAATCGCCTATTGCAGCACTCAGATTTTAACCTGGTGCCG	*	NM:i:0	ms:i:400	AS:i:304	nn:i:0	ts:A:-	tp:A:P	cm:i:120	s1:i:369	s2:i:0	de:f:0	MD:Z:400	rl:i:0	el:Z:chr1:360:1720:0

test_data/gDNA_reads.fa

@@ -0,0 +1,12 @@
+>perfect_read.fwd	el:Z:chr1:180:280:1	em:i:0
+TACGCCACACGGGCTACACTCTCGCCTTCTCGTCGCAACTACGAGCTGGACTATCGGCCGAGAGGATCTAACACGAGAAGTACTTGCCGGCAATCCCTAA
+>perfect_read.rev	el:Z:chr1:180:280:0	em:i:0
+TTAGGGATTGCCGGCAAGTACTTCTCGTGTTAGATCCTCTCGGCCGATAGTCCAGCTCGTAGTTGCGACGAGAAGGCGAGAGTGTAGCCCGTGTGGCGTA
+>imperfect_read.fwd	el:Z:chr1:180:280:1	em:i:5
+TACGCCACACAGGCTACACTCTCGCCTTCTCGTCGCAACTACGCGCTGGACTATCCGCCGAGAGGATCTAACACGTGAAGTACTTGCCGGCATTCCCTAA
+>unmappable_read	el:Z:chr1:180:280:1	em:i:99
+GTTTGTGAGGATCGCGGCGATGGTTTAGGGACGATAAGTATTATCGGAAGGGGCATATTCGTGTCGGAAGGTGAAGATGCCGAAAAAGCAATCGAAGCCG
+>perfect_inv_duplicate	el:Z:chr1:540:640:1	em:i:0
+GAAATACGGGTCTCTGGTTTGACATAAAGGTCCAACTGTAATAACTGATTTTATCTGTGGGTGATGCGTTTCTCGGACAACCACGACCGCGCCCAGACTT
+>split_read	
+ATGCCTAGAAGTGTGTGATCGCATTGCTGCCAAGTATTCGATGCATCTGTTACCCAGAGGTGCTCCTCACTACAGCCAGGTCATGGACTTCTTCTCAGGACTACCCACCTGTTTCATGATCCCCCCTTTGTGAACAATAAACTTAGTAAACATTTTTACGATTAAATGTTTAACTCCTAC

test_data/gDNA_vs_genome.sam

@@ -0,0 +1,12 @@
+@HD	VN:1.6	SO:queryname
+@SQ	SN:chr1	LN:1720
+@PG	ID:minimap2	PN:minimap2	VN:2.24-r1122	CL:minimap2 -ay --MD --cs test_data/genome.fa test_data/gDNA_reads.fa
+@PG	ID:samtools	PN:samtools	PP:minimap2	VN:1.16.1	CL:samtools sort -n -o test_data/gDNA_vs_genome.sam
+imperfect_read.fwd	0	chr1	181	27	100M	*	0	0	TACGCCACACAGGCTACACTCTCGCCTTCTCGTCGCAACTACGCGCTGGACTATCCGCCGAGAGGATCTAACACGTGAAGTACTTGCCGGCATTCCCTAA	*	NM:i:5	ms:i:170	AS:i:170	nn:i:0	tp:A:P	cm:i:6	s1:i:56	s2:i:0	de:f:0.05	MD:Z:10G32A11G19A16A7	rl:i:0	el:Z:chr1:180:280:1	em:i:5
+perfect_inv_duplicate	0	chr1	541	36	100M	*	0	0	GAAATACGGGTCTCTGGTTTGACATAAAGGTCCAACTGTAATAACTGATTTTATCTGTGGGTGATGCGTTTCTCGGACAACCACGACCGCGCCCAGACTT	*	NM:i:0	ms:i:200	AS:i:200	nn:i:0	tp:A:P	cm:i:15	s1:i:85	s2:i:71	de:f:0	MD:Z:100	rl:i:0	el:Z:chr1:540:640:1	em:i:0
+perfect_inv_duplicate	272	chr1	721	0	100M	*	0	0	*	*	NM:i:2	ms:i:188	AS:i:188	nn:i:0	tp:A:S	cm:i:11	s1:i:71	de:f:0.02	MD:Z:8T65C25	rl:i:0	el:Z:chr1:540:640:1	em:i:0
+perfect_read.fwd	0	chr1	181	60	100M	*	0	0	TACGCCACACGGGCTACACTCTCGCCTTCTCGTCGCAACTACGAGCTGGACTATCGGCCGAGAGGATCTAACACGAGAAGTACTTGCCGGCAATCCCTAA	*	NM:i:0	ms:i:200	AS:i:200	nn:i:0	tp:A:P	cm:i:17	s1:i:84	s2:i:0	de:f:0	MD:Z:100	rl:i:0	el:Z:chr1:180:280:1	em:i:0
+perfect_read.rev	16	chr1	181	60	100M	*	0	0	TACGCCACACGGGCTACACTCTCGCCTTCTCGTCGCAACTACGAGCTGGACTATCGGCCGAGAGGATCTAACACGAGAAGTACTTGCCGGCAATCCCTAA	*	NM:i:0	ms:i:200	AS:i:200	nn:i:0	tp:A:P	cm:i:17	s1:i:84	s2:i:0	de:f:0	MD:Z:100	rl:i:0	el:Z:chr1:180:280:0	em:i:0
+split_read	0	chr1	1	60	100M80S	*	0	0	ATGCCTAGAAGTGTGTGATCGCATTGCTGCCAAGTATTCGATGCATCTGTTACCCAGAGGTGCTCCTCACTACAGCCAGGTCATGGACTTCTTCTCAGGACTACCCACCTGTTTCATGATCCCCCCTTTGTGAACAATAAACTTAGTAAACATTTTTACGATTAAATGTTTAACTCCTAC	*	NM:i:0	ms:i:200	AS:i:200	nn:i:0	tp:A:P	cm:i:13	s1:i:90	s2:i:0	de:f:0	SA:Z:chr1,821,-,80M100S,52,0;	MD:Z:100	rl:i:0
+split_read	2064	chr1	821	52	80M100H	*	0	0	GTAGGAGTTAAACATTTAATCGTAAAAATGTTTACTAAGTTTATTGTTCACAAAGGGGGGATCATGAAACAGGTGGGTAG	*	NM:i:0	ms:i:160	AS:i:160	nn:i:0	tp:A:P	cm:i:11	s1:i:70	s2:i:0	de:f:0	SA:Z:chr1,1,+,100M80S,60,0;	MD:Z:80	rl:i:0
+unmappable_read	4	*	0	0	*	*	0	0	GTTTGTGAGGATCGCGGCGATGGTTTAGGGACGATAAGTATTATCGGAAGGGGCATATTCGTGTCGGAAGGTGAAGATGCCGAAAAAGCAATCGAAGCCG	*	rl:i:0	el:Z:chr1:180:280:1	em:i:99

test_data/genome.fa

@@ -0,0 +1,2 @@
+>chr1
+ATGCCTAGAAGTGTGTGATCGCATTGCTGCCAAGTATTCGATGCATCTGTTACCCAGAGGTGCTCCTCACTACAGCCAGGTCATGGACTTCTTCTCAGGAGTAATTTGCGCTGCGGAAAACGGCTGATGGGGAGTCGACCTACCTTAATATCTCCGAGGTTGCCCTCACAAATGGCGTAGTACGCCACACGGGCTACACTCTCGCCTTCTCGTCGCAACTACGAGCTGGACTATCGGCCGAGAGGATCTAACACGAGAAGTACTTGCCGGCAATCCCTAAGTACCTAGGCTCTCGGTCCACTATGACGCAGGACAGGGTTCAGTTAAAAGGCCTCTTCATGCGGTCTTAAGACCTTATAGTTATAGCCATCGCCTAGGCCATTAAAATTAGACGGTAACCTTCTCGCATAAACAAGTACGTGTGCATTAGAGTACAGTACTCACGGCTGCTTCAGGAGTAGTACATAGATGTCCGTTAGACTCTCACGTGTTTGCCACTCAGTTCCGACATGTTTATACGCAATCCTATTGAACGACTAGGAAATACGGGTCTCTGGTTTGACATAAAGGTCCAACTGTAATAACTGATTTTATCTGTGGGTGATGCGTTTCTCGGACAACCACGACCGCGCCCAGACTTGTACTCCGTTCGATGGTGACAGAAGGGTTATGGAGGGAGTGATTATCGCTTTTTACGGAAGGAAAGAGCCGGGTTTGTAGAAGTCTGGTCGCGGTCGTGGTTGTCCGAGAAACGCATCACCCACAGATAAAATCAGTTATTACAGTTGGACCTTCATGTCAAACCAGAGACCCGTATTTCGTAGGAGTTAAACATTTAATCGTAAAAATGTTTACTAAGTTTATTGTTCACAAAGGGGGGATCATGAAACAGGTGGGTAGAAATCGCACATACTGCGTCGTGCAATGCCGGGCGCTAACGGCTCAATATCACGCTGCGTCACTATGGCTACCCCAAAGCGGGGGGGGCATCGACGGGCTGGTATATTGGGCTACGTCCAAAATTTAACGCCAGCAGTACTGTCCGAAGGAGCGTATCTCCTGAGTCTGCAGATGCAGTAGTTTGATTTGAGCTCCATTACCCTACAATTAGAACACTGGCAACATTTGGGCGTTGAGCGGTCTTCCGTGTCGCTCGATCCGCTGGAACTTGGCAACCACAGTACCCAGGACCAGGGGTTGTCTTGTGGCTGTCTAAGGTCGGTGCTATACTGTGGCAACAGTTCTTTCCTACATGAATAGCTCTAAACTACATGTGGTATGGCTCATAAGATCATGCGGATCGTGGCACTGCTTTCGGCCACGTTAGAGCCGCTGTGCTCGAAGATTGGGACCTACCAACGTAAGACCTGTCCCTTTTGGTCCTCAGCTTCGTGAGGGACTACCAACTACAGCTGTGCTCTTGGCTGTCTCAAGCACTAGTGTCAATTGACAAAGATCCGTCAATAGCTATGCATCCCACAACACTAAAGCCACCCGTCTGCACGATCTGCAGCATCACCTGAAGACAATAATATAAAGGGTAGAGGTCACTGGCATTTTCCATGCTGTTACGTACGTAGCAGCGCTATAGCCAGAGGTGCGTCATAGAATAGGAATTAGCATAGTGGCAGTATCGAGACGACGACCGTTAAAGAATTTCGGAGGCTGGGTGCTGTACTGTAATCGCCTATTGCAGCACTCAGATTTTAACCTGGTGCCG