Page MenuHomec4science
Files F87431085

main_test.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Oct 12, 15:56

main_test.cpp
View Options

#include <iostream>
#include <string>
#include <Matrix2D.hpp>
#include <Matrix4D.hpp>
#include <CorrelationMatrixCreator.hpp>
#
using namespace std ;
template<class T>
std::ostream& operator << (std::ostream& stream, const std::vector<T>& v)
{ for(const auto x : v)
{ stream << x << " " ; }
return stream ;
}
int main()
{ // path
std::string bed_file = "/local/groux/scATAC-seq/data/10xgenomics_PBMC_5k/atac_v1_pbmc_5k_peaks_rmsk_sampled.bed" ;
std::string bam_file = "/local/groux/scATAC-seq/data/10xgenomics_PBMC_5k/atac_v1_pbmc_5k_possorted_filtered_30-84bp.bam" ;
std::string bai_file = "/local/groux/scATAC-seq/data/10xgenomics_PBMC_5k/atac_v1_pbmc_5k_possorted_filtered_30-84bp.bam.bai" ;
std::string prob_file = "/local/groux/scATAC-seq/results/10xgenomics_PBMC_5k_peaks_classification_6/peaks_rmsk_sampled_sequences_1kb_23class_prob.mat4d" ;
// posterior prob
std::cerr << "loading posterior prob" << std::endl ;
Matrix4D<double> prob(prob_file) ;
size_t n_row = prob.get_dim()[0] ;
size_t n_class = prob.get_dim()[1] ;
size_t n_shift = prob.get_dim()[2] ;
size_t n_flip = prob.get_dim()[3] ;
bool flip = n_flip == 2 ;
std::cerr << "posterior prob " << prob.get_dim() << std::endl ;
// class prob
std::cerr << "computing class probabilities" << std::endl ;
std::vector<double> prob_colsum(n_class, 0.) ;
double tot = 0. ;
for(size_t i=0; i<n_row; i++)
{ for(size_t j=0; j<n_class; j++)
{ for(size_t k=0; k<n_shift; k++)
{ for(size_t l=0; l<n_flip; l++)
{ double p = prob(i,j,k,l) ;
prob_colsum[j] += p ;
tot += p ;
}
}
}
}
for(auto& p : prob_colsum)
{ p /= tot ; }
// these are from3 and to3
int from1 = -1000 ;
int to1 = 1000 ;
// class of interest
size_t class_k = 17 ;
class_k -= 1 ;
// get a wider matrix to realign it
std::cerr << "creating data2 matrix" << std::endl ;
// extend limits
int ext = n_shift - 1 ;
int ext_r = ext / 2 ;
int ext_l = ext - ext_r ;
int from2 = from1 - ext_l ;
int to2 = to1 + ext_r ;
CorrelationMatrixCreator mc(bed_file,
bam_file,
bai_file,
from2,
to2,
1,
CorrelationMatrixCreator::READ_ATAC) ;
Matrix2D<int> data2 = mc.create_matrix() ;
size_t n_col2 = data2.get_ncol() ;
std::cerr << "data2 matrix " << data2.get_dim() << std::endl ;
// realign matrix
std::cerr << "computing data3 matrix" << std::endl ;
size_t n_col3 = to1 - from1 + 1 ;
Matrix2D<double> data3(n_row,
n_col3,
0.) ;
std::cerr << "data3 matrix " << data3.get_dim() << std::endl ;
for(size_t i=0; i<n_row; i++)
{ std::cerr << i << std::endl ;
for(size_t s=0; s<n_shift; s++)
{ // --------------- forward ---------------
int from_dat2_fw = s ;
int to_dat2_fw = from_dat2_fw + n_col3 - 1 ;
for(int j_dat2_fw=from_dat2_fw, j_dat3_fw=0;
j_dat2_fw<=to_dat2_fw;
j_dat2_fw++, j_dat3_fw++)
{ data3(i,j_dat3_fw) +=
(prob(i,class_k,s,0) *
data2(i,j_dat2_fw)) /
prob_colsum[class_k] ;
}
// --------------- reverse ---------------
if(flip)
{ int from_dat2_rev = n_col2 - 1 - s ;
int to_dat2_rev = from_dat2_rev - (n_col3 - 1) ;
for(int j_dat2_rev=from_dat2_rev, j_dat3_fw=0;
j_dat2_rev >= to_dat2_rev;
j_dat2_rev--, j_dat3_fw++)
{ data3(i,j_dat3_fw) +=
(prob(i,class_k,s,1) *
data2(i,j_dat2_rev)) /
prob_colsum[class_k] ;
}
}
}
}
// clean memory
prob = Matrix4D<double>() ;
// convert to integer
Matrix2D<int> data4(data3.get_nrow(),
data3.get_ncol()) ;
for(size_t i=0; i< data4.get_nrow(); i++)
{ for(size_t j=0; j<data4.get_ncol(); j++)
{ data4(i,j) = (int)data3(i,j) ; }
}
std::cout << data4 << std::endl ;
return EXIT_SUCCESS ;
}

Event Timeline

c4science · Help