switch RF and FR int values; RF is now 1 and FR is now 2 to match wit…

…h first-strand/second-strand descriptors

src/cis-splice-effects/cis_splice_effects_identifier.cc

@@ -38,10 +38,10 @@ void CisSpliceEffectsIdentifier::usage(ostream& out) {
     out << "Options:" << endl;
     out << "\t"   << "-o STR Output file containing the aberrant splice junctions with annotations. [STDOUT]" << endl;
     out << "\t\t" << "-v STR Output file containing variants annotated as splice relevant (VCF format)." << endl;
+    out << "\t\t" << "-j STR Output file containing the aberrant junctions in BED12 format." << endl;
     out << "\t\t" << "-w INT\tWindow size in b.p to identify splicing events in. "
         << "\t\t\t" << "The tool identifies events in variant.start +/- w basepairs."
         << "\t\t\t" << "Default behaviour is to look at the window between previous and next exons." << endl;
-    out << "\t\t" << "-j STR Output file containing the aberrant junctions in BED12 format." << endl;
     out << "\t\t" << "-e INT\tMaximum distance from the start/end of an exon "
                        "\t\t\tto annotate a variant as relevant to splicing, the variant "
                        "\t\t\tis in exonic space, i.e a coding variant. [3]" << endl;

src/cis-splice-effects/cis_splice_effects_identifier.h

@@ -92,7 +92,7 @@ class CisSpliceEffectsIdentifier {
                                        all_intronic_space_(false),
                                        all_exonic_space_(false),
                                        skip_single_exon_genes_(true),
-                                       strandness_(2) {}
+                                       strandness_(1) {}
         //Destructor
         ~CisSpliceEffectsIdentifier() {
             if(ofs_.is_open()) {

src/junctions/junctions_extractor.cc

@@ -97,6 +97,7 @@ int JunctionsExtractor::usage(ostream& out) {
     out << "\t\t" << "-o FILE\tThe file to write output to. [STDOUT]" << endl;
     out << "\t\t" << "-r STR\tThe region to identify junctions "
                      "in \"chr:start-end\" format. Entire BAM by default." << endl;
+    out << "\t\t" << "-s INT\tStrand specificity of RNA library preparation. (0 = unstranded, 1 = first-strand/RF, 2, = second-strand/FR) [1]" << endl;
     out << endl;
     return 0;
 }
@@ -227,10 +228,10 @@ void JunctionsExtractor::set_junction_strand_flag(bam1_t *aln, Junction& j1) {
     int mate_reversed = (flag >> 5) % 2;
     int first_in_pair = (flag >> 6) % 2;
     int second_in_pair = (flag >> 7) % 2;
-    // strandness_ is 0 for unstranded, 1 for FR, and 2 for RF
+    // strandness_ is 0 for unstranded, 1 for RF, and 2 for FR
     int bool_strandness = strandness_ - 1;
-    int first_strand = bool_strandness ^ first_in_pair ^ reversed;
-    int second_strand = bool_strandness ^ second_in_pair ^ mate_reversed;
+    int first_strand = !bool_strandness ^ first_in_pair ^ reversed;
+    int second_strand = !bool_strandness ^ second_in_pair ^ mate_reversed;
     char strand;
     if (first_strand){
         strand = '+';
@@ -251,22 +252,13 @@ void JunctionsExtractor::set_junction_strand_flag(bam1_t *aln, Junction& j1) {
 //Get the strand
 void JunctionsExtractor::set_junction_strand(bam1_t *aln, Junction& j1) {
     // if unstranded data
-    //cerr << "strandness is " << strandness_ << endl;
     if (strandness_ > 0){
         return set_junction_strand_flag(aln, j1);
     } else {
         return set_junction_strand_XS(aln, j1);
     }
 }
 
-// // test
-// void JunctionsExtractor::set_junction_strand(bam1_t *aln, Junction& j1) {
-//     // if unstranded data
-//     cerr << "strandness is " << strandness_ << endl;
-//     set_junction_strand_flag(aln, j1);
-//     return set_junction_strand_XS(aln, j1);
-// }
-
 //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;
@@ -278,15 +270,6 @@ int JunctionsExtractor::parse_alignment_into_junctions(bam_hdr_t *header, bam1_t
     string chr(header->target_name[chr_id]);
     uint32_t *cigar = bam_get_cigar(aln);
 
-    /*
-    //Skip duplicates
-    int flag = aln->core.flag;
-    if(flag & 1024) {
-        cerr << "Skipping read_pos " << read_pos << " flag " << flag << endl;
-        return 0;
-    }
-    */
-
     Junction j1;
     j1.chrom = chr;
     j1.start = read_pos; //maintain start pos of junction

src/junctions/junctions_extractor.h

@@ -145,7 +145,7 @@ class JunctionsExtractor {
         string output_file_;
         //Region to identify junctions, in "chr:start-end" format
         string region_;
-        //strandness off data; 0 = unstranded, 1 = FR, 2 = RF
+        //strandness of data; 0 = unstranded, 1 = RF, 2 = FR
         int strandness_;
     public:
         //Default constructor
@@ -155,14 +155,13 @@ class JunctionsExtractor {
             min_intron_length_ = 70;
             max_intron_length_ = 500000;
             junctions_sorted_ = false;
-            strandness_ = 2;
+            strandness_ = 1;
             bam_ = "NA";
             output_file_ = "NA";
             region_ = ".";
         }
         JunctionsExtractor(string bam1, string region1, int strandness1) : bam_(bam1), region_(region1), strandness_(strandness1) {
             //cerr << "param constructor called" << endl;
-            //cerr << "strandness1 is " << strandness1 << endl;
             min_anchor_length_ = 8;
             min_intron_length_ = 70;
             max_intron_length_ = 500000;