Skip to content

Commit 72b35e6

Browse files
gkarthikgkarthik
committed
Added test for unpaired reads
1 parent ff768bb commit 72b35e6

File tree

5 files changed

+155
-2
lines changed

5 files changed

+155
-2
lines changed

‎data/test.sim.merged.sorted.bam‎

822 Bytes
Binary file not shown.

‎data/test.sim.merged.sorted.bam.bai‎

96 Bytes
Binary file not shown.

‎data/test_merged.bed‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
gb:NC_009942|Organism:West 68 90 WNV_400_1_LEFT 1 +
2+
gb:NC_009942|Organism:West 68 92 WNV_400_2_LEFT 1 +
3+
gb:NC_009942|Organism:West 270 291 WNV_400_3_LEFT 1 +
4+
gb:NC_009942|Organism:West 269 291 WNV_400_3_LEFT 1 +
5+
gb:NC_009942|Organism:West 361 381 WNV_400_4_LEFT_alt 1 +
6+
gb:NC_009942|Organism:West 152 173 WNV_400_5_LEFT_alt 1 +
7+
gb:NC_009942|Organism:West 170 201 WNV_400_6_LEFT_alt 1 +

‎tests/Makefile.am‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@ LIBS = -lhts -lz -lpthread
22

33
CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror
44

5-
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
6-
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
5+
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
6+
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
77
check_primer_trim_SOURCES = test_primer_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp
88
check_quality_trim_SOURCES = check_quality_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp
99
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
1515
check_primer_bed_SOURCES = test_primer_bed.cpp ../src/primer_bed.cpp
1616
check_getmasked_SOURCES = test_getmasked.cpp ../src/get_masked_amplicons.cpp ../src/primer_bed.cpp
1717
check_removereads_SOURCES = test_removereads.cpp ../src/remove_reads_from_amplicon.cpp ../src/primer_bed.cpp ../src/trim_primer_quality.cpp
18+
check_unpaired_trim_SOURCES = test_unpaired_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp

‎tests/test_unpaired_trim.cpp‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
#include<iostream>
2+
#include "../src/trim_primer_quality.h"
3+
#include "../src/primer_bed.h"
4+
#include "htslib/sam.h"
5+
6+
int main(){
7+
int success = 0;
8+
std::string bam = "../data/test.sim.merged.sorted.bam";
9+
std::vector<primer> primers = populate_from_file("../data/test_merged.bed");
10+
std::string region_;
11+
samFile *in = hts_open(bam.c_str(), "r");
12+
hts_idx_t *idx = sam_index_load(in, bam.c_str());
13+
if (idx == NULL) {
14+
if (sam_index_build2(bam.c_str(), 0, 0) < 0) {
15+
std::cerr << ("Unable to open BAM/SAM index.") << std::endl;
16+
return -1;
17+
} else {
18+
idx = sam_index_load(in, bam.c_str());
19+
if (idx == NULL) {
20+
std::cerr << "Unable to create BAM/SAM index." << std::endl;
21+
return -1;
22+
}
23+
}
24+
}
25+
bam_hdr_t *header = sam_hdr_read(in);
26+
region_.assign(header->target_name[0]);
27+
std::string temp(header->text);
28+
hts_itr_t *iter = NULL;
29+
iter = sam_itr_querys(idx, header, region_.c_str());
30+
bam1_t *aln = bam_init1();
31+
cigar_ t;
32+
uint32_t *cigar;
33+
int primer_ctr = 0;
34+
int forward_primer_indices[] = {1, 1, 5, 5, 6};
35+
int rev_primer_indices[] = {3, -1, -1, -1, 4};
36+
uint8_t cigar_flag[5][16] = {
37+
{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
38+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH},
39+
{BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP},
40+
{BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH},
41+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP}
42+
};
43+
uint32_t cigar_len[5][16] = {
44+
{3, 15, 2, 1, 14, 1, 56, 1, 6, 3, 93, 2, 1, 3, 2, 8},
45+
{1, 19, 1, 56, 1, 6, 3, 99, 2, 65},
46+
{8, 13, 1, 6, 3, 99, 2, 62, 3},
47+
{9, 10, 3, 97, 2, 26},
48+
{18, 71, 2, 89, 18}
49+
};
50+
uint32_t read_start_pos[5] = {
51+
94, 92, 173, 175, 201
52+
};
53+
uint8_t condense_cigar_flag[5][14] = {
54+
{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},
55+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH},
56+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP},
57+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH},
58+
{BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CSOFT_CLIP}
59+
};
60+
uint32_t condense_cigar_len[5][13] = {
61+
{18, 2, 1, 14,1, 56, 1, 6, 3, 93, 2, 1, 13},
62+
{1, 19, 1, 56, 1, 6, 3, 99, 2, 65},
63+
{31, 99, 2, 62, 3},
64+
{22, 97, 2, 26},
65+
{18, 71, 2, 89, 18}
66+
};
67+
unsigned int fwd_overlapping_primer_sizes[] = {2, 1, 1, 1, 1};
68+
unsigned int rev_overlapping_primer_sizes[] = {2, 0, 0, 0, 1};
69+
std::vector<primer> overlapping_primers;
70+
primer cand_primer;
71+
while(sam_itr_next(in, iter, aln) >= 0) {
72+
if((aln->core.flag&BAM_FUNMAP) != 0){
73+
continue;
74+
}
75+
std::cout << bam_get_qname(aln) << std::endl;
76+
std::cout << std::endl << "Forward" << std::endl;
77+
get_overlapping_primers(aln, primers, overlapping_primers, false);
78+
// Forward primer
79+
if(overlapping_primers.size() > 0){
80+
cand_primer = get_max_end(overlapping_primers);
81+
t = primer_trim(aln, cand_primer.get_end() + 1, false);
82+
aln->core.pos += t.start_pos;
83+
replace_cigar(aln, t.nlength, t.cigar);
84+
if(overlapping_primers.size() != fwd_overlapping_primer_sizes[primer_ctr]){
85+
success = -1;
86+
std::cout << "Overlapping primer sizes for " << bam_get_qname(aln) <<". Expected: " << fwd_overlapping_primer_sizes[primer_ctr] << ". Got: " << overlapping_primers.size() << std::endl;
87+
}
88+
if(cand_primer.get_indice() != forward_primer_indices[primer_ctr]){
89+
success = -1;
90+
std::cout << "Primer indice wrong. Expected: " << forward_primer_indices[primer_ctr] << ". Got: " << cand_primer.get_indice() << std::endl;
91+
}
92+
}
93+
// Reverse primer
94+
std::cout << std::endl << "Reverse" << std::endl;
95+
get_overlapping_primers(aln, primers, overlapping_primers, true);
96+
if(overlapping_primers.size() > 0){
97+
cand_primer = get_min_start(overlapping_primers);
98+
t = primer_trim(aln, cand_primer.get_start() - 1, true);
99+
replace_cigar(aln, t.nlength, t.cigar);
100+
if(overlapping_primers.size() != rev_overlapping_primer_sizes[primer_ctr]){
101+
success = -1;
102+
std::cout << "Overlapping primer sizes for " << bam_get_qname(aln) <<". Expected: " << rev_overlapping_primer_sizes[primer_ctr] << ". Got: " << overlapping_primers.size() << std::endl;
103+
}
104+
if(cand_primer.get_indice() != rev_primer_indices[primer_ctr]){
105+
success = -1;
106+
std::cout << "Primer indice wrong. Expected: " << rev_primer_indices[primer_ctr] << ". Got: " << cand_primer.get_indice() << std::endl;
107+
}
108+
}
109+
if(aln->core.pos != read_start_pos[primer_ctr]){
110+
success = -1;
111+
std::cout << "Start pos didn't match" << std::endl;
112+
}
113+
replace_cigar(aln, t.nlength, t.cigar);
114+
cigar = bam_get_cigar(aln);
115+
for (uint i = 0; i < t.nlength; ++i){
116+
if(((cigar[i]) & BAM_CIGAR_MASK) != cigar_flag[primer_ctr][i]){
117+
success = -1;
118+
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;
119+
}
120+
if((((cigar[i]) >> BAM_CIGAR_SHIFT)) != cigar_len[primer_ctr][i]){
121+
success = -1;
122+
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;
123+
}
124+
}
125+
// Condense cigar
126+
std::cout << std::endl << "Condensing cigar ... " << std::endl;
127+
t = condense_cigar(t.cigar, t.nlength);
128+
replace_cigar(aln, t.nlength, t.cigar);
129+
cigar = bam_get_cigar(aln);
130+
for (uint i = 0; i < t.nlength; ++i){
131+
if(((cigar[i]) & BAM_CIGAR_MASK) != condense_cigar_flag[primer_ctr][i]){
132+
success = -1;
133+
std::cout << "Cigar flag didn't match! Expected " << condense_cigar_flag[primer_ctr][i] << " " << "Got " << ((cigar[i]) & BAM_CIGAR_MASK) << std::endl;
134+
}
135+
if((((cigar[i]) >> BAM_CIGAR_SHIFT)) != condense_cigar_len[primer_ctr][i]){
136+
success = -1;
137+
std::cout << "Cigar length didn't match after condense! Expected " << condense_cigar_len[primer_ctr][i] << " " << "Got " << ((cigar[i]) >> BAM_CIGAR_SHIFT) << std::endl;
138+
}
139+
}
140+
primer_ctr++;
141+
}
142+
// Check if primers found at all
143+
success = (primer_ctr > 0) ? success : -1;
144+
return success;
145+
}

0 commit comments

Comments
 (0)