incorporate changes from https://www.nature.com/articles/srep39734 to…

config/parameters/aug_cmdln_parameters.json

@@ -717,6 +717,11 @@
         "development": true,
         "description": ""
     },
+    {
+        "name": "/IntronModel/allow_dss_consensus_ga",
+        "development": true,
+        "description": ""
+    },
     {
         "name": "/IntronModel/ass_motif_memory",
         "development": true,
@@ -769,7 +774,17 @@
         "description": ""
     },
     {
-        "name": "/IntronModel/non_gt_dss_prob",
+        "name": "/IntronModel/gt_dss_prob",
+        "development": true,
+        "description": ""
+    },
+    {
+        "name": "/IntronModel/gc_dss_prob",
+        "development": true,
+        "description": ""
+    },
+    {
+        "name": "/IntronModel/ga_dss_prob",
         "development": true,
         "description": ""
     },

docs/tutorial/data/ppx.cfg

@@ -9,5 +9,4 @@ gff3		true
 print_blocks	true
 
 /IntronModel/allow_dss_consensus_gc	true
-/IntronModel/non_gt_dss_prob	0.01
-
+/IntronModel/gc_dss_prob	0.01

docs/tutorial/ppx.html

@@ -329,7 +329,7 @@ <h2>4. RUN AUGUSTUS-PPX</h2>
 print_blocks    true
 
 /IntronModel/allow_dss_consensus_gc     true
-/IntronModel/non_gt_dss_prob    0.01
+/IntronModel/gc_dss_prob    0.01
 </pre>
 
 Here, we do not use sampling (sampling is not supported in the protein extension), 

docs/tutorial2015/ppx.html

@@ -329,7 +329,7 @@ <h2>4. RUN AUGUSTUS-PPX</h2>
 print_blocks    true
 
 /IntronModel/allow_dss_consensus_gc     true
-/IntronModel/non_gt_dss_prob    0.01
+/IntronModel/gc_dss_prob    0.01
 </pre>
 
 Here, we do not use sampling (sampling is not supported in the protein extension), 

include/geneticcode.hh

@@ -32,25 +32,38 @@ using namespace std;
 	return strncmp(dna, PATTERN, COUNT) == 0;	\
     }
 // 'gt'
-#define DSS_SEQUENCE "gt"
-#define RDSS_SEQUENCE "ac"
-DECLARE_ON(onDSS,     DSS_SEQUENCE, 2)
-DECLARE_ON(onRDSS,    RDSS_SEQUENCE, 2)
+#define GT_DSS_SEQUENCE "gt"
+#define GT_RDSS_SEQUENCE "ac"
+DECLARE_ON(onGtDSS,     GT_DSS_SEQUENCE, 2)
+DECLARE_ON(onGtRDSS,    GT_RDSS_SEQUENCE, 2)
 
 // 'gc'
-#define ALT_DSS_SEQUENCE "gc"
-#define ALT_RDSS_SEQUENCE "gc"
-DECLARE_ON(onAltDSS,  ALT_DSS_SEQUENCE, 2)
-DECLARE_ON(onAltRDSS, ALT_RDSS_SEQUENCE, 2)
+#define GC_DSS_SEQUENCE "gc"
+#define GC_RDSS_SEQUENCE "gc"
+DECLARE_ON(onGcDSS,  GC_DSS_SEQUENCE, 2)
+DECLARE_ON(onGcRDSS, GC_RDSS_SEQUENCE, 2)
 
-// 'gt' or 'gc'
+// added for GA-AG splice site [https://www.nature.com/articles/srep39734, YONG LI]
+//
+// 'ga'
+#define GA_DSS_SEQUENCE "ga"
+#define GA_RDSS_SEQUENCE "tc"
+DECLARE_ON(onGaDSS,  GA_DSS_SEQUENCE, 2)
+DECLARE_ON(onGaRDSS, GA_RDSS_SEQUENCE, 2)
+
+
+// changed for GA-AG splice site [https://www.nature.com/articles/srep39734, YONG LI]
+//
+// 'gt' or 'gc' or 'ga'
 inline bool onGenDSS(const char* dna) {
     return 
-	onDSS(dna) || (Constant::dss_gc_allowed && onAltDSS(dna));
+        onGtDSS(dna) || (Constant::dss_gc_allowed && onGcDSS(dna)) || (Constant::dss_ga_allowed && onGaDSS(dna));
+
 }
 inline bool onGenRDSS(const char* dna) {
     return 
-	onRDSS(dna) || (Constant::dss_gc_allowed && onAltRDSS(dna));
+        onGtRDSS(dna) || (Constant::dss_gc_allowed && onGcRDSS(dna)) || (Constant::dss_ga_allowed && onGaRDSS(dna));
+
 }
 
 // 'ag'

include/intronmodel.hh

@@ -123,7 +123,9 @@ private:
     static Double          dssneighborfactor; // taken from the prob of the neighbor patterns
     static Integer         ass_motif_memory;  // order of the markov chain in the ass upstream motif
     static Integer         ass_motif_radius;  // radius of the pooling window in the ass upstream motif
-    static double          non_gt_dss_prob;
+    static double          gt_dss_prob;       // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+    static double          gc_dss_prob;       // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+    static double          ga_dss_prob;       // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
     static double          non_ag_ass_prob;
     static SnippetProbs    *snippetProbs, *rSnippetProbs;  
     static bool            initAlgorithmsCalled, haveSnippetProbs;

include/types.hh

@@ -359,6 +359,8 @@ public:
     static bool CRFtrainCDS;
     static bool CRFtrainUTR;
     static bool dss_gc_allowed;
+    // added for GA-AG splice site [https://www.nature.com/articles/srep39734 YONG LI]
+    static bool dss_ga_allowed;
     static Boolean tieIgenicIntron; // whether to tie igenic model parameters to intron model parameters, i.e. use just one content model, that of the intron
     static Boolean exoncands;
     static Boolean proteinOutput;

src/etraining.cc

@@ -73,6 +73,7 @@ int main( int argc, char* argv[] ){
 	}
 
 	Constant::dss_gc_allowed = true; // default value for training!
+        Constant::dss_ga_allowed = true; // added for GA-AG splice site [https://www.nature.com/articles/srep39734, YONG LI]
 	Constant::init();
 	BaseCount::init();  // replaces train.init()
 	GeneticCode::init();

src/exoncand.cc

@@ -192,7 +192,7 @@ void findExonCands(map<int_fast64_t, ExonCandidate*> &ecs, map<int_fast64_t, Exo
             }
         }
         // positions of all reverse DSS "ac"
-        if (withSplicing && onRDSS(dna+i) && (i + Constant::dss_whole_size() - Constant::dss_end < n)) {
+        if (withSplicing && onGtRDSS(dna+i) && (i + Constant::dss_whole_size() - Constant::dss_end < n)) {
             p = IntronModel::dSSProb(i - Constant::dss_end, false);
             if (p >= dssminprob) {
                 ssWithScore.first = i;
@@ -249,7 +249,7 @@ void findExonCands(map<int_fast64_t, ExonCandidate*> &ecs, map<int_fast64_t, Exo
         }
 
         // computing initial exons on the forward strand with at least startcodon plus base
-        if (withSplicing && onDSS(dna + i) && (i + Constant::dss_whole_size() - Constant::dss_start  < n)) {
+        if (withSplicing && onGtDSS(dna + i) && (i + Constant::dss_whole_size() - Constant::dss_start  < n)) {
             p = IntronModel::dSSProb(i - Constant::dss_start,true);
             for (frame=0; frame<=2; frame++) {
                 ritStart=ritStart_cur;
@@ -584,17 +584,19 @@ bool ExonCandidate::correctType(const char* dna, int dnalen){
 
     bool correctType = false;
 
+    // add support for GA/GC DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+    //
     switch(type){
     case singleGene:
         if(onStart(dna+begin) && GeneticCode::isStopcodon(dna+end-2))
             correctType=true;
         break;
     case initial_0: case initial_1: case initial_2:
-        if(onStart(dna+begin) && (rightTruncated || onDSS(dna+end+1) || onAltDSS(dna+end+1)))
+        if(onStart(dna+begin) && (rightTruncated || onGtDSS(dna+end+1) || onGcDSS(dna+end+1) || onGaDSS(dna+end+1)))
             correctType=true;
         break;
     case internal_0: case internal_1: case internal_2:
-        if((leftTruncated || onASS(dna+begin-2)) && (rightTruncated || onDSS(dna+end+1) || onAltDSS(dna+end+1)))
+        if((leftTruncated || onASS(dna+begin-2)) && (rightTruncated || onGtDSS(dna+end+1) || onGcDSS(dna+end+1) || onGaDSS(dna+end+1)))
             correctType=true;
         break;
     case terminal_exon:
@@ -606,11 +608,11 @@ bool ExonCandidate::correctType(const char* dna, int dnalen){
             correctType=true;
         break;
     case rinitial_exon:
-        if(onRStart(dna+end-2) && (leftTruncated || onRDSS(dna+begin-2) || onAltRDSS(dna+begin-2)))
+        if(onRStart(dna+end-2) && (leftTruncated || onGtRDSS(dna+begin-2) || onGcRDSS(dna+begin-2) || onGaRDSS(dna+begin-2)))
             correctType=true;
         break;
     case rinternal_0: case rinternal_1: case rinternal_2:
-        if((rightTruncated || onRASS(dna+end+1)) && (leftTruncated || onRDSS(dna+begin-2) || onAltRDSS(dna+begin-2)))
+        if((rightTruncated || onRASS(dna+end+1)) && (leftTruncated || onGtRDSS(dna+begin-2) || onGcRDSS(dna+begin-2) || onGaRDSS(dna+begin-2)))
             correctType=true;
         break;
     case  rterminal_0: case rterminal_1: case rterminal_2:

src/extrinsicinfo.cc

@@ -21,7 +21,8 @@
 #include <climits>
 #include <iterator> // for ostream_iterator
 
-Bitmask     SequenceFeatureCollection::validDSS = A_SET_FLAG(Seq2Int(2)("gt")) | A_SET_FLAG(Seq2Int(2)("gc"));
+// Add support for GA site - [https://www.nature.com/articles/srep39734, BC]
+Bitmask     SequenceFeatureCollection::validDSS = A_SET_FLAG(Seq2Int(2)("gt")) | A_SET_FLAG(Seq2Int(2)("gc")) | A_SET_FLAG(Seq2Int(2)("ga"));
 Bitmask     SequenceFeatureCollection::validASS = A_SET_FLAG(Seq2Int(2)("ag"));
 set<string> SequenceFeatureCollection::validHintedSites;
 
@@ -2595,7 +2596,7 @@ void FeatureCollection::printAccuracyForSequenceSet(const AnnoSequence* annoseqs
 
 	    // 'donor splice site' feature
 	    num_correct_pos[dssF]   += dssAS.size();
-	    num_incorrect_pos[dssF] += 2*(howOftenOccursIt(curannoseq->sequence, DSS_SEQUENCE)) - dssAS.size();
+	    num_incorrect_pos[dssF] += 2*(howOftenOccursIt(curannoseq->sequence, GT_DSS_SEQUENCE)) - dssAS.size();
 
 	    flist = sfc.getFeatureList(dssF);
 	    flist.sort();

src/intronmodel.cc

@@ -41,7 +41,9 @@ vector<Double>  IntronModel::lenDist;
 Integer         IntronModel::gesbasen = 0;
 Integer         IntronModel::ass_motif_memory = 3;
 Integer         IntronModel::ass_motif_radius = 3;
-double          IntronModel::non_gt_dss_prob = 0.001;
+double          IntronModel::gt_dss_prob = 1;      // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+double          IntronModel::gc_dss_prob = 0.001;  // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+double          IntronModel::ga_dss_prob = 0;      // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
 double          IntronModel::non_ag_ass_prob = 0.001;
 double          IntronModel::geoProb = 0.9997447;
 
@@ -146,7 +148,10 @@ void IntronModel::init() {
     }catch( ProjectError &e) { 
 	cerr << e.getMessage();
     }
-    Properties::assignProperty("/IntronModel/non_gt_dss_prob", non_gt_dss_prob);
+    // add support for GT/GC/GA DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+    Properties::assignProperty("/IntronModel/gt_dss_prob", gt_dss_prob);
+    Properties::assignProperty("/IntronModel/gc_dss_prob", gc_dss_prob);
+    Properties::assignProperty("/IntronModel/ga_dss_prob", ga_dss_prob);
     Properties::assignProperty("/IntronModel/non_ag_ass_prob", non_ag_ass_prob);
 
     ass_upwindow_size = Constant::ass_upwindow_size;
@@ -1211,26 +1216,35 @@ Double IntronModel::dSSProb(int base, bool forwardStrand){
     if (!forwardStrand && (got = memoR.find(base), got != memoR.end()))
     	return got->second;
 
-    bool nonGT;
+    // add support for GA/GC DSS probabilities [https://www.sciencedirect.com/science/article/pii/S0960982213006878, BC]
+    bool isGtDss;
+    bool isGcDss;
+    bool isGaDss;
     if (forwardStrand) { // forward strand
 	int dsspos = base + Constant::dss_start;
 	if (!isPossibleDSS(dsspos))
 	    return 0;
-	nonGT = !onDSS(sequence + dsspos);
+	isGtDss = onGtDSS(sequence + dsspos);
+	isGcDss = onGcDSS(sequence + dsspos);
+	isGaDss = onGaDSS(sequence + dsspos);
 	strncpy(astr, sequence + base, Constant::dss_start);
 	strncpy(astr + Constant::dss_start, sequence + dsspos + DSS_MIDDLE, Constant::dss_end);
     } else { // reverse complement
 	int dsspos = base + Constant::dss_end;
 	if (!isPossibleRDSS(dsspos + 1))
 	    return 0;
-	nonGT = !onRDSS(sequence + dsspos);
+	isGtDss = onGtRDSS(sequence + dsspos);
+	isGcDss = onGcRDSS(sequence + dsspos);
+	isGaDss = onGaRDSS(sequence + dsspos);
 	putReverseComplement(astr, sequence + dsspos + DSS_MIDDLE, Constant::dss_start);
 	putReverseComplement(astr + Constant::dss_start, sequence + base, Constant::dss_end);
     }
     astr[Constant::dss_size()] = '\0';
     try {
 	Double dssprob = dssprobs[ s2i(astr) ];
-	if (nonGT) dssprob *= non_gt_dss_prob;
+	if (isGtDss) dssprob *= gt_dss_prob;
+	if (isGcDss) dssprob *= gc_dss_prob;
+	if (isGaDss) dssprob *= ga_dss_prob;
 	if (dssBinProbs.nbins >= 1) {
 	    int idx = dssBinProbs.getIndex(dssprob);
 	    if (inCRFTraining && (countEnd < 0 || (base >= countStart && base <= countEnd)))

src/introntrain.cc

@@ -270,6 +270,7 @@ void IntronModel::processDSS( const char* dna, int pos){
     if (!onGenDSS(dna+pos+1))
 	throw IntronModelError("DSS error! Expected 'gt'" +
 			       string(Constant::dss_gc_allowed ?  " or 'gc'" : "") +
+                               string(Constant::dss_ga_allowed ?  " or 'ga'" : "") +    // added for GA-AG splice site [https://www.nature.com/articles/srep39734 YONG LI]
 			       ", but found '" + string(dna, pos+1, 2) + "' at position " +
 			       itoa(pos+1) + ".");
     else {

src/types.cc

@@ -66,6 +66,8 @@ bool Constant::CRFtrainSS = true;
 bool Constant::CRFtrainTIS = true;
 bool Constant::CRFtrainUTR = false;
 bool Constant::dss_gc_allowed = false;
+// added for GA-AG splice site [https://www.nature.com/articles/srep39734, YONG LI]
+bool Constant::dss_ga_allowed = false;
 Boolean Constant::tieIgenicIntron = true;
 Boolean Constant::proteinOutput = true;
 Boolean Constant::codSeqOutput = false;
@@ -425,6 +427,8 @@ void Constant::init(){
     Properties::assignProperty(UTR_KEY, utr_option_on);
     Properties::assignProperty(NONCODING_KEY, nc_option_on);
     Properties::assignProperty("/IntronModel/allow_dss_consensus_gc", dss_gc_allowed);
+    // added for GA-AG splice site [https://www.nature.com/articles/srep39734, YONG LI]
+    Properties::assignProperty("/IntronModel/allow_dss_consensus_ga", dss_ga_allowed);
     Properties::assignProperty("tieIgenicIntron", tieIgenicIntron);
     Properties::assignProperty("protein", proteinOutput);
     Properties::assignProperty("codingseq", codSeqOutput);