remove_reads_from_amplicon.cpp
#include "primer_bed.h"
#include "htslib/sam.h"
#include "htslib/bgzf.h"
#include<iostream>
int rmv_reads_from_amplicon(std::string bam, std::string region_, std::string bam_out, uint16_t* amplicon, int amp_n){
for (int i = 0; i < amp_n; ++i){
std::cout << "Amplicon: " << amplicon[i] << " ";
}
bam_out += ".bam";
std::cout << std::endl;
std::cout << bam_out << std::endl;
if(bam.empty()){
std::cout << "Bam in empty" << std::endl;
return 0;
}
//open BAM for reading
samFile *in = hts_open(bam.c_str(), "r");
BGZF *out = bgzf_open(bam_out.c_str(), "w");
if(in == NULL) {
std::cout << ("Unable to open BAM/SAM file.") << std::endl;
return -1;
}
//Load the index
hts_idx_t *idx = sam_index_load(in, bam.c_str());
if(idx == NULL) {
std::cout << ("Unable to open BAM/SAM index.") << std::endl; // TODO: Generate index
return -1;
}
//Get the header
bam_hdr_t *header = sam_hdr_read(in);
if(bam_hdr_write(out, header) < 0){
std::cout << "Unable to write BAM header to path." << std::endl;
sam_close(in);
return -1;
}
if(header == NULL) {
sam_close(in);
std::cout << "Unable to open BAM/SAM header." << std::endl;
}
if (region_.empty()){
std::cout << "Number of references: " << header->n_targets << std::endl;
for (int i = 0; i < header->n_targets; ++i){
std::cout << "Reference Name: " << header->target_name[i] << std::endl;
std::cout << "Reference Length: " << header->target_len[i] << std::endl;
if(i==0){
region_.assign(header->target_name[i]);
}
}
std::cout << "Using Region: " << region_ << std::endl;
}
std::string temp(header->text);
std::string sortFlag ("SO:coordinate");
if (temp.find(sortFlag)) {
std::cout << "Sorted By Coordinate" << std::endl; // Sort by coordinate
} else {
std::cout << "Not sorted" << std::endl;
}
//Initialize iterator
hts_itr_t *iter = NULL;
//Move the iterator to the region we are interested in
iter = sam_itr_querys(idx, header, region_.c_str());
if(header == NULL || iter == NULL) {
sam_close(in);
std::cout << "Unable to iterate to region within BAM/SAM." << std::endl;
return -1;
}
//Initiate the alignment record
bam1_t *aln = bam_init1();
int ctr = 0;
bool w;
while(sam_itr_next(in, iter, aln) >= 0) {
uint8_t* a = bam_aux_get(aln, "xa");
w = true;
if(a != 0){
for (int i = 0; i < amp_n; ++i){
if(bam_aux2i(a) == amplicon[i]){
w = false;
}
}
}
if(w){
if(bam_write1(out, aln) < 0){
std::cout << "Not able to write to BAM" << std::endl;
hts_itr_destroy(iter);
hts_idx_destroy(idx);
bam_destroy1(aln);
bam_hdr_destroy(header);
sam_close(in);
bgzf_close(out);
return -1;
};
}
ctr++;
if(ctr % 100000 == 0){
std::cout << ctr << std::endl;
}
}
hts_itr_destroy(iter);
hts_idx_destroy(idx);
bam_destroy1(aln);
bam_hdr_destroy(header);
sam_close(in);
bgzf_close(out);
return 0;
}
// int main_old(int argc, char* argv[]){
// std::string bam = std::string(argv[1]);
// std::string region_;
// std::string bam_out = std::string(argv[2]);
// argc -= 3;
// std::cout << argc << std::endl;
// int16_t amplicon[argc];
// for (int i = 0; i < argc; ++i){
// amplicon[i] = atoi(argv[i+3]); // TODO: Convert to int16_t
// }
// for (int i = 0; i < sizeof(amplicon)/sizeof(*amplicon); ++i){
// std::cout << "Amplicon: " << amplicon[i] << " ";
// }
// std::cout << std::endl;
// std::cout << bam_out << std::endl;
// if(!bam.empty()) {
// //open BAM for reading
// samFile *in = hts_open(bam.c_str(), "r");
// BGZF *out = bgzf_open(bam_out.c_str(), "w");
// if(in == NULL) {
// throw std::runtime_error("Unable to open BAM/SAM file.");
// }
// //Load the index
// hts_idx_t *idx = sam_index_load(in, bam.c_str());
// if(idx == NULL) {
// throw std::runtime_error("Unable to open BAM/SAM index."); // Generate index
// }
// //Get the header
// bam_hdr_t *header = sam_hdr_read(in);
// bam_hdr_write(out, header);
// if(header == NULL) {
// sam_close(in);
// throw std::runtime_error("Unable to open BAM header.");
// }
// if (region_.empty()){
// std::cout << "Number of references: " << header->n_targets << std::endl;
// for (int i = 0; i < header->n_targets; ++i){
// std::cout << "Reference Name: " << header->target_name[i] << std::endl;
// std::cout << "Reference Length: " << header->target_len[i] << std::endl;
// if(i==0){
// region_.assign(header->target_name[i]);
// }
// }
// std::cout << "Using Region: " << region_ << std::endl;
// }
// std::string temp(header->text);
// std::string sortFlag ("SO:coordinate");
// if (temp.find(sortFlag)) {
// std::cout << "Sorted By Coordinate" << std::endl; // Sort by coordinate
// } else {
// std::cout << "Not sorted" << std::endl;
// }
// //Initialize iterator
// hts_itr_t *iter = NULL;
// //Move the iterator to the region we are interested in
// iter = sam_itr_querys(idx, header, region_.c_str());
// if(header == NULL || iter == NULL) {
// sam_close(in);
// throw std::runtime_error("Unable to iterate to region within BAM.");
// }
// //Initiate the alignment record
// bam1_t *aln = bam_init1();
// int ctr = 0;
// bool w;
// while(sam_itr_next(in, iter, aln) >= 0) {
// // std::cout << "Query Length: " << bam_cigar2qlen(aln->core.n_cigar, cigar) << std::endl;
// // std::cout << "Read Length: " << bam_cigar2rlen(aln->core.n_cigar, cigar) << std::endl;
// // std::cout << "Query Start: " << get_query_start(aln) << std::endl;
// // std::cout << "Query End: " << get_query_end(aln) << std::endl
// uint8_t* a = bam_aux_get(aln, "xa");
// w = true;
// for (int i = 0; i < sizeof(amplicon)/sizeof(*amplicon); ++i){
// if(bam_aux2i(a) == amplicon[i]){
// w = false;
// }
// }
// if(w)
// bam_write1(out, aln);
// ctr++;
// if(ctr % 100000 == 0){
// std::cout << ctr << std::endl;
// }
// }
// hts_itr_destroy(iter);
// hts_idx_destroy(idx);
// bam_destroy1(aln);
// bam_hdr_destroy(header);
// sam_close(in);
// bgzf_close(out);
// }
// return 0;
// }