From 97c35c48fcb531a6b22bb44682245250c49a6f69 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 27 Jul 2022 18:10:22 -0700 Subject: [PATCH 01/22] adding rtrim function for trailing char in primer files --- src/primer_bed.cpp | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index 9fa99867..7b2d551f 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -233,22 +233,38 @@ std::vector get_primers(std::vector p, unsigned int pos){ return primers_with_mismatches; } +//function to trim trailing spaces +std::string& rtrim(std::string& str, const std::string& chars = "\t\n\v\f\r ") +{ + str.erase(str.find_last_not_of(chars) + 1); + return str; +} + // Assumes unique primer names in BED file +//returns the index of a primer by name int get_primer_indice(std::vector p, std::string name){ + //iterate through the primers from the bed file for(std::vector::iterator it = p.begin(); it != p.end(); ++it) { - if(it->get_name().compare(name) == 0){ + //check if the two strings are the same + if(it->get_name().compare(rtrim(name)) == 0){ return it - p.begin(); } } return -1; } +//function using the tab seperated primer pair file int populate_pair_indices(std::vector &primers, std::string path){ + /* + * @param path the path to the primer pair file + */ + //load primer pair file std::ifstream fin(path.c_str()); std::string line, cell, p1,p2; std::stringstream line_stream; std::vector::iterator it; int32_t indice; + //iterate the primer pair file line by line while (std::getline(fin, line)){ line_stream << line; std::getline(line_stream, cell, '\t'); @@ -259,7 +275,9 @@ int populate_pair_indices(std::vector &primers, std::string path){ line_stream.clear(); if(!p1.empty() && !p2.empty()){ for(it = primers.begin(); it != primers.end(); ++it) { + //search for primer name in pair file if (it->get_name() == p1) { + //make sure it's pair exists indice = get_primer_indice(primers, p2); if (indice != -1) it->set_pair_indice(indice); @@ -273,6 +291,8 @@ int populate_pair_indices(std::vector &primers, std::string path){ std::cout << "Primer pair for " << p2 << " not found in BED file." << std::endl; } } + } else { + std::cout << "Primer pair is empty." << std::endl; } } return 0; @@ -287,4 +307,4 @@ primer get_min_start(std::vector primers){ primer get_max_end(std::vector primers){ auto minmax_start = std::minmax_element(primers.begin(), primers.end(), [] (primer lhs, primer rhs) {return lhs.get_end() < rhs.get_end();}); return *(minmax_start.second); -} \ No newline at end of file +} From 6b9ed985a862dff5f01445b060f2e42df8809cac Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 12:44:50 -0700 Subject: [PATCH 02/22] adding new variable to primer_trim and adding call to populate primer pairs --- src/trim_primer_quality.cpp | 67 +++++++++++++++++++++---------------- 1 file changed, 39 insertions(+), 28 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index 37d53e6d..ae56c22b 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -160,7 +160,11 @@ void print_cigar(uint32_t *cigar, int nlength){ std::cout << std::endl; } -cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_rev = false){ +cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, int32_t end_pos=-1, bool unpaired_rev = false){ + /* + * @param end_pos: the start of the opposite primer + */ + std::cout << end_pos << "\n"; uint32_t *ncigar = (uint32_t*) malloc(sizeof(uint32_t) * (r->core.n_cigar + 1)), // Maximum edit is one more element with soft mask *cigar = bam_get_cigar(r); uint32_t i = 0, j = 0; @@ -452,6 +456,11 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, if(!pair_info.empty()){ amplicons = populate_amplicons(pair_info, primers); } + + //this call is technically redundant but it doesn't stick when + //called in populate _amplicons + populate_pair_indices(primers, pair_info); + std::cout << "Amplicons detected: " << std::endl; amplicons.inOrder(); if(bam.empty()){ @@ -562,35 +571,37 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, } isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); // if reverse strand - if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired - get_overlapping_primers(aln, primers, overlapping_primers); - if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) - primer_trimmed = true; - if(bam_is_rev(aln)){ // Reverse read - cand_primer = get_min_start(overlapping_primers); // fetch reverse primer (?) - t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false); - } else { // Forward read - cand_primer = get_max_end(overlapping_primers); // fetch forward primer (?) - t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); - aln->core.pos += t.start_pos; + if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired + get_overlapping_primers(aln, primers, overlapping_primers); + if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) + primer_trimmed = true; + if(bam_is_rev(aln)){ // Reverse read + cand_primer = get_min_start(overlapping_primers); // fetch reverse primer (?) + + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false); + } else { // Forward read + cand_primer = get_max_end(overlapping_primers); // fetch forward primer (?) + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); + aln->core.pos += t.start_pos; + } + replace_cigar(aln, t.nlength, t.cigar); + free(t.cigar); + // Add count to primer + cit = std::find(primers.begin(), primers.end(), cand_primer); + if(cit != primers.end()) + cit->add_read_count(1); } + t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming + if(bam_is_rev(aln)) // if reverse strand + aln->core.pos = t.start_pos; + condense_cigar(&t); + // aln->core.pos += t.start_pos; replace_cigar(aln, t.nlength, t.cigar); - free(t.cigar); - // Add count to primer - cit = std::find(primers.begin(), primers.end(), cand_primer); - if(cit != primers.end()) - cit->add_read_count(1); - } - t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming - if(bam_is_rev(aln)) // if reverse strand - aln->core.pos = t.start_pos; - condense_cigar(&t); - // aln->core.pos += t.start_pos; - replace_cigar(aln, t.nlength, t.cigar); - } else { // Unpaired reads: Might be stitched reads - if(abs(aln->core.isize) <= abs(aln->core.l_qseq)){ - failed_frag_size++; - } + } else { // Unpaired reads: Might be stitched reads + if(abs(aln->core.isize) <= abs(aln->core.l_qseq)){ + failed_frag_size++; + } + // Forward primer get_overlapping_primers(aln, primers, overlapping_primers, false); if(overlapping_primers.size() > 0){ From f1915505dca8280eaa07cc0c9b539c62558b04c8 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 14:28:32 -0700 Subject: [PATCH 03/22] adding print code for primers to help debug --- src/primer_bed.cpp | 17 +++++++++++++++++ src/primer_bed.h | 2 ++ 2 files changed, 19 insertions(+) diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index 7b2d551f..9772f788 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -308,3 +308,20 @@ primer get_max_end(std::vector primers){ auto minmax_start = std::minmax_element(primers.begin(), primers.end(), [] (primer lhs, primer rhs) {return lhs.get_end() < rhs.get_end();}); return *(minmax_start.second); } + +void print_all_primer_info(std::vector primers){ + std::vector::iterator it; + for(it = primers.begin(); it != primers.end(); ++it) { + std::cout<< "Primer start: " << it->get_start() << std::endl; + std::cout<< "Primer end: " << it->get_end() << std::endl; + std::cout<< "Indice: " << it->get_indice() << std::endl; + std::cout<< "Pair indice: " << it->get_pair_indice() << std::endl; + } +} + +void print_primer_info(primer primer){ + std::cout<< "Primer start: " << primer.get_start() << std::endl; + std::cout<< "Primer end: " << primer.get_end() << std::endl; + std::cout<< "Indice: " << primer.get_indice() << std::endl; + std::cout<< "Pair indice: " << primer.get_pair_indice() << std::endl; +} diff --git a/src/primer_bed.h b/src/primer_bed.h index c8aa71ca..91e2f876 100644 --- a/src/primer_bed.h +++ b/src/primer_bed.h @@ -51,6 +51,8 @@ std::vector populate_from_file(std::string path); std::vector get_primers(std::vector p, unsigned int pos); int get_primer_indice(std::vector p, std::string name); int populate_pair_indices(std::vector &primers, std::string path); +void print_primer_info(primer primers); +void print_all_primer_info(std::vector primers); primer get_min_start(std::vector primers); primer get_max_end(std::vector primers); From 07dff38e050e6939ec2effd048b2c72c2c5c18a0 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 14:40:33 -0700 Subject: [PATCH 04/22] changing the return types of the primer pair populate function to retain pair indice --- src/primer_bed.cpp | 12 +++++++----- src/primer_bed.h | 2 +- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index 9772f788..70c664ae 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -254,10 +254,12 @@ int get_primer_indice(std::vector p, std::string name){ } //function using the tab seperated primer pair file -int populate_pair_indices(std::vector &primers, std::string path){ +void populate_pair_indices(std::vector &primers, std::string path){ /* - * @param path the path to the primer pair file + * @param primers: the primer vector to add pair info to + * @param path: the path to the primer pair file */ + //load primer pair file std::ifstream fin(path.c_str()); std::string line, cell, p1,p2; @@ -279,10 +281,11 @@ int populate_pair_indices(std::vector &primers, std::string path){ if (it->get_name() == p1) { //make sure it's pair exists indice = get_primer_indice(primers, p2); - if (indice != -1) + if (indice != -1){ it->set_pair_indice(indice); - else + }else{ std::cout << "Primer pair for " << p1 << " not found in BED file." <get_name() == p2){ indice = get_primer_indice(primers, p1); if(indice != -1) @@ -295,7 +298,6 @@ int populate_pair_indices(std::vector &primers, std::string path){ std::cout << "Primer pair is empty." << std::endl; } } - return 0; } primer get_min_start(std::vector primers){ diff --git a/src/primer_bed.h b/src/primer_bed.h index 91e2f876..6918a38c 100644 --- a/src/primer_bed.h +++ b/src/primer_bed.h @@ -50,7 +50,7 @@ std::vector populate_from_file(std::string path, int32_t offset); std::vector populate_from_file(std::string path); std::vector get_primers(std::vector p, unsigned int pos); int get_primer_indice(std::vector p, std::string name); -int populate_pair_indices(std::vector &primers, std::string path); +void populate_pair_indices(std::vector &primers, std::string path); void print_primer_info(primer primers); void print_all_primer_info(std::vector primers); primer get_min_start(std::vector primers); From 93dfdac834b5d143d1761d58d783581b1e5dcd6e Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 14:45:05 -0700 Subject: [PATCH 05/22] adding name to print primer statements --- src/primer_bed.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index 70c664ae..dbab7091 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -314,6 +314,7 @@ primer get_max_end(std::vector primers){ void print_all_primer_info(std::vector primers){ std::vector::iterator it; for(it = primers.begin(); it != primers.end(); ++it) { + std::cout << "Primer name: " << it->get_name() << std::endl; std::cout<< "Primer start: " << it->get_start() << std::endl; std::cout<< "Primer end: " << it->get_end() << std::endl; std::cout<< "Indice: " << it->get_indice() << std::endl; @@ -322,6 +323,7 @@ void print_all_primer_info(std::vector primers){ } void print_primer_info(primer primer){ + std::cout << "Primer name: " << primer.get_name() << std::endl; std::cout<< "Primer start: " << primer.get_start() << std::endl; std::cout<< "Primer end: " << primer.get_end() << std::endl; std::cout<< "Indice: " << primer.get_indice() << std::endl; From fc8a3be6b1ff00ba6171b77134851f8c3901ed50 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 15:56:09 -0700 Subject: [PATCH 06/22] add function to find primer pair --- src/primer_bed.cpp | 16 ++++++++++++++++ src/primer_bed.h | 1 + 2 files changed, 17 insertions(+) diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index dbab7091..f53e541a 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -311,6 +311,22 @@ primer get_max_end(std::vector primers){ return *(minmax_start.second); } +//give a primer index, find its pair and return the primer +primer fetch_primer_pair(int16_t index, std::vector primers){ + /* + * @param index: the pair index for the primer of interest + * @return pair_primer: the primer object that's pairs with the primer of interest + */ + std::vector::iterator it; + primer pair_primer; + for(it = primers.begin(); it != primers.end(); ++it) { + if(it->get_indice() == index){ + return *it; + } + } + return(pair_primer); +} + void print_all_primer_info(std::vector primers){ std::vector::iterator it; for(it = primers.begin(); it != primers.end(); ++it) { diff --git a/src/primer_bed.h b/src/primer_bed.h index 6918a38c..04a5d152 100644 --- a/src/primer_bed.h +++ b/src/primer_bed.h @@ -50,6 +50,7 @@ std::vector populate_from_file(std::string path, int32_t offset); std::vector populate_from_file(std::string path); std::vector get_primers(std::vector p, unsigned int pos); int get_primer_indice(std::vector p, std::string name); +primer fetch_primer_pair(int16_t index, std::vector primers); void populate_pair_indices(std::vector &primers, std::string path); void print_primer_info(primer primers); void print_all_primer_info(std::vector primers); From 3ce6ee715927a6d59f67afae7bb6c853daa255fb Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 29 Jul 2022 16:01:23 -0700 Subject: [PATCH 07/22] passing appropriate information to primer trim function --- src/trim_primer_quality.cpp | 82 +++++++++++++++++++++++++++---------- 1 file changed, 61 insertions(+), 21 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index ae56c22b..ff85ab87 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -160,16 +160,19 @@ void print_cigar(uint32_t *cigar, int nlength){ std::cout << std::endl; } -cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, int32_t end_pos=-1, bool unpaired_rev = false){ +cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_rev = false, int32_t end_pos=-1){ /* + * @param r: alignment object * @param end_pos: the start of the opposite primer */ - std::cout << end_pos << "\n"; uint32_t *ncigar = (uint32_t*) malloc(sizeof(uint32_t) * (r->core.n_cigar + 1)), // Maximum edit is one more element with soft mask - *cigar = bam_get_cigar(r); + *cigar = bam_get_cigar(r); uint32_t i = 0, j = 0; int max_del_len = 0, cig, temp, del_len = 0; bool reverse = false; + + //std::cout << "start pos " << new_pos << " end pos " << end_pos << "\n"; + if((r->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired and isize > read length if (bam_is_rev(r)){ // If -ve strand (?) max_del_len = bam_cigar2qlen(r->core.n_cigar, bam_get_cigar(r)) - get_pos_on_query(cigar, r->core.n_cigar, new_pos, r->core.pos) - 1; @@ -457,10 +460,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, amplicons = populate_amplicons(pair_info, primers); } - //this call is technically redundant but it doesn't stick when - //called in populate _amplicons - populate_pair_indices(primers, pair_info); - std::cout << "Amplicons detected: " << std::endl; amplicons.inOrder(); if(bam.empty()){ @@ -547,41 +546,68 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, uint32_t failed_frag_size = 0; uint32_t unmapped_counter = 0; uint32_t amplicon_flag_ctr = 0; + uint32_t pair_start = -1; //the start of the paired primer primer cand_primer; + primer pair; //the pair to the cand_primer std::vector overlapping_primers; std::vector::iterator cit; bool primer_trimmed = false; + + //redundant since this should have been called earlier, but somehow necessary? + populate_pair_indices(primers, pair_info); + //Iterate through reads while(sam_itr_next(in, iter, aln) >= 0) { unmapped_flag = false; primer_trimmed = false; + pair_start = -1; + get_overlapping_primers(aln, primers, overlapping_primers); + //print_all_primer_info(primers); + if((aln->core.flag&BAM_FUNMAP) == 0){ // If mapped // if primer pair info provided, check if read correctly overlaps with atleast one amplicon if(!pair_info.empty()){ amplicon_flag = amplicon_filter(amplicons, aln); if(!amplicon_flag){ - if (keep_for_reanalysis) { // -k (keep) option - aln->core.flag |= BAM_FQCFAIL; - if (bam_write1(out, aln) < 0) { retval = -1; goto error; } + if (keep_for_reanalysis) { // -k (keep) option + aln->core.flag |= BAM_FQCFAIL; + if (bam_write1(out, aln) < 0) { retval = -1; goto error; } } amplicon_flag_ctr++; continue; - } + } } - isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); - // if reverse strand + isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); + // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired get_overlapping_primers(aln, primers, overlapping_primers); - if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) + if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) primer_trimmed = true; - if(bam_is_rev(aln)){ // Reverse read - cand_primer = get_min_start(overlapping_primers); // fetch reverse primer (?) - - t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false); + if(bam_is_rev(aln)){ // Reverse read + cand_primer = get_min_start(overlapping_primers); // fetch reverse primer (?) + if(cand_primer.get_pair_indice() != -1){ + pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); + //we returned a non empty object + if(pair.get_start() != 0){ + pair_start = pair.get_end(); + } + } + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false, pair_start); } else { // Forward read cand_primer = get_max_end(overlapping_primers); // fetch forward primer (?) - t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); + if(cand_primer.get_pair_indice() != -1){ + pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); + //we returned a non empty object + if(pair.get_start() != 0){ + pair_start = pair.get_start(); + } + } + //TEST LINES + //print_primer_info(cand_primer); + //print_primer_info(pair); + + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false, pair_start); aln->core.pos += t.start_pos; } replace_cigar(aln, t.nlength, t.cigar); @@ -607,7 +633,14 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, if(overlapping_primers.size() > 0){ primer_trimmed = true; cand_primer = get_max_end(overlapping_primers); - t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); + if(cand_primer.get_pair_indice() != -1){ + pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); + //we returned a non empty object + if(pair.get_start() != 0){ + pair_start = pair.get_start(); + } + } + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false, pair_start); // Update read's left-most coordinate aln->core.pos += t.start_pos; replace_cigar(aln, t.nlength, t.cigar); @@ -621,7 +654,14 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, if(overlapping_primers.size() > 0){ primer_trimmed = true; cand_primer = get_min_start(overlapping_primers); - t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true); + if(cand_primer.get_pair_indice() != -1){ + pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); + //we returned a non empty object + if(pair.get_start() != 0){ + pair_start = pair.get_end(); + } + } + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true, pair_start); replace_cigar(aln, t.nlength, t.cigar); // Add count to primer cit = std::find(primers.begin(), primers.end(), cand_primer); From b5d75f58d25f98a4cd1b9c621b730cf16b947d3d Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 2 Aug 2022 08:21:51 -0700 Subject: [PATCH 08/22] reverting back the populate primers function and instead passing by reference --- src/interval_tree.cpp | 2 +- src/interval_tree.h | 2 +- src/primer_bed.cpp | 31 ++++++++++++++++++------------- src/primer_bed.h | 2 +- 4 files changed, 21 insertions(+), 16 deletions(-) diff --git a/src/interval_tree.cpp b/src/interval_tree.cpp index 1f1098cf..86e37f5f 100755 --- a/src/interval_tree.cpp +++ b/src/interval_tree.cpp @@ -83,7 +83,7 @@ void IntervalTree::inOrder(ITNode *root){ // A stand-alone function to create a tree containing the coordinates of each amplicon // based on user-specified primer pairs -IntervalTree populate_amplicons(std::string pair_info_file, std::vector primers){ +IntervalTree populate_amplicons(std::string pair_info_file, std::vector &primers){ int amplicon_start = -1; int amplicon_end = -1; IntervalTree tree = IntervalTree(); diff --git a/src/interval_tree.h b/src/interval_tree.h index 63742437..0ccbf3db 100755 --- a/src/interval_tree.h +++ b/src/interval_tree.h @@ -48,6 +48,6 @@ class IntervalTree{ void inOrder() {inOrder(_root);} }; -IntervalTree populate_amplicons(std::string pair_info_file, std::vector primers); +IntervalTree populate_amplicons(std::string pair_info_file, std::vector &primers); #endif diff --git a/src/primer_bed.cpp b/src/primer_bed.cpp index f53e541a..a1f0e315 100644 --- a/src/primer_bed.cpp +++ b/src/primer_bed.cpp @@ -254,7 +254,7 @@ int get_primer_indice(std::vector p, std::string name){ } //function using the tab seperated primer pair file -void populate_pair_indices(std::vector &primers, std::string path){ +int populate_pair_indices(std::vector &primers, std::string path){ /* * @param primers: the primer vector to add pair info to * @param path: the path to the primer pair file @@ -275,29 +275,34 @@ void populate_pair_indices(std::vector &primers, std::string path){ std::getline(line_stream, cell, '\t'); p2 = cell; line_stream.clear(); + + p1 = rtrim(p1); + p2 = rtrim(p2); + if(!p1.empty() && !p2.empty()){ for(it = primers.begin(); it != primers.end(); ++it) { - //search for primer name in pair file - if (it->get_name() == p1) { - //make sure it's pair exists - indice = get_primer_indice(primers, p2); - if (indice != -1){ - it->set_pair_indice(indice); - }else{ - std::cout << "Primer pair for " << p1 << " not found in BED file." <get_name() == p2){ - indice = get_primer_indice(primers, p1); + //search for primer name in pair file + if (it->get_name() == p1) { + //make sure it's pair exists + indice = get_primer_indice(primers, p2); + if (indice != -1){ + it->set_pair_indice(indice); + }else{ + std::cout << "Primer pair for " << p1 << " not found in BED file." <get_name() == p2){ + indice = get_primer_indice(primers, p1); if(indice != -1) it->set_pair_indice(indice); else std::cout << "Primer pair for " << p2 << " not found in BED file." << std::endl; - } + } } } else { std::cout << "Primer pair is empty." << std::endl; } } + return 0; } primer get_min_start(std::vector primers){ diff --git a/src/primer_bed.h b/src/primer_bed.h index 04a5d152..ef209b10 100644 --- a/src/primer_bed.h +++ b/src/primer_bed.h @@ -51,7 +51,7 @@ std::vector populate_from_file(std::string path); std::vector get_primers(std::vector p, unsigned int pos); int get_primer_indice(std::vector p, std::string name); primer fetch_primer_pair(int16_t index, std::vector primers); -void populate_pair_indices(std::vector &primers, std::string path); +int populate_pair_indices(std::vector &primers, std::string path); void print_primer_info(primer primers); void print_all_primer_info(std::vector primers); primer get_min_start(std::vector primers); From 7fc709ed581ab1c1db9290b9439d63b122d5769a Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 2 Aug 2022 08:50:30 -0700 Subject: [PATCH 09/22] saving for branch switch --- src/trim_primer_quality.cpp | 80 ++++++++++++++++++++++++------------- 1 file changed, 53 insertions(+), 27 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index ff85ab87..a736f2dc 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -2,13 +2,18 @@ #define round_int(x,total) ((int) (0.5 + ((float)x / float(total)) * 10000))/(float)100 +//called from primer_trim int32_t get_pos_on_query(uint32_t *cigar, uint32_t ncigar, int32_t pos, int32_t ref_start){ + /* + * @param cigar : cigar + * @param pos : starting position + */ int cig; int32_t n; int32_t ql = 0, rl = ref_start; for (uint32_t i = 0; i < ncigar; ++i){ - cig = bam_cigar_op(cigar[i]); - n = bam_cigar_oplen(cigar[i]); + cig = bam_cigar_op(cigar[i]); // cigar op + n = bam_cigar_oplen(cigar[i]); //cigar len if (bam_cigar_type(cig) & 2) { // Reference consuming if (pos <= rl + n) { if (bam_cigar_type(cig) & 1) // Query consuming @@ -164,15 +169,18 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r /* * @param r: alignment object * @param end_pos: the start of the opposite primer + * @param end_pos_2: the end of the opposite primer */ uint32_t *ncigar = (uint32_t*) malloc(sizeof(uint32_t) * (r->core.n_cigar + 1)), // Maximum edit is one more element with soft mask *cigar = bam_get_cigar(r); uint32_t i = 0, j = 0; int max_del_len = 0, cig, temp, del_len = 0; bool reverse = false; - - //std::cout << "start pos " << new_pos << " end pos " << end_pos << "\n"; - + std::cout << "start pos " << new_pos << " end pos " << end_pos << " r core " << r->core.pos << "\n"; + + //totality of what we're working with + std::cout << "cig length " << bam_cigar2qlen(r->core.n_cigar, bam_get_cigar(r)) << "\n"; + std::cout << "quer pos " << get_pos_on_query(cigar, r->core.n_cigar, new_pos, r->core.pos) << "\n"; if((r->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired and isize > read length if (bam_is_rev(r)){ // If -ve strand (?) max_del_len = bam_cigar2qlen(r->core.n_cigar, bam_get_cigar(r)) - get_pos_on_query(cigar, r->core.n_cigar, new_pos, r->core.pos) - 1; @@ -191,13 +199,32 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r } } max_del_len = (max_del_len > 0) ? max_del_len : 0; // For cases where reads spans only primer region - int32_t n, start_pos = 0, ref_add = 0; + int32_t n, start_pos = 0, ref_add = 0, ref_track = 0; bool pos_start = false; + bool soft_clip_all = false; del_len = max_del_len; - + ref_track = new_pos; //track absolute position + int direction = 0; + //if reverse the direction we count changes + if(reverse) direction = -1; + else direction = 1; + + //make it so that everytime we add a cigar, we add to start pos. if we exceed the start of the paired primer we chop + //iterate the cigar string while(i < r->core.n_cigar){ + if(soft_clip_all){ + ncigar[j] = bam_cigar_gen(bam_cigar_oplen(cigar[i]), BAM_CSOFT_CLIP); + i++; + j++; + continue; + } if (del_len == 0 && pos_start){ // No more bases on query to soft clip + ref_track += (bam_cigar_oplen(cigar[i]) * direction); ncigar[j] = cigar[i]; + std::cout << "ref track " << ref_track << " " << bam_cigar_op(cigar[i]) << " " << bam_cigar_oplen(cigar[i]) << std::endl; + if(ref_track > end_pos){ + soft_clip_all = true; + } i++; j++; continue; @@ -211,6 +238,7 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r //add the condition that we have a leading deletion if(cig == 2 && del_len == 0){ + ref_track += (bam_cigar_oplen(cigar[i]) * direction); ncigar[j] = cigar[i]; pos_start = true; i++; @@ -222,14 +250,17 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r if ((bam_cigar_type(cig) & 1)){ // Consumes Query if(del_len >= n ){ - ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); + ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); + ref_track += (n * direction); } else if (del_len < n && del_len > 0){ - ncigar[j] = bam_cigar_gen(del_len, BAM_CSOFT_CLIP); + ncigar[j] = bam_cigar_gen(del_len, BAM_CSOFT_CLIP); + ref_track += (del_len * direction); } else if (del_len == 0) { // Adding insertions before start position of read - ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); - j++; - i++; - continue; + ref_track += (n * direction); + ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); + j++; + i++; + continue; } j++; ref_add = std::min(del_len, n); @@ -237,8 +268,9 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r n = std::max(n - del_len, 0); del_len = std::max(del_len - temp, 0); if(n > 0){ - ncigar[j] = bam_cigar_gen(n, cig); - j++; + ref_track += (n * direction); + ncigar[j] = bam_cigar_gen(n, cig); + j++; } //add the condition that we have a leading deletion if(cig == 2 && del_len == 0){ @@ -256,18 +288,18 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r //deletions consume the reference but not the query, //insertions consume the query but not the reference if((bam_cigar_type(cig) & 2)) { // Consumes reference but not query - //std::cout << "ref add " << ref_add << " " << start_pos << " " << cig << "\n"; start_pos += ref_add; } i++; } - - /*uint32_t p=0; + + /* + uint32_t p=0; while(p < j){ std::cout << bam_cigar_op(ncigar[p]) << " " << bam_cigar_oplen(ncigar[p]) <<"\n"; p++; - }*/ - //std::cout << "start pos " << start_pos << "\n"; + } + */ if(reverse){ reverse_cigar(ncigar, j); } @@ -553,9 +585,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::vector::iterator cit; bool primer_trimmed = false; - //redundant since this should have been called earlier, but somehow necessary? - populate_pair_indices(primers, pair_info); - //Iterate through reads while(sam_itr_next(in, iter, aln) >= 0) { unmapped_flag = false; @@ -603,10 +632,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, pair_start = pair.get_start(); } } - //TEST LINES - //print_primer_info(cand_primer); - //print_primer_info(pair); - + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false, pair_start); aln->core.pos += t.start_pos; } From 6695b5675378b2fbbf7d3dbceccaea07d7fe2e12 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 5 Aug 2022 14:20:49 -0700 Subject: [PATCH 10/22] adding function to get overlapping primers for unpaired reads from pair info file --- src/trim_primer_quality.cpp | 205 ++++++++++++++++++------------------ 1 file changed, 101 insertions(+), 104 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index a736f2dc..c46db74f 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -165,22 +165,20 @@ void print_cigar(uint32_t *cigar, int nlength){ std::cout << std::endl; } -cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_rev = false, int32_t end_pos=-1){ +//the tricky function to change, edit with caution +cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_rev = false){ /* * @param r: alignment object - * @param end_pos: the start of the opposite primer - * @param end_pos_2: the end of the opposite primer + * @param new_pos: the start pos on query */ + uint32_t *ncigar = (uint32_t*) malloc(sizeof(uint32_t) * (r->core.n_cigar + 1)), // Maximum edit is one more element with soft mask *cigar = bam_get_cigar(r); uint32_t i = 0, j = 0; int max_del_len = 0, cig, temp, del_len = 0; bool reverse = false; - std::cout << "start pos " << new_pos << " end pos " << end_pos << " r core " << r->core.pos << "\n"; - - //totality of what we're working with - std::cout << "cig length " << bam_cigar2qlen(r->core.n_cigar, bam_get_cigar(r)) << "\n"; - std::cout << "quer pos " << get_pos_on_query(cigar, r->core.n_cigar, new_pos, r->core.pos) << "\n"; + //test lines + if((r->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired and isize > read length if (bam_is_rev(r)){ // If -ve strand (?) max_del_len = bam_cigar2qlen(r->core.n_cigar, bam_get_cigar(r)) - get_pos_on_query(cigar, r->core.n_cigar, new_pos, r->core.pos) - 1; @@ -199,32 +197,12 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r } } max_del_len = (max_del_len > 0) ? max_del_len : 0; // For cases where reads spans only primer region - int32_t n, start_pos = 0, ref_add = 0, ref_track = 0; + int32_t n, start_pos = 0, ref_add = 0; bool pos_start = false; - bool soft_clip_all = false; del_len = max_del_len; - ref_track = new_pos; //track absolute position - int direction = 0; - //if reverse the direction we count changes - if(reverse) direction = -1; - else direction = 1; - - //make it so that everytime we add a cigar, we add to start pos. if we exceed the start of the paired primer we chop - //iterate the cigar string while(i < r->core.n_cigar){ - if(soft_clip_all){ - ncigar[j] = bam_cigar_gen(bam_cigar_oplen(cigar[i]), BAM_CSOFT_CLIP); - i++; - j++; - continue; - } if (del_len == 0 && pos_start){ // No more bases on query to soft clip - ref_track += (bam_cigar_oplen(cigar[i]) * direction); ncigar[j] = cigar[i]; - std::cout << "ref track " << ref_track << " " << bam_cigar_op(cigar[i]) << " " << bam_cigar_oplen(cigar[i]) << std::endl; - if(ref_track > end_pos){ - soft_clip_all = true; - } i++; j++; continue; @@ -238,7 +216,6 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r //add the condition that we have a leading deletion if(cig == 2 && del_len == 0){ - ref_track += (bam_cigar_oplen(cigar[i]) * direction); ncigar[j] = cigar[i]; pos_start = true; i++; @@ -247,20 +224,18 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r } ref_add = n; - + if ((bam_cigar_type(cig) & 1)){ // Consumes Query if(del_len >= n ){ - ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); - ref_track += (n * direction); - } else if (del_len < n && del_len > 0){ - ncigar[j] = bam_cigar_gen(del_len, BAM_CSOFT_CLIP); - ref_track += (del_len * direction); + ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); + } + else if (del_len < n && del_len > 0){ //this doesn't work for reverse reads + ncigar[j] = bam_cigar_gen(del_len, BAM_CSOFT_CLIP); } else if (del_len == 0) { // Adding insertions before start position of read - ref_track += (n * direction); - ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); - j++; - i++; - continue; + ncigar[j] = bam_cigar_gen(n, BAM_CSOFT_CLIP); + j++; + i++; + continue; } j++; ref_add = std::min(del_len, n); @@ -268,12 +243,11 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r n = std::max(n - del_len, 0); del_len = std::max(del_len - temp, 0); if(n > 0){ - ref_track += (n * direction); - ncigar[j] = bam_cigar_gen(n, cig); - j++; + ncigar[j] = bam_cigar_gen(n, cig); + j++; } - //add the condition that we have a leading deletion - if(cig == 2 && del_len == 0){ + //add the condition that we have a leading deletion + if(cig == 2 && del_len == 0){ ncigar[j] = cigar[i]; pos_start = true; i++; @@ -292,14 +266,15 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r } i++; } - + + //test lines /* uint32_t p=0; while(p < j){ std::cout << bam_cigar_op(ncigar[p]) << " " << bam_cigar_oplen(ncigar[p]) <<"\n"; p++; - } - */ + }*/ + if(reverse){ reverse_cigar(ncigar, j); } @@ -312,6 +287,7 @@ cigar_ primer_trim(bam1_t *r, bool &isize_flag, int32_t new_pos, bool unpaired_r }; } + void replace_cigar(bam1_t *b, uint32_t n, uint32_t *cigar){ if (n != b->core.n_cigar) { int o = b->core.l_qname + b->core.n_cigar * 4; @@ -370,6 +346,7 @@ std::vector insertionSort(std::vector primers, uint32_t n){ return primers; } + // For paired reads void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector &overlapped_primers){ overlapped_primers.clear(); @@ -397,12 +374,48 @@ void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector //std::cout << it->get_start() << " " << start_pos << "\n"; overlapped_primers.push_back(*it); } - //print_primers(overlapped_primers); - //std::cout << "break\n"; } + +//function to get overlapping primers for unpaired reads, using the actual read direction +void get_overlapping_primers_better(bam1_t* r, std::vector primers, std::vector &overlapped_primers, bool unpaired_rev_gt, bool return_rev){ + /* + * @param unpaired_rev_gt : whether or not the read is reverse + * @param return_rev : whether we're looking for the forward or reverse primers + */ + overlapped_primers.clear(); + uint32_t start_pos = -1; + char strand = '+'; + if(unpaired_rev_gt){ + start_pos = bam_endpos(r) - 1; + strand = '-'; + } else { + start_pos = r->core.pos; + } + + //we either have a forward read && are looking for forward primers or we have a reverse read && are looking for reverse primers + if (unpaired_rev_gt == return_rev){ + for(std::vector::iterator it = primers.begin(); it != primers.end(); ++it) { + if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() ||it->get_strand() == 0)) + overlapped_primers.push_back(*it); + } + } else { + //in this case we need to invert what we find by indice + for(std::vector::iterator it = primers.begin(); it != primers.end(); ++it) { + if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() ||it->get_strand() == 0)){ + //if we have a paired primer + if(it->get_pair_indice() != -1){ + overlapped_primers.push_back(fetch_primer_pair(it->get_pair_indice(), primers)); + } + } + } + } +} + + // For unpaired reads void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector &overlapped_primers, bool unpaired_rev){ + std::cout << "get overlapping primers " << unpaired_rev << std::endl; overlapped_primers.clear(); uint32_t start_pos = -1; char strand = '+'; @@ -606,34 +619,21 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, amplicon_flag_ctr++; continue; } - } + } + isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired - get_overlapping_primers(aln, primers, overlapping_primers); + get_overlapping_primers(aln, primers, overlapping_primers); if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) primer_trimmed = true; if(bam_is_rev(aln)){ // Reverse read cand_primer = get_min_start(overlapping_primers); // fetch reverse primer (?) - if(cand_primer.get_pair_indice() != -1){ - pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); - //we returned a non empty object - if(pair.get_start() != 0){ - pair_start = pair.get_end(); - } - } - t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false, pair_start); + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false); + } else { // Forward read - cand_primer = get_max_end(overlapping_primers); // fetch forward primer (?) - if(cand_primer.get_pair_indice() != -1){ - pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); - //we returned a non empty object - if(pair.get_start() != 0){ - pair_start = pair.get_start(); - } - } - - t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false, pair_start); + cand_primer = get_max_end(overlapping_primers); // fetch forward primer (?) + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); aln->core.pos += t.start_pos; } replace_cigar(aln, t.nlength, t.cigar); @@ -653,53 +653,50 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, if(abs(aln->core.isize) <= abs(aln->core.l_qseq)){ failed_frag_size++; } - - // Forward primer - get_overlapping_primers(aln, primers, overlapping_primers, false); + + //handle the unpaired reads + //forward primer unpaired read + if (!pair_info.empty()){ + get_overlapping_primers_better(aln, primers, overlapping_primers, bam_is_rev(aln), false); + }else{ + get_overlapping_primers(aln, primers, overlapping_primers, false); + } if(overlapping_primers.size() > 0){ primer_trimmed = true; - cand_primer = get_max_end(overlapping_primers); - if(cand_primer.get_pair_indice() != -1){ - pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); - //we returned a non empty object - if(pair.get_start() != 0){ - pair_start = pair.get_start(); - } - } - t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false, pair_start); - // Update read's left-most coordinate + cand_primer = get_max_end(overlapping_primers); + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); + // Update read's left-most coordinate aln->core.pos += t.start_pos; replace_cigar(aln, t.nlength, t.cigar); // Add count to primer cit = std::find(primers.begin(), primers.end(), cand_primer); if(cit != primers.end()) - cit->add_read_count(1); - } - // Reverse primer - get_overlapping_primers(aln, primers, overlapping_primers, true); + cit->add_read_count(1); + } + + //reverse primer unpaired read + if (!pair_info.empty()){ + get_overlapping_primers_better(aln, primers, overlapping_primers, bam_is_rev(aln), true); + }else{ + get_overlapping_primers(aln, primers, overlapping_primers, true); + } if(overlapping_primers.size() > 0){ primer_trimmed = true; cand_primer = get_min_start(overlapping_primers); - if(cand_primer.get_pair_indice() != -1){ - pair = fetch_primer_pair(cand_primer.get_pair_indice(), primers); - //we returned a non empty object - if(pair.get_start() != 0){ - pair_start = pair.get_end(); - } - } - t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true, pair_start); - replace_cigar(aln, t.nlength, t.cigar); + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true); + replace_cigar(aln, t.nlength, t.cigar); // Add count to primer cit = std::find(primers.begin(), primers.end(), cand_primer); if(cit != primers.end()) cit->add_read_count(1); - } - t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming - if(bam_is_rev(aln)) // if reverse strand - aln->core.pos = t.start_pos; - condense_cigar(&t); - replace_cigar(aln, t.nlength, t.cigar); - } + } + t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming + if(bam_is_rev(aln)) // if reverse strand + aln->core.pos = t.start_pos; + condense_cigar(&t); + replace_cigar(aln, t.nlength, t.cigar); + } //end handling unpaired reads + if(primer_trimmed){ primer_trim_count++; } From 99ccc6ed74c5e1e6d2a935b3815fe09cf661d995 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Fri, 5 Aug 2022 14:30:33 -0700 Subject: [PATCH 11/22] removing extra unused variable --- src/trim_primer_quality.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index c46db74f..b341982f 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -587,13 +587,11 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, uint32_t primer_trim_count = 0, no_primer_counter = 0, low_quality = 0; bool unmapped_flag = false; bool amplicon_flag = false; + primer cand_primer; bool isize_flag = true; uint32_t failed_frag_size = 0; uint32_t unmapped_counter = 0; uint32_t amplicon_flag_ctr = 0; - uint32_t pair_start = -1; //the start of the paired primer - primer cand_primer; - primer pair; //the pair to the cand_primer std::vector overlapping_primers; std::vector::iterator cit; bool primer_trimmed = false; @@ -602,8 +600,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, while(sam_itr_next(in, iter, aln) >= 0) { unmapped_flag = false; primer_trimmed = false; - pair_start = -1; - + get_overlapping_primers(aln, primers, overlapping_primers); //print_all_primer_info(primers); From 29b6a307c1228dd96b5a6b0f5e3208fbd0aa0f21 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 30 Aug 2022 13:20:37 -0700 Subject: [PATCH 12/22] adding test data, template for test code --- data/primer_pair_test/test_primer_pair.bed | 17 +++++++++++++++ .../test_primer_pair_trim.tsv | 7 ++++++ .../test_primer_pair_trim.untrimmed.bam | Bin 0 -> 731 bytes tests/test_primer_trim_pairs.cpp | 20 ++++++++++++++++++ 4 files changed, 44 insertions(+) create mode 100755 data/primer_pair_test/test_primer_pair.bed create mode 100644 data/primer_pair_test/test_primer_pair_trim.tsv create mode 100644 data/primer_pair_test/test_primer_pair_trim.untrimmed.bam create mode 100755 tests/test_primer_trim_pairs.cpp diff --git a/data/primer_pair_test/test_primer_pair.bed b/data/primer_pair_test/test_primer_pair.bed new file mode 100755 index 00000000..8094d7ed --- /dev/null +++ b/data/primer_pair_test/test_primer_pair.bed @@ -0,0 +1,17 @@ +MN908947.3 3507 3531 nCoV-2019_11_RIGHT 1 - +MN908947.3 3460 3482 nCoV-2019_12_LEFT 2 + +MN908947.3 3826 3853 nCoV-2019_12_RIGHT 2 - +MN908947.3 3771 3795 nCoV-2019_13_LEFT 1 + +MN908947.3 4142 4164 nCoV-2019_13_RIGHT 1 - +MN908947.3 4054 4077 nCoV-2019_14_LEFT 2 + +MN908947.3 4044 4068 nCoV-2019_14_LEFT_alt4 2 + +MN908947.3 4428 4450 nCoV-2019_14_RIGHT 2 - +MN908947.3 4402 4424 nCoV-2019_14_RIGHT_alt2 2 - +MN908947.3 4294 4321 nCoV-2019_15_LEFT 1 + +MN908947.3 4295 4322 nCoV-2019_15_LEFT_alt1_ssi1 1 + +MN908947.3 4296 4322 nCoV-2019_15_LEFT_alt1 1 + +MN908947.3 4674 4696 nCoV-2019_15_RIGHT 1 - +MN908947.3 4666 4689 nCoV-2019_15_RIGHT_alt3 1 - +MN908947.3 4636 4658 nCoV-2019_16_LEFT 2 + +MN908947.3 4995 5017 nCoV-2019_16_RIGHT 2 - +MN908947.3 4939 4966 nCoV-2019_17_LEFT 1 + diff --git a/data/primer_pair_test/test_primer_pair_trim.tsv b/data/primer_pair_test/test_primer_pair_trim.tsv new file mode 100644 index 00000000..b861573a --- /dev/null +++ b/data/primer_pair_test/test_primer_pair_trim.tsv @@ -0,0 +1,7 @@ +nCoV-2019_12_LEFT nCoV-2019_12_RIGHT +nCoV-2019_13_LEFT nCoV-2019_13_RIGHT +nCoV-2019_14_LEFT nCoV-2019_14_RIGHT +nCoV-2019_14_LEFT_alt4 nCoV-2019_14_RIGHT_alt2 +nCoV-2019_15_LEFT nCoV-2019_15_RIGHT +nCoV-2019_15_LEFT_alt1 nCoV-2019_15_RIGHT_alt3 +nCoV-2019_16_LEFT nCoV-2019_16_RIGHT \ No newline at end of file diff --git a/data/primer_pair_test/test_primer_pair_trim.untrimmed.bam b/data/primer_pair_test/test_primer_pair_trim.untrimmed.bam new file mode 100644 index 0000000000000000000000000000000000000000..53b687af426a6ff3ce6de1ab7c58e4dcdc0ee2aa GIT binary patch literal 731 zcmV<10wnz(iwFb&00000{{{d;LjnN40=<<_kJ3OCz=!CCHK~+zO_N@V*#pU%-I;;X zcH)7OM#^Tv-6b011yP9EL|FpteM@@i(ZnxcvIno;)KB5}(0EYsj5DP`3IEz|;!7S( z-^@$q_x?@Es3txF0E(8Gd|D&u04D2?$z`X1(Hpc|JH`5GvR)(AnhrHRuO8%*$2Ed= z9p-k5XXRwYBu_fs=6SE(YV_L8R=?5g_uFq@zMNX?OzuA`xS8$qHg1(5T|otcP!1xH^0}|5 zZcVmpr>{2hY}yPQiNpPCiTfIq1Gca6xL<_hPDQBm4RUknS2O4lwXwvnfzDHCg%!C1 zS&^fHuEvGtrnhISX%LuE^1e7Pp=?KnO+4A0-YVc zu_=_qP1lkvbCS=-1rR4u;QRo|k{2PFl1DG^pFZGs3iyCRC+G-f*&X`=)@K@1P+H#v z=QKql3eMw0 +#include +#include "../src/trim_primer_quality.h" +#include "../src/interval_tree.h" +#include "../src/primer_bed.h" +#include "htslib/sam.h" +#include "../src/primer_bed.h" + +int main(){ + int num_tests = 1; + int success = 0; + std::vector::iterator it; + std::vector primers = populate_from_file("../data/primer_pair_test/test_primer_pair.bed"); + populate_pair_indices(primers, "../data/primer_pair_test/test_primer_pair_trim.tsv"); + std::string bam = "../data/primer/pair_test/test_primer_pair_trim.untrimmed.bam"; + + + //success += flag; + return (num_tests == success) ? 0 : -1; +} From 3c8792ba3494c4aeef956c98266f042b041633b7 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 31 Aug 2022 16:27:08 -0700 Subject: [PATCH 13/22] updating the primer pair test the call trimming --- tests/test_primer_trim_pairs.cpp | 110 ++++++++++++++++++++++++++++++- 1 file changed, 107 insertions(+), 3 deletions(-) diff --git a/tests/test_primer_trim_pairs.cpp b/tests/test_primer_trim_pairs.cpp index ea80429b..9129e238 100755 --- a/tests/test_primer_trim_pairs.cpp +++ b/tests/test_primer_trim_pairs.cpp @@ -9,12 +9,116 @@ int main(){ int num_tests = 1; int success = 0; - std::vector::iterator it; - std::vector primers = populate_from_file("../data/primer_pair_test/test_primer_pair.bed"); - populate_pair_indices(primers, "../data/primer_pair_test/test_primer_pair_trim.tsv"); std::string bam = "../data/primer/pair_test/test_primer_pair_trim.untrimmed.bam"; + std::string pair_info = "../data/primer_pair_test/test_primer_pair_trim.tsv" + //populate primers and amplicons from pair file + std::vector primers = populate_from_file("../data/primer_pair_test/test_primer_pair.bed"); + int max_primer_len = 0; + max_primer_len = get_bigger_primer(primers); + std::string region_; + //open bam file and build and index if needed + 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; + } + } + } + //read in header + 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; + //all this needs to be modified + int primer_indices[] = {0, 0, 7, 7, 6}; + uint8_t cigar_flag[5][11] = { + {BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CDEL, BAM_CPAD, BAM_CDEL, BAM_CMATCH}, + {BAM_CMATCH, BAM_CPAD, BAM_CDEL, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP}, + {BAM_CMATCH, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP, BAM_CSOFT_CLIP}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH} + }; + uint32_t read_start_pos[5] = { + 30, 29, 231, 249, 274 //also change + }; + uint8_t condense_cigar_flag[5][8] = { + {BAM_CSOFT_CLIP, BAM_CMATCH}, + {BAM_CSOFT_CLIP, BAM_CDEL, BAM_CPAD, BAM_CDEL, BAM_CMATCH}, //change this + {BAM_CMATCH, BAM_CPAD, BAM_CDEL, BAM_CSOFT_CLIP}, //change this + {BAM_CMATCH, BAM_CSOFT_CLIP}, + {BAM_CSOFT_CLIP, BAM_CMATCH, BAM_CDEL, BAM_CMATCH, BAM_CPAD, BAM_CMATCH, BAM_CPAD, BAM_CMATCH} + }; + uint32_t condense_cigar_len[5][8] = {{12, 139}, {33, 1,1,1,114}, {121, 1, 4, 23}, {103, 49}, {23, 18,3, 57, 1, 10, 1, 39}}; + unsigned int overlapping_primer_sizes[] = {0, 2, 2, 0, 0, 0, 0, 2, 2, 1}; + int ctr = 0; + std::vector overlapping_primers; + primer cand_primer; + bool isize_flag = false; + while(sam_itr_next(in, iter, aln) >= 0) { + if((aln->core.flag&BAM_FUNMAP) != 0){ + continue; + } + isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); + std::cout << bam_get_qname(aln) << std::endl; + get_overlapping_primers(aln, primers, overlapping_primers); + if(overlapping_primers.size() != overlapping_primer_sizes[ctr]){ + success = -1; + std::cout << "Overlapping primer sizes for " << bam_get_qname(aln) <<". Expected: " << overlapping_primer_sizes[ctr] << ". Got: " << overlapping_primers.size() << std::endl; + } + if(overlapping_primers.size() > 0){ + print_cigar(bam_get_cigar(aln), aln->core.n_cigar); + if(bam_is_rev(aln)){ + cand_primer = get_min_start(overlapping_primers); + t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, false); + } else { + cand_primer = get_max_end(overlapping_primers); + t = primer_trim(aln, isize_flag, cand_primer.get_end() + 1, false); + aln->core.pos += t.start_pos; + } + if(cand_primer.get_indice() != primer_indices[primer_ctr]){ + success = -1; + std::cout << "Primer indice wrong. Expected: " << primer_indices[primer_ctr] << ". Got: " << cand_primer.get_indice() << std::endl; + } + // Check start pos + if(aln->core.pos != (int) read_start_pos[primer_ctr]){ + success = -1; + std::cout << "Incorrect start position" << std::endl; + } + // Replace cigar + 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]){ + std::cout << bam_cigar_op(cigar[i]) << " " << bam_cigar_oplen(cigar[i]) << " i " << i << "\n"; + success = -1; + std::cout << "Cigar flag didn't match for " << cand_primer.get_indice() << " ! 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; + } + } //success += flag; return (num_tests == success) ? 0 : -1; + + //clear memory and close files + bam_destroy1(aln); + bam_itr_destroy(iter); + sam_hdr_destroy(header); + hts_idx_destroy(idx); + hts_close(in); } From 80116220fd57490d524e3016f43d95e1130f4a00 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Mon, 23 Jan 2023 12:05:58 -0800 Subject: [PATCH 14/22] removing extra function trim_primers_better and switching the way isize is calculated --- src/trim_primer_quality.cpp | 67 +++++++------------------------------ 1 file changed, 13 insertions(+), 54 deletions(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index b341982f..339c694e 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -359,7 +359,6 @@ void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector start_pos = r->core.pos; } - //print_primers(primers); //sort it first std::vector test = insertionSort(primers, primers.size()); //std::cout << test.size() << "\n"; @@ -377,45 +376,9 @@ void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector } -//function to get overlapping primers for unpaired reads, using the actual read direction -void get_overlapping_primers_better(bam1_t* r, std::vector primers, std::vector &overlapped_primers, bool unpaired_rev_gt, bool return_rev){ - /* - * @param unpaired_rev_gt : whether or not the read is reverse - * @param return_rev : whether we're looking for the forward or reverse primers - */ - overlapped_primers.clear(); - uint32_t start_pos = -1; - char strand = '+'; - if(unpaired_rev_gt){ - start_pos = bam_endpos(r) - 1; - strand = '-'; - } else { - start_pos = r->core.pos; - } - - //we either have a forward read && are looking for forward primers or we have a reverse read && are looking for reverse primers - if (unpaired_rev_gt == return_rev){ - for(std::vector::iterator it = primers.begin(); it != primers.end(); ++it) { - if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() ||it->get_strand() == 0)) - overlapped_primers.push_back(*it); - } - } else { - //in this case we need to invert what we find by indice - for(std::vector::iterator it = primers.begin(); it != primers.end(); ++it) { - if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() ||it->get_strand() == 0)){ - //if we have a paired primer - if(it->get_pair_indice() != -1){ - overlapped_primers.push_back(fetch_primer_pair(it->get_pair_indice(), primers)); - } - } - } - } -} - // For unpaired reads void get_overlapping_primers(bam1_t* r, std::vector primers, std::vector &overlapped_primers, bool unpaired_rev){ - std::cout << "get overlapping primers " << unpaired_rev << std::endl; overlapped_primers.clear(); uint32_t start_pos = -1; char strand = '+'; @@ -595,15 +558,15 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::vector overlapping_primers; std::vector::iterator cit; bool primer_trimmed = false; + std::string test = "NB552570:188:HL75JAFX3:2:21105:7297:5549"; //Iterate through reads while(sam_itr_next(in, iter, aln) >= 0) { unmapped_flag = false; primer_trimmed = false; - + get_overlapping_primers(aln, primers, overlapping_primers); - //print_all_primer_info(primers); - + if((aln->core.flag&BAM_FUNMAP) == 0){ // If mapped // if primer pair info provided, check if read correctly overlaps with atleast one amplicon if(!pair_info.empty()){ @@ -617,11 +580,14 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, continue; } } - + + //isize is insert size + //l_qseq is the query length isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); + isize_flag = (abs(aln->core.isize) - max_primer_len) > 0; // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired - get_overlapping_primers(aln, primers, overlapping_primers); + get_overlapping_primers(aln, primers, overlapping_primers); if(overlapping_primers.size() > 0){ // If read starts before overlapping regions (?) primer_trimmed = true; if(bam_is_rev(aln)){ // Reverse read @@ -651,13 +617,9 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, failed_frag_size++; } - //handle the unpaired reads + //handle the unpaired reads //forward primer unpaired read - if (!pair_info.empty()){ - get_overlapping_primers_better(aln, primers, overlapping_primers, bam_is_rev(aln), false); - }else{ - get_overlapping_primers(aln, primers, overlapping_primers, false); - } + get_overlapping_primers(aln, primers, overlapping_primers, false); if(overlapping_primers.size() > 0){ primer_trimmed = true; cand_primer = get_max_end(overlapping_primers); @@ -672,12 +634,9 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, } //reverse primer unpaired read - if (!pair_info.empty()){ - get_overlapping_primers_better(aln, primers, overlapping_primers, bam_is_rev(aln), true); - }else{ - get_overlapping_primers(aln, primers, overlapping_primers, true); - } - if(overlapping_primers.size() > 0){ + get_overlapping_primers(aln, primers, overlapping_primers, true); + + if(overlapping_primers.size() > 0){ primer_trimmed = true; cand_primer = get_min_start(overlapping_primers); t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true); From 75eb684569b9718bd837eb13d873adbf04b6e11a Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 24 Jan 2023 15:54:42 -0800 Subject: [PATCH 15/22] revert to the original way isize_flag is calculated --- src/trim_primer_quality.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index 339c694e..df13fbd3 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -584,7 +584,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, //isize is insert size //l_qseq is the query length isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); - isize_flag = (abs(aln->core.isize) - max_primer_len) > 0; // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired get_overlapping_primers(aln, primers, overlapping_primers); From 69f0ed78a28c33e0cebd40177bcf08acb75a8984 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 24 Jan 2023 16:22:40 -0800 Subject: [PATCH 16/22] changing loc of declaration of variable for github action macos --- src/call_consensus_pileup.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/call_consensus_pileup.cpp b/src/call_consensus_pileup.cpp index 192f1563..528acc04 100644 --- a/src/call_consensus_pileup.cpp +++ b/src/call_consensus_pileup.cpp @@ -32,7 +32,7 @@ ret_t get_consensus_allele(std::vector ad, uint8_t min_qual, double thre std::vector nuc_pos; allele tmp_a; char n; - uint32_t max_l = 0, max_depth = 0, tmp_depth = 0, cur_depth = 0, total_max_depth = 0, gap_depth = 0, total_indel_depth = 0; + uint32_t max_l = 0, max_depth = 0, tmp_depth = 0, total_max_depth = 0, gap_depth = 0, total_indel_depth = 0; uint8_t ambg_n = 1, ctr = 0; double q = 0, tq = 0, cur_threshold = 0; std::vector::iterator it; @@ -54,7 +54,7 @@ ret_t get_consensus_allele(std::vector ad, uint8_t min_qual, double thre tq = 0; max_depth = 0; tmp_depth = 0; - cur_depth = 0; + uint32_t cur_depth = 0; // prev_depth = 0; ctr = 1; it = ad.begin(); From f823ae1d643e3bd99ce5447ea4fa8217ae2f798d Mon Sep 17 00:00:00 2001 From: cmaceves Date: Tue, 24 Jan 2023 16:27:25 -0800 Subject: [PATCH 17/22] changing the declaration of cur_depth back, suppressing warning in compiler --- src/Makefile.am | 2 +- src/call_consensus_pileup.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index 0b55a8db..ab8c4508 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ LIBS = -lhts -lz -lpthread -CXXFLAGS = -v -g -std=c++11 -Wall -Wextra -Werror +CXXFLAGS = -v -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-parameter # this lists the binaries to produce, the (non-PHONY, binary) targets in # the previous manual Makefile diff --git a/src/call_consensus_pileup.cpp b/src/call_consensus_pileup.cpp index 528acc04..83619caa 100644 --- a/src/call_consensus_pileup.cpp +++ b/src/call_consensus_pileup.cpp @@ -32,7 +32,7 @@ ret_t get_consensus_allele(std::vector ad, uint8_t min_qual, double thre std::vector nuc_pos; allele tmp_a; char n; - uint32_t max_l = 0, max_depth = 0, tmp_depth = 0, total_max_depth = 0, gap_depth = 0, total_indel_depth = 0; + uint32_t max_l = 0, max_depth = 0, cur_depth = 0, tmp_depth = 0, total_max_depth = 0, gap_depth = 0, total_indel_depth = 0; uint8_t ambg_n = 1, ctr = 0; double q = 0, tq = 0, cur_threshold = 0; std::vector::iterator it; @@ -54,7 +54,7 @@ ret_t get_consensus_allele(std::vector ad, uint8_t min_qual, double thre tq = 0; max_depth = 0; tmp_depth = 0; - uint32_t cur_depth = 0; + cur_depth = 0; // prev_depth = 0; ctr = 1; it = ad.begin(); From 73fb315f42f98bb59d1258574a9c530fa54700f8 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 25 Jan 2023 13:02:02 -0800 Subject: [PATCH 18/22] Revert "revert to the original way isize_flag is calculated" This reverts commit 75eb684569b9718bd837eb13d873adbf04b6e11a. --- src/trim_primer_quality.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index df13fbd3..339c694e 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -584,6 +584,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, //isize is insert size //l_qseq is the query length isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); + isize_flag = (abs(aln->core.isize) - max_primer_len) > 0; // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired get_overlapping_primers(aln, primers, overlapping_primers); From 983d233cf99f17b2edfd6bbe3832fe40cf70005f Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 25 Jan 2023 13:06:45 -0800 Subject: [PATCH 19/22] attempt to silence specific warning on macos build --- src/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.am b/src/Makefile.am index ab8c4508..7f075282 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ LIBS = -lhts -lz -lpthread -CXXFLAGS = -v -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-parameter +CXXFLAGS = -v -g -std=c++11 -Wall -Wno-unused-but-set-variable -Wextra -Werror # this lists the binaries to produce, the (non-PHONY, binary) targets in # the previous manual Makefile From 66d4d9f5f1e274a7db2cedf8d0c4ea61f9cee963 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 25 Jan 2023 13:08:11 -0800 Subject: [PATCH 20/22] going back to the original way isize was called --- src/trim_primer_quality.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/trim_primer_quality.cpp b/src/trim_primer_quality.cpp index 339c694e..df13fbd3 100755 --- a/src/trim_primer_quality.cpp +++ b/src/trim_primer_quality.cpp @@ -584,7 +584,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, //isize is insert size //l_qseq is the query length isize_flag = (abs(aln->core.isize) - max_primer_len) > abs(aln->core.l_qseq); - isize_flag = (abs(aln->core.isize) - max_primer_len) > 0; // if reverse strand if((aln->core.flag&BAM_FPAIRED) != 0 && isize_flag){ // If paired get_overlapping_primers(aln, primers, overlapping_primers); From 1c8ef6ce250aa4f858ab7b5029cf660f5834fcc7 Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 25 Jan 2023 13:15:49 -0800 Subject: [PATCH 21/22] silencing the unused variable warning --- src/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.am b/src/Makefile.am index 7f075282..418fc7ed 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ LIBS = -lhts -lz -lpthread -CXXFLAGS = -v -g -std=c++11 -Wall -Wno-unused-but-set-variable -Wextra -Werror +CXXFLAGS = -v -g -std=c++11 -Wall -Wno-unused-variable -Wextra -Werror # this lists the binaries to produce, the (non-PHONY, binary) targets in # the previous manual Makefile From cbac14b5cb440093be5e47710f23707eb7b147dc Mon Sep 17 00:00:00 2001 From: cmaceves Date: Wed, 25 Jan 2023 13:27:40 -0800 Subject: [PATCH 22/22] changing workflow to use stable mac version and removing trial flag from Makefile that was meant to correct compile error --- .github/workflows/ccpp.yml | 2 +- src/Makefile.am | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml index 954f063f..b6b707ad 100644 --- a/.github/workflows/ccpp.yml +++ b/.github/workflows/ccpp.yml @@ -57,7 +57,7 @@ jobs: build-macos: - runs-on: [ macos-latest ] + runs-on: [ macos-11 ] steps: - uses: actions/checkout@master diff --git a/src/Makefile.am b/src/Makefile.am index 418fc7ed..0b55a8db 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ LIBS = -lhts -lz -lpthread -CXXFLAGS = -v -g -std=c++11 -Wall -Wno-unused-variable -Wextra -Werror +CXXFLAGS = -v -g -std=c++11 -Wall -Wextra -Werror # this lists the binaries to produce, the (non-PHONY, binary) targets in # the previous manual Makefile