diff --git a/data/test.sim.merged.sorted.bam b/data/test.sim.merged.sorted.bam new file mode 100644 index 00000000..7fb71e0a Binary files /dev/null and b/data/test.sim.merged.sorted.bam differ diff --git a/data/test.sim.merged.sorted.bam.bai b/data/test.sim.merged.sorted.bam.bai new file mode 100644 index 00000000..e4c78c21 Binary files /dev/null and b/data/test.sim.merged.sorted.bam.bai differ diff --git a/data/test_merged.bed b/data/test_merged.bed new file mode 100644 index 00000000..382128cd --- /dev/null +++ b/data/test_merged.bed @@ -0,0 +1,7 @@ +gb:NC_009942|Organism:West 68 90 WNV_400_1_LEFT 1 + +gb:NC_009942|Organism:West 68 92 WNV_400_2_LEFT 1 + +gb:NC_009942|Organism:West 270 291 WNV_400_3_LEFT 1 + +gb:NC_009942|Organism:West 269 291 WNV_400_3_LEFT 1 + +gb:NC_009942|Organism:West 361 381 WNV_400_4_LEFT_alt 1 + +gb:NC_009942|Organism:West 152 173 WNV_400_5_LEFT_alt 1 + +gb:NC_009942|Organism:West 170 201 WNV_400_6_LEFT_alt 1 + diff --git a/tests/Makefile.am b/tests/Makefile.am index 64fc6c48..75959a09 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -2,8 +2,8 @@ LIBS = -lhts -lz -lpthread CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror -TESTS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants -check_PROGRAMS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants +TESTS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim +check_PROGRAMS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_SOURCES = test_primer_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp check_quality_trim_SOURCES = check_quality_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp check_consensus_SOURCES = test_call_consensus_from_plup.cpp ../src/call_consensus_pileup.cpp ../src/allele_functions.cpp @@ -15,3 +15,4 @@ check_common_variants_SOURCES = test_common_variants.cpp ../src/get_common_varia check_primer_bed_SOURCES = test_primer_bed.cpp ../src/primer_bed.cpp check_getmasked_SOURCES = test_getmasked.cpp ../src/get_masked_amplicons.cpp ../src/primer_bed.cpp check_removereads_SOURCES = test_removereads.cpp ../src/remove_reads_from_amplicon.cpp ../src/primer_bed.cpp ../src/trim_primer_quality.cpp +check_unpaired_trim_SOURCES = test_unpaired_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp diff --git a/tests/test_unpaired_trim.cpp b/tests/test_unpaired_trim.cpp new file mode 100644 index 00000000..c113a2df --- /dev/null +++ b/tests/test_unpaired_trim.cpp @@ -0,0 +1,145 @@ +#include +#include "../src/trim_primer_quality.h" +#include "../src/primer_bed.h" +#include "htslib/sam.h" + +int main(){ + int success = 0; + std::string bam = "../data/test.sim.merged.sorted.bam"; + std::vector primers = populate_from_file("../data/test_merged.bed"); + std::string region_; + samFile *in = hts_open(bam.c_str(), "r"); + hts_idx_t *idx = sam_index_load(in, bam.c_str()); + if (idx == NULL) { + if (sam_index_build2(bam.c_str(), 0, 0) < 0) { + std::cerr << ("Unable to open BAM/SAM index.") << std::endl; + return -1; + } else { + idx = sam_index_load(in, bam.c_str()); + if (idx == NULL) { + std::cerr << "Unable to create BAM/SAM index." << std::endl; + return -1; + } + } + } + bam_hdr_t *header = sam_hdr_read(in); + region_.assign(header->target_name[0]); + std::string temp(header->text); + hts_itr_t *iter = NULL; + iter = sam_itr_querys(idx, header, region_.c_str()); + bam1_t *aln = bam_init1(); + cigar_ t; + uint32_t *cigar; + int primer_ctr = 0; + int forward_primer_indices[] = {1, 1, 5, 5, 6}; + int rev_primer_indices[] = {3, -1, -1, -1, 4}; + uint8_t cigar_flag[5][16] = { + {BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CINS, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP}, // 3S15M5D2M1D14M1I56M1P6M3P93M2D4M2I8M + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP}, + {BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP} + }; + uint32_t cigar_len[5][16] = { + {3, 15, 2, 1, 14, 1, 56, 1, 6, 3, 93, 2, 1, 3, 2, 8}, + {1, 19, 1, 56, 1, 6, 3, 99, 2, 65}, + {8, 13, 1, 6, 3, 99, 2, 62, 3}, + {9, 10, 3, 97, 2, 26}, + {18, 71, 2, 89, 18} + }; + uint32_t read_start_pos[5] = { + 94, 92, 173, 175, 201 + }; + uint8_t condense_cigar_flag[5][14] = { + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CINS, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CSOFT_CLIP}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP} + }; + uint32_t condense_cigar_len[5][13] = { + {18, 2, 1, 14,1, 56, 1, 6, 3, 93, 2, 1, 13}, + {1, 19, 1, 56, 1, 6, 3, 99, 2, 65}, + {31, 99, 2, 62, 3}, + {22, 97, 2, 26}, + {18, 71, 2, 89, 18} + }; + unsigned int fwd_overlapping_primer_sizes[] = {2, 1, 1, 1, 1}; + unsigned int rev_overlapping_primer_sizes[] = {2, 0, 0, 0, 1}; + std::vector overlapping_primers; + primer cand_primer; + while(sam_itr_next(in, iter, aln) >= 0) { + if((aln->core.flag&BAM_FUNMAP) != 0){ + continue; + } + std::cout << bam_get_qname(aln) << std::endl; + std::cout << std::endl << "Forward" << std::endl; + get_overlapping_primers(aln, primers, overlapping_primers, false); + // Forward primer + if(overlapping_primers.size() > 0){ + cand_primer = get_max_end(overlapping_primers); + t = primer_trim(aln, cand_primer.get_end() + 1, false); + aln->core.pos += t.start_pos; + replace_cigar(aln, t.nlength, t.cigar); + if(overlapping_primers.size() != fwd_overlapping_primer_sizes[primer_ctr]){ + success = -1; + std::cout << "Overlapping primer sizes for " << bam_get_qname(aln) <<". Expected: " << fwd_overlapping_primer_sizes[primer_ctr] << ". Got: " << overlapping_primers.size() << std::endl; + } + if(cand_primer.get_indice() != forward_primer_indices[primer_ctr]){ + success = -1; + std::cout << "Primer indice wrong. Expected: " << forward_primer_indices[primer_ctr] << ". Got: " << cand_primer.get_indice() << std::endl; + } + } + // Reverse primer + std::cout << std::endl << "Reverse" << std::endl; + get_overlapping_primers(aln, primers, overlapping_primers, true); + if(overlapping_primers.size() > 0){ + cand_primer = get_min_start(overlapping_primers); + t = primer_trim(aln, cand_primer.get_start() - 1, true); + replace_cigar(aln, t.nlength, t.cigar); + if(overlapping_primers.size() != rev_overlapping_primer_sizes[primer_ctr]){ + success = -1; + std::cout << "Overlapping primer sizes for " << bam_get_qname(aln) <<". Expected: " << rev_overlapping_primer_sizes[primer_ctr] << ". Got: " << overlapping_primers.size() << std::endl; + } + if(cand_primer.get_indice() != rev_primer_indices[primer_ctr]){ + success = -1; + std::cout << "Primer indice wrong. Expected: " << rev_primer_indices[primer_ctr] << ". Got: " << cand_primer.get_indice() << std::endl; + } + } + if(aln->core.pos != read_start_pos[primer_ctr]){ + success = -1; + std::cout << "Start pos didn't match" << std::endl; + } + replace_cigar(aln, t.nlength, t.cigar); + cigar = bam_get_cigar(aln); + for (uint i = 0; i < t.nlength; ++i){ + if(((cigar[i]) & BAM_CIGAR_MASK) != cigar_flag[primer_ctr][i]){ + success = -1; + std::cout << "Cigar flag didn't match for " << bam_get_qname(aln) << " ! Expected " << (uint) cigar_flag[primer_ctr][i] << " " << "Got " << ((cigar[i]) & BAM_CIGAR_MASK) << std::endl; + } + if((((cigar[i]) >> BAM_CIGAR_SHIFT)) != cigar_len[primer_ctr][i]){ + success = -1; + std::cout << "Cigar length didn't match for " << bam_get_qname(aln) << " ! Expected " << (uint) cigar_len[primer_ctr][i] << " " << "Got " << ((cigar[i]) >> BAM_CIGAR_SHIFT) << std::endl; + } + } + // Condense cigar + std::cout << std::endl << "Condensing cigar ... " << std::endl; + t = condense_cigar(t.cigar, t.nlength); + replace_cigar(aln, t.nlength, t.cigar); + cigar = bam_get_cigar(aln); + for (uint i = 0; i < t.nlength; ++i){ + if(((cigar[i]) & BAM_CIGAR_MASK) != condense_cigar_flag[primer_ctr][i]){ + success = -1; + std::cout << "Cigar flag didn't match! Expected " << condense_cigar_flag[primer_ctr][i] << " " << "Got " << ((cigar[i]) & BAM_CIGAR_MASK) << std::endl; + } + if((((cigar[i]) >> BAM_CIGAR_SHIFT)) != condense_cigar_len[primer_ctr][i]){ + success = -1; + std::cout << "Cigar length didn't match after condense! Expected " << condense_cigar_len[primer_ctr][i] << " " << "Got " << ((cigar[i]) >> BAM_CIGAR_SHIFT) << std::endl; + } + } + primer_ctr++; + } + // Check if primers found at all + success = (primer_ctr > 0) ? success : -1; + return success; +}