Merge branch 'master' into improved_compare_junctions

Author: kcotto

src/cis-splice-effects/cis_splice_effects_identifier.cc

@@ -40,6 +40,7 @@ void CisSpliceEffectsIdentifier::usage(ostream& out) {
     out << "\t\t" << "-j STR\tOutput file containing the aberrant junctions in BED12 format." << endl;
     out << "\t\t" << "-s INT\tStrand specificity of RNA library preparation \n"
         << "\t\t\t " << "(0 = unstranded, 1 = first-strand/RF, 2, = second-strand/FR). REQUIRED" << endl;
+    out << "\t\t" << "-t STR\tTag used in bam to label strand. [XS]" << endl;
     out << "\t\t" << "-a INT\tMinimum anchor length. Junctions which satisfy a minimum \n"
         << "\t\t\t " << "anchor length on both sides are reported. [8]" << endl;
     out << "\t\t" << "-m INT\tMinimum intron length. [70]" << endl;
@@ -103,7 +104,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:IShs:a:m:M:")) != -1) {
+    while((c = getopt(argc, argv, "o:w:v:j:e:Ei:ISht:s:a:m:M:")) != -1) {
         switch(c) {
             case 'o':
                 output_file_ = string(optarg);
@@ -138,6 +139,9 @@ void CisSpliceEffectsIdentifier::parse_options(int argc, char* argv[]) {
             case 's':
                 strandness_ = atoi(optarg);
                 break;
+            case 't':
+                strand_tag_ = string(optarg);
+                break;
             case 'a':
                 min_anchor_length_ = atoi(optarg);
                 break;
@@ -246,7 +250,7 @@ void CisSpliceEffectsIdentifier::identify() {
             if(write_annotated_variants_)
                 va.write_annotation_output(v1);
             //Extract junctions near this variant
-            JunctionsExtractor je1(bam_, variant_region, strandness_, min_anchor_length_, min_intron_length_, max_intron_length_);
+            JunctionsExtractor je1(bam_, variant_region, strandness_, strand_tag_, min_anchor_length_, min_intron_length_, max_intron_length_);
             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

@@ -80,6 +80,8 @@ class CisSpliceEffectsIdentifier {
         bool skip_single_exon_genes_;
         //strandness of data
         int strandness_;
+        //tag used in BAM to denote strand, default "XS"
+        string strand_tag_;
         //Minimum anchor length for junctions
         //Junctions need atleast this many bp overlap
         // on both ends.
@@ -101,6 +103,7 @@ class CisSpliceEffectsIdentifier {
                                        all_exonic_space_(false),
                                        skip_single_exon_genes_(true),
                                        strandness_(-1),
+                                       strand_tag_("XS"),
                                        min_anchor_length_(8),
                                        min_intron_length_(70),
                                        max_intron_length_(500000) {}

src/junctions/junctions_extractor.cc

@@ -41,7 +41,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:s:")) != -1) {
+    while((c = getopt(argc, argv, "ha:m:M:o:r:t:s:")) != -1) {
         switch(c) {
             case 'a':
                 min_anchor_length_ = atoi(optarg);
@@ -58,6 +58,9 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
             case 'r':
                 region_ = string(optarg);
                 break;
+            case 't':
+                strand_tag_ = string(optarg);
+                break;
             case 'h':
                 usage(help_ss);
                 throw common::cmdline_help_exception(help_ss.str());
@@ -104,6 +107,8 @@ int JunctionsExtractor::usage(ostream& out) {
         << "\t\t\t " << "in \"chr:start-end\" format. Entire BAM by default." << endl;
     out << "\t\t" << "-s INT\tStrand specificity of RNA library preparation \n"
         << "\t\t\t " << "(0 = unstranded, 1 = first-strand/RF, 2, = second-strand/FR). REQUIRED" << endl;
+    out << "\t\t" << "-t STR\tTag used in bam to label strand. [XS]" << endl;
+        
     out << endl;
     return 0;
 }
@@ -224,7 +229,7 @@ void JunctionsExtractor::print_all_junctions(ostream& out) {
 
 //Get the strand from the XS aux tag
 void JunctionsExtractor::set_junction_strand_XS(bam1_t *aln, Junction& j1) {
-    uint8_t *p = bam_aux_get(aln, "XS");
+    uint8_t *p = bam_aux_get(aln, strand_tag_.c_str());
     if(p != NULL) {
         char strand = bam_aux2A(p);
         strand ? j1.strand = string(1, strand) : j1.strand = string(1, '?');

src/junctions/junctions_extractor.h

@@ -155,6 +155,8 @@ class JunctionsExtractor {
         string region_;
         //strandness of data; 0 = unstranded, 1 = RF, 2 = FR
         int strandness_;
+        //tag used in BAM to denote strand, default "XS"
+        string strand_tag_;
     public:
         //Default constructor
         JunctionsExtractor() {
@@ -163,12 +165,13 @@ class JunctionsExtractor {
             max_intron_length_ = 500000;
             junctions_sorted_ = false;
             strandness_ = -1;
+            strand_tag_ = "XS";
             bam_ = "NA";
             output_file_ = "NA";
             region_ = ".";
         }
-        JunctionsExtractor(string bam1, string region1, int strandness1, uint32_t min_anchor_length1, uint32_t min_intron_length1, uint32_t max_intron_length1) : 
-            bam_(bam1), region_(region1), strandness_(strandness1), min_anchor_length_(min_anchor_length1), min_intron_length_(min_anchor_length1), max_intron_length_(max_intron_length1) {
+        JunctionsExtractor(string bam1, string region1, int strandness1, string strand_tag1, uint32_t min_anchor_length1, uint32_t min_intron_length1, uint32_t max_intron_length1) : 
+            bam_(bam1), region_(region1), strandness_(strandness1), strand_tag_(strand_tag1), min_anchor_length_(min_anchor_length1), min_intron_length_(min_anchor_length1), max_intron_length_(max_intron_length1) {
             junctions_sorted_ = false;
             output_file_ = "NA";
         }