115 lines
4.7 KiB
Diff
115 lines
4.7 KiB
Diff
diff -pruN orig/caps.cpp new/caps.cpp
|
|
--- orig/caps.cpp 2012-12-15 17:13:23.000000000 +0200
|
|
+++ new/caps.cpp 2020-09-09 23:07:46.080566000 +0300
|
|
@@ -14,7 +14,7 @@
|
|
#include <gsl/gsl_statistics.h>
|
|
#include<sys/time.h>
|
|
#include<iomanip>
|
|
-
|
|
+#include <bits/stdc++.h>
|
|
|
|
|
|
|
|
@@ -69,6 +69,8 @@
|
|
const gsl_rng_type * T;
|
|
gsl_rng *r;
|
|
|
|
+vector<double> totaltempnew;
|
|
+double alphathresh = 0;
|
|
int main(int argc, char *argv[]){
|
|
|
|
|
|
@@ -543,16 +545,27 @@ int main(int argc, char *argv[]){
|
|
|
|
|
|
print_splash(output);
|
|
+ OUTPUT << "\n\File1: " << files[i] << endl;
|
|
vec1.print_to_fasta(output.c_str());
|
|
+ OUTPUT << "\n\nFile2: " << files[j] << endl;
|
|
vec2.print_to_fasta(output.c_str());
|
|
int length1 = vec1.sequences[0].length();
|
|
int length2 = vec2.sequences[0].length();
|
|
|
|
+ OUTPUT << "\n\nLength1: " << length1 << endl;
|
|
+ OUTPUT << "Length2: " << length2 << endl;
|
|
|
|
|
|
if(tree_in ==0){
|
|
tree1 = create_input_tree(vec1.names, vec1.sequences);
|
|
tree2 = create_input_tree(vec2.names, vec2.sequences);
|
|
+
|
|
+ // Output the CAPS generated trees to the .out file of each pair
|
|
+ string temptre1 = TreeTemplateTools::treeToParenthesis(*tree1, true);
|
|
+ string temptre2 = TreeTemplateTools::treeToParenthesis(*tree2, true);
|
|
+ OUTPUT << "\n" << endl;
|
|
+ OUTPUT << "CAPS generated tree 1: " << temptre1 << endl;
|
|
+ OUTPUT << "CAPS generated tree 2: " << temptre2 << endl;
|
|
}/*else if(tree_in ==1 && variable==1){
|
|
|
|
std::auto_ptr<DistanceMatrix> DS;
|
|
@@ -666,6 +679,7 @@ int main(int argc, char *argv[]){
|
|
int value = floor(((totaltemp.size())*(1-(threshval))))+1;
|
|
|
|
threshold = totaltemp[value];
|
|
+ totaltempnew = totaltemp;
|
|
|
|
|
|
/*=======================================================*/
|
|
@@ -870,6 +884,30 @@ int Chi_squared (int num_pairs, int num_
|
|
|
|
} /* ----- end of function Chi_squared ----- */
|
|
|
|
+/*
|
|
+ * === FUNCTION ======================================================================
|
|
+ * Name: find_alpha
|
|
+ * Description: Find the index of an element in a vector totaltemp
|
|
+ * Help from: https://www.geeksforgeeks.org/how-to-find-index-of-a-given-element-in-a-vector-in-cpp/
|
|
+ * https://stackoverflow.com/questions/8647635/elegant-way-to-find-closest-value-in-a-vector-from-above
|
|
+ * Author: Petar Petrov, University of Turku (Finland); pebope@utu.fi
|
|
+ * =====================================================================================
|
|
+ */
|
|
+double getIndex(std::vector<double> const& v, double K)
|
|
+{
|
|
+ auto const it = std::lower_bound(v.begin(), v.end(), fabs(K));
|
|
+ //auto it = std::upper_bound(v.begin(), v.end(), fabs(K));
|
|
+
|
|
+ if (it != v.end()) {
|
|
+ int index = distance(v.begin(), it);
|
|
+ alphathresh = (((int)1+(double)v.size()-(int)index)/(double)v.size());
|
|
+ return alphathresh;
|
|
+ //cerr << index << "\t" << alphathresh << endl;
|
|
+ }
|
|
+ else {
|
|
+ cerr << "ELEMENT NOT FOUND!" << endl;
|
|
+ }
|
|
+}
|
|
|
|
|
|
|
|
@@ -890,9 +928,9 @@ int print_inter(vector<double>& Correl1,
|
|
output << endl << endl;
|
|
|
|
output << "Coevolving Pairs of amino acid sites\n";
|
|
- output << "=============================================================================\n";
|
|
- output << "Col1(real)\tCol2(real)\tDmean1\t\tDmean2\t\tCorrelation\tBootstrap value\n\n";
|
|
- output << "=============================================================================\n";
|
|
+ output << "================================================================================================================================\n";
|
|
+ output << "Col1(real)\tCol2(real)\tDmean1\t\tDmean2\t\tCorrelation\tBootstrap value\tP-value1\tP-value2\tMean P-value\tCorrelation1\tCorrelation2\n\n";
|
|
+ output << "================================================================================================================================\n";
|
|
|
|
//double mean = average_vec<double>(Correl);
|
|
//double SD = SD_vf(Correl, mean);
|
|
@@ -951,9 +989,11 @@ int print_inter(vector<double>& Correl1,
|
|
|
|
// }
|
|
|
|
+ double Alpha1 = getIndex(totaltempnew, Correl1[cor]);
|
|
+ double Alpha2 = getIndex(totaltempnew, Correl2[cor]);
|
|
//if(bootval>=bootcut && re1<=8 && re2<=8 ){
|
|
if(bootval>=bootcut){
|
|
- output << i+1 << "(" << i-gaps1+1 << ")\t\t" << j+1 << "(" << (j+1)-gaps2 << ")\t\t" << averDi << "\t\t" << averDj << "\t\t" << (Correl1[cor]+Correl2[cor])/2 << "\t" << bootval << endl;
|
|
+ output << i+1 << "(" << i-gaps1+1 << ")\t\t" << j+1 << "(" << (j+1)-gaps2 << ")\t\t" << averDi << "\t\t" << averDj << "\t" << (Correl1[cor]+Correl2[cor])/2 << "\t" << bootval << "\t" << Alpha1 << "\t" << Alpha2 << "\t" << (Alpha1+Alpha2)/2 << "\t" << Correl1[cor] << "\t" << Correl2[cor] << endl;
|
|
signif.push_back(((Correl1[cor]+Correl2[cor])/2));
|
|
++pairs;
|
|
vector<int> tem;
|