feat: added per-read flag and map quality filtering

 - closes griffithlab#183, flag based filtering
 - closes griffithlab#189, mapping quality based filtering
 - option '-F' filters reads containing any of these flags
 - option '-f' filters reads not containing all these flags
 - option '-q' filters reads below this mapping quality

Author: TimD1

src/cis-splice-effects/cis_splice_effects_identifier.cc

@@ -46,6 +46,9 @@ void CisSpliceEffectsIdentifier::usage(ostream& out) {
         << "\t\t\t " << "anchor length on both sides are reported. [8]" << endl;
     out << "\t\t" << "-m INT\tMinimum intron length. [70]" << endl;
     out << "\t\t" << "-M INT\tMaximum intron length. [500000]" << endl;
+    out << "\t\t" << "-f INT\tOnly use alignments where all flag bits set here are set. [0]" << endl;
+    out << "\t\t" << "-F INT\tOnly use alignments where no flag bits set here are set. [0]" << endl;
+    out << "\t\t" << "-q INT\tOnly use alignments with this mapping quality or above. [0]" << endl;
     out << "\t\t" << "-w INT\tWindow size in b.p to identify splicing events in.\n" 
         << "\t\t\t " << "The tool identifies events in variant.start +/- w basepairs.\n"
         << "\t\t\t " << "Default behaviour is to look at the window between previous and next exons." << endl;
@@ -113,7 +116,7 @@ void CisSpliceEffectsIdentifier::parse_options(int argc, char* argv[]) {
     optind = 1; //Reset before parsing again.
     stringstream help_ss;
     char c;
-    while((c = getopt(argc, argv, "o:w:v:j:e:Ei:ISht:s:a:m:M:b:C")) != -1) {
+    while((c = getopt(argc, argv, "o:w:v:j:e:Ei:ISht:s:a:m:M:f:F:q:b:C")) != -1) {
         switch(c) {
             case 'o':
                 output_file_ = string(optarg);
@@ -170,6 +173,15 @@ void CisSpliceEffectsIdentifier::parse_options(int argc, char* argv[]) {
             case 'M':
                 max_intron_length_ = atoi(optarg);
                 break;
+            case 'f':
+                require_flags_ = atoi(optarg);
+                break;
+            case 'F':
+                filter_flags_ = atoi(optarg);
+                break;
+            case 'q':
+                min_map_qual_ = atoi(optarg);
+                break;
             case 'b':
                 output_barcodes_file_ = string(optarg);
                 break;
@@ -285,7 +297,9 @@ void CisSpliceEffectsIdentifier::identify() {
             } else {
                 ref_to_pass = "NA";
             }
-            JunctionsExtractor je1(bam_, variant_region, strandness_, strand_tag_, min_anchor_length_, min_intron_length_, max_intron_length_, ref_to_pass);
+            JunctionsExtractor je1(bam_, variant_region, strandness_, 
+                    strand_tag_, min_anchor_length_, min_intron_length_, max_intron_length_, 
+                    filter_flags_, require_flags_, min_map_qual_, ref_to_pass);
             je1.identify_junctions_from_BAM();
             vector<Junction> junctions = je1.get_all_junctions();
             //Add all the junctions to the unique set

src/cis-splice-effects/cis_splice_effects_identifier.h

@@ -87,13 +87,19 @@ class CisSpliceEffectsIdentifier {
         //tag used in BAM to denote strand, default "XS"
         string strand_tag_;
         //Minimum anchor length for junctions
-        //Junctions need atleast this many bp overlap
+        //Junctions need at least this many bp overlap
         // on both ends.
         uint32_t min_anchor_length_;
         //Minimum length of an intron, i.e min junction width
         uint32_t min_intron_length_;
         //Maximum length of an intron, i.e max junction width
         uint32_t max_intron_length_;
+        //filter reads containing any of these flags
+        uint16_t filter_flags_;
+        // filter reads not containing all of these flags
+        uint16_t require_flags_;
+        // filter reads below the minimum mapping quality
+        uint8_t min_map_qual_;
         //whether to override strand of extracted junctions using intron-motif method
         bool override_strand_with_canonical_intron_motif_;
     public:
@@ -114,6 +120,9 @@ class CisSpliceEffectsIdentifier {
                                        min_anchor_length_(8),
                                        min_intron_length_(70),
                                        max_intron_length_(500000),
+                                       filter_flags_(0),
+                                       require_flags_(0),
+                                       min_map_qual_(0),
                                        override_strand_with_canonical_intron_motif_(false) {}
         //Destructor
         ~CisSpliceEffectsIdentifier() {

src/junctions/junctions_extractor.cc

@@ -43,7 +43,7 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
     optind = 1; //Reset before parsing again.
     int c;
     stringstream help_ss;
-    while((c = getopt(argc, argv, "ha:m:M:o:r:t:s:b:")) != -1) {
+    while((c = getopt(argc, argv, "ha:m:M:f:F:q:o:r:t:s:b:")) != -1) {
         switch(c) {
             case 'h':
                 usage(help_ss);
@@ -57,6 +57,15 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
             case 'M':
                 max_intron_length_ = atoi(optarg);
                 break;
+            case 'f':
+                require_flags_ = atoi(optarg);
+                break;
+            case 'F':
+                filter_flags_ = atoi(optarg);
+                break;
+            case 'q':
+                min_map_qual_ = atoi(optarg);
+                break;
             case 'o':
                 output_file_ = string(optarg);
                 break;
@@ -108,10 +117,17 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
             throw runtime_error("Strandness mode 'intron-motif' requires a fasta file!\n\n");
         }
     }
+    if ( (require_flags_ & filter_flags_) != 0) {
+        usage();
+        throw runtime_error("No overlap allowed between '-f' and '-F' options (same flag filtered and required)!\n\n");
+    }
 
     cerr << "Minimum junction anchor length: " << min_anchor_length_ << endl;
     cerr << "Minimum intron length: " << min_intron_length_ << endl;
     cerr << "Maximum intron length: " << max_intron_length_ << endl;
+    cerr << "Require alignment flags: " << require_flags_ << endl;
+    cerr << "Filter alignment flags: " << filter_flags_ << endl;
+    cerr << "Minimum alignment mapping quality: " << int(min_map_qual_) << endl;
     cerr << "Alignment: " << bam_ << endl;
     cerr << "Output file: " << output_file_ << endl;
     if (output_barcodes_file_ != "NA"){
@@ -130,6 +146,9 @@ int JunctionsExtractor::usage(ostream& out) {
         << "\t\t\t " << "anchor length on both sides are reported. [8]" << endl;
     out << "\t\t" << "-m INT\tMinimum intron length. [70]" << endl;
     out << "\t\t" << "-M INT\tMaximum intron length. [500000]" << endl;
+    out << "\t\t" << "-f INT\tOnly use alignments where all flag bits set here are set. [0]" << endl;
+    out << "\t\t" << "-F INT\tOnly use alignments where no flag bits set here are set. [0]" << endl;
+    out << "\t\t" << "-q INT\tOnly use alignments with this mapping quality or above. [0]" << endl;
     out << "\t\t" << "-o FILE\tThe file to write output to. [STDOUT]" << endl;
     out << "\t\t" << "-r STR\tThe region to identify junctions \n"
         << "\t\t\t " << "in \"chr:start-end\" format. Entire BAM by default." << endl;
@@ -358,7 +377,7 @@ void JunctionsExtractor::set_junction_strand(bam1_t *aln, Junction& j1, string i
     return;
 }
 
-//Get the the barcode
+//Get the barcode
 void JunctionsExtractor::set_junction_barcode(bam1_t *aln, Junction& j1) {
     uint8_t *p = bam_aux_get(aln, barcode_tag_.c_str());
     if(p != NULL) {
@@ -373,6 +392,14 @@ void JunctionsExtractor::set_junction_barcode(bam1_t *aln, Junction& j1) {
     }
 }
 
+//Filters alignments based on filtering flags and mapping quality
+bool JunctionsExtractor::filter_alignment(bam_hdr_t *header, bam1_t *aln) {
+    if ((aln->core.flag & filter_flags_) != 0) return true;
+    if ((aln->core.flag | require_flags_) != aln->core.flag) return true;
+    if (aln->core.qual < min_map_qual_) return true;
+    return false;
+}
+
 //Parse junctions from the read and store in junction map
 int JunctionsExtractor::parse_alignment_into_junctions(bam_hdr_t *header, bam1_t *aln) {
     int n_cigar = aln->core.n_cigar;
@@ -523,6 +550,7 @@ int JunctionsExtractor::identify_junctions_from_BAM() {
         //Initiate the alignment record
         bam1_t *aln = bam_init1();
         while(sam_itr_next(in, iter, aln) >= 0) {
+            if (filter_alignment(header, aln)) continue;
             parse_alignment_into_junctions(header, aln);
         }
         hts_itr_destroy(iter);

src/junctions/junctions_extractor.h

@@ -180,12 +180,21 @@ class JunctionsExtractor {
         string strand_tag_;
         //tag used in BAM to denote single cell barcode
         string barcode_tag_;
+        //filter reads containing any of these flags
+        uint16_t filter_flags_;
+        // filter reads not containing all of these flags
+        uint16_t require_flags_;
+        // filter reads below the minimum mapping quality
+        uint8_t min_map_qual_;
     public:
         //Default constructor
         JunctionsExtractor() {
             min_anchor_length_ = 8;
             min_intron_length_ = 70;
             max_intron_length_ = 500000;
+            filter_flags_ = 0;
+            require_flags_ = 0;
+            min_map_qual_ = 0;
             junctions_sorted_ = false;
             strandness_ = -1;
             strand_tag_ = "XS";
@@ -204,14 +213,20 @@ class JunctionsExtractor {
                 uint32_t min_anchor_length1, 
                 uint32_t min_intron_length1, 
                 uint32_t max_intron_length1, 
+                uint16_t filter_flags, 
+                uint16_t require_flags, 
+                uint8_t min_map_qual, 
                 string ref1) : 
             bam_(bam1), 
             region_(region1), 
             strandness_(strandness1), 
             strand_tag_(strand_tag1), 
             min_anchor_length_(min_anchor_length1), 
-            min_intron_length_(min_intron_length1), 
+            min_intron_length_(min_anchor_length1), 
             max_intron_length_(max_intron_length1), 
+            filter_flags_(filter_flags), 
+            require_flags_(require_flags), 
+            min_map_qual_(min_map_qual), 
             ref_(ref1) {
             junctions_sorted_ = false;
             output_file_ = "NA";
@@ -241,6 +256,8 @@ class JunctionsExtractor {
         void create_junctions_vector();
         //Pull out the cigar string from the read
         int parse_read(bam_hdr_t *header, bam1_t *aln);
+        //Returns whether alignment should be filtered from junction analysis
+        bool filter_alignment(bam_hdr_t *header, bam1_t *aln);
         //Parse junctions from the read and store in junction map
         int parse_cigar_into_junctions(string chr, int read_pos,
                                        uint32_t *cigar, int n_cigar);