From 06305387fb1ee0effe999b83262311d0c757d08b Mon Sep 17 00:00:00 2001
From: Incardona Pietro <incardon@mpi-cbg.de>
Date: Sat, 20 Apr 2024 19:17:53 +0200
Subject: [PATCH] Latest changes to CNVcaller
---
example/Sequencing/CNVcaller/functions.hpp | 105 +++++------
example/Sequencing/CNVcaller/testing.cpp | 192 +++++++++++++++++++--
2 files changed, 235 insertions(+), 62 deletions(-)
diff --git a/example/Sequencing/CNVcaller/functions.hpp b/example/Sequencing/CNVcaller/functions.hpp
index c917f3a73..3cabc347c 100644
--- a/example/Sequencing/CNVcaller/functions.hpp
+++ b/example/Sequencing/CNVcaller/functions.hpp
@@ -2374,8 +2374,8 @@ void collect_reads_for_call_evidence_forward(std::string & bam_file,
int k = 0;
int s = 0;
-// for (int k = -padding ; k <= padding ; k++) {
-// for (int s = -padding ; s <= padding ; s++) {
+ for (int k = -padding ; k <= padding ; k++) {
+ for (int s = -padding ; s <= padding ; s++) {
int base_jump = calls.get(cid).start+k - start_call;
@@ -2449,8 +2449,8 @@ void collect_reads_for_call_evidence_forward(std::string & bam_file,
goto found_read;
}
- // }
- //}
+ }
+ }
found_read:;
}
@@ -2851,7 +2851,7 @@ size_t add_32bases_pattern_with_variability(int chrom,
{
if ( !bcf_sr_has_line(sr,i) ) continue;
bcf1_t *rec = bcf_sr_get_line(sr, i);
- if (rec->pos >= start && rec->pos <= start+32) {
+ if (rec->pos >= start && rec->pos < start+32) {
for (int j=1; j<rec->n_allele; j++)
{
@@ -2862,7 +2862,7 @@ size_t add_32bases_pattern_with_variability(int chrom,
variants.last().alt = std::string(rec->d.allele[j]);
int n_float = 1;
- void ** ptr;
+ void * ptr[1];
ptr[0] = &variants.last().freq;
// Get the frequency from AF in info of rec
@@ -2938,7 +2938,7 @@ size_t add_32bases_pattern_with_variability(int chrom,
return bases_pattern;
}
-size_t add_32bases_pattern_with_variability_range_shift(int chrom,
+void add_32bases_pattern_with_variability_range_shift(int chrom,
int start,
int stop,
int shift,
@@ -2949,8 +2949,8 @@ size_t add_32bases_pattern_with_variability_range_shift(int chrom,
int padding_left,
int padding_right) {
- for (int i = start ; i < stop ; i+= shift) {
- add_32bases_pattern_with_variability(chrom,i,chr_sequences,map,database,call_id,padding_left+shift,padding_right);
+ for (int i = start ; i <= stop ; i+= shift) {
+ add_32bases_pattern_with_variability(chrom,i,chr_sequences,map,database,call_id,padding_left+(i-start),padding_right);
}
}
@@ -3213,8 +3213,8 @@ Eigen::ArrayXd Evidence_matrix_log(openfpm::vector<openfpm::vector<char>> & refe
int i = 0;
int j = 0;
-// for (int i = -padding ; i <= padding ; i++) {
-// for (int j = -padding ; j <= padding ; j++) {
+ for (int i = -padding ; i <= padding ; i++) {
+ for (int j = -padding ; j <= padding ; j++) {
int id = (i+padding)*(2*padding+1) + (j+padding);
int pos_event;
@@ -3228,72 +3228,77 @@ Eigen::ArrayXd Evidence_matrix_log(openfpm::vector<openfpm::vector<char>> & refe
}
get_sequence(reference,chr,pos_event,pos_event+seq.size()+(call.stop+j-call.start-i),call.start+i,call.stop+j,seq1);
- std::cout << "--------------- " << i << "," << j << " --------------------------------- event: " << id << std::endl;
events(id) = calculate_event_probability(seq2,seq1,seq3,call.stop-call.start,chr,pos,3e-3,chr_ins,reference,bp_insert_prob);
- std::cout << "EVENT ID: " << id << std::endl;
// Print
- for (int j = 0 ; j < seq.size() ; j++) {
- if (j == distance_from_seed) {
- std::cout << "*" << seq2[j];
- } else {
- std::cout << seq2[j];
- }
- }
- std::cout << std::endl;
+ if (events(id) > -60.0) {
+ std::cout << "--------------- " << i << "," << j << " --------------------------------- event: " << id << std::endl;
+ std::cout << "EVENT ID: " << id << std::endl;
- for (int j = 0 ; j < seq.size() ; j++) {
- if (j == distance_from_seed) {
- std::cout << "*" << seq1[j];
- } else {
- std::cout << seq1[j];
+ for (int j = 0 ; j < seq.size() ; j++) {
+ if (j == distance_from_seed) {
+ std::cout << "*" << seq2[j];
+ } else {
+ std::cout << seq2[j];
+ }
}
- }
- std::cout << std::endl;
+ std::cout << std::endl;
- std::cout << " Probability: " << events(id) << std::endl;
+ for (int j = 0 ; j < seq.size() ; j++) {
+ if (j == distance_from_seed) {
+ std::cout << "*" << seq1[j];
+ } else {
+ std::cout << seq1[j];
+ }
+ }
+ std::cout << std::endl;
+ std::cout << " Probability: " << events(id) << std::endl;
+ }
// Create a file
-// }
-// }
+ }
+ }
// Event 3 unaligned
// for all unaligned points (max 100)
for (int i = 0 ; i < un_seed.size() ; i++) {
get_sequence(reference,un_seed.get<0>(i),un_seed.get<1>(i),un_seed.get<1>(i)+32,-1,-1,seq1);
- std::cout << "UNALIGNED SEED: " << std::string(seq1) << std::endl;
get_sequence(reference,un_seed.get<0>(i),un_seed.get<1>(i)-distance_from_seed,un_seed.get<1>(i)-distance_from_seed+seq.size(),-1,-1,seq1);
int id = (2*padding+1)*(2*padding+1) + i;
- std::cout << "------------------- un_seed: " << i << " ----------------------------- event: " << id << std::endl;
events(id) = calculate_event_probability(seq2,seq1,seq3,-1,chr,un_seed.get<1>(i)-distance_from_seed+1,3e-3,chr_ins,reference,bp_insert_prob);
- std::cout << id << " SEED position: " << un_seed.get<1>(i) << " SEED ID: " << un_seed.get<2>(i) << std::endl;
+ if (events(id) > -60.0) {
+ std::cout << "------------------- un_seed: " << i << " ----------------------------- event: " << id << std::endl;
+ std::cout << id << " SEED position: " << un_seed.get<1>(i) << " SEED ID: " << un_seed.get<2>(i) << std::endl;
+
+ // Print
- // Print
+ std::cout << "UNALIGNED SEED: " << std::string(seq1) << std::endl;
- for (int j = 0 ; j < seq.size() ; j++) {
- if (j == distance_from_seed) {
- std::cout << "*" << seq2[j];
- } else {
- std::cout << seq2[j];
+ for (int j = 0 ; j < seq.size() ; j++) {
+ if (j == distance_from_seed) {
+ std::cout << "*" << seq2[j];
+ } else {
+ std::cout << seq2[j];
+ }
}
- }
- std::cout << std::endl;
+ std::cout << std::endl;
- for (int j = 0 ; j < seq.size() ; j++) {
- if (j == distance_from_seed) {
- std::cout << "*" << seq1[j];
- } else {
- std::cout << seq1[j];
+ for (int j = 0 ; j < seq.size() ; j++) {
+ if (j == distance_from_seed) {
+ std::cout << "*" << seq1[j];
+ } else {
+ std::cout << seq1[j];
+ }
}
- }
- std::cout << std::endl;
+ std::cout << std::endl;
- std::cout << " Probability: " << events(id) << std::endl;
+ std::cout << " Probability: " << events(id) << std::endl;
+ }
}
return events;
diff --git a/example/Sequencing/CNVcaller/testing.cpp b/example/Sequencing/CNVcaller/testing.cpp
index 4823a4bd6..e9e196d80 100644
--- a/example/Sequencing/CNVcaller/testing.cpp
+++ b/example/Sequencing/CNVcaller/testing.cpp
@@ -108,11 +108,14 @@ BOOST_AUTO_TEST_CASE( CNV_caller_test_add_32bases_pattern_with_variability_range
reference.load("reference_vector");
- // chromosome 0
- add_32bases_pattern_with_variability_range_shift(2,195452030,19545262,8,reference,map,"test/gnomad_variants.vcf.gz",0,-1,NO_PADDING);
+ // chromosome 3 // 195452062
+ add_32bases_pattern_with_variability_range_shift(2,195452030,195452062,8,reference,map,"test/gnomad_variants.vcf.gz",0,0,NO_PADDING);
- BOOST_REQUIRE_EQUAL(map.size(),8);
+ BOOST_REQUIRE_EQUAL(map.size(),8+8+8+8+4);
+ // From 30
+
+ {
// check that these seed are presents
std::string s1("TCATCACCCCGTCACGGGCCTCAGAGAGCAGC"); // reference
std::string s2("TCATCACCCCCTCACGGGCCTCAGAGAGCAGC"); // variants
@@ -133,6 +136,47 @@ BOOST_AUTO_TEST_CASE( CNV_caller_test_add_32bases_pattern_with_variability_range
size_t s7_s = CNVcaller::convert_to_bit_pattern(s7.c_str());
size_t s8_s = CNVcaller::convert_to_bit_pattern(s8.c_str());
+ BOOST_REQUIRE(map.find(s1_s) != map.end());
+ BOOST_REQUIRE(map.find(s2_s) != map.end());
+ BOOST_REQUIRE(map.find(s3_s) != map.end());
+ BOOST_REQUIRE(map.find(s4_s) != map.end());
+ BOOST_REQUIRE(map.find(s5_s) != map.end());
+ BOOST_REQUIRE(map.find(s6_s) != map.end());
+ BOOST_REQUIRE(map.find(s7_s) != map.end());
+ BOOST_REQUIRE(map.find(s8_s) != map.end());
+
+ BOOST_REQUIRE(map.find(s1_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s2_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s3_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s4_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s5_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s6_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s7_s)->second.get(0).padding_left == 0);
+ BOOST_REQUIRE(map.find(s8_s)->second.get(0).padding_left == 0);
+ }
+
+ // from 38
+
+ // check that these seed are presents
+ {
+ std::string s1("CCGTCACGGGCCTCAGAGAGCAGCGCCTCTTC"); // reference
+ std::string s2("CCCTCACGGGCCTCAGAGAGCAGCGCCTCTTC");
+ std::string s3("CCGTCACGGGTCTCAGAGAGCAGCGCCTCTTC");
+ std::string s4("CCGTCACGGGCCTCAGAGCAGCGCCTCTTCCG"); // <--- indel
+ std::string s5("CCGTCACGGGCCTCAGAGAGCAGTGCCTCTTC");
+
+ std::string s6("CCCTCACGGGTCTCAGAGAGCAGCGCCTCTTC"); // 1-2
+ std::string s7("CCCTCACGGGCCTCAGAGAGCAGTGCCTCTTC"); // 1-4
+ std::string s8("CCGTCACGGGTCTCAGAGAGCAGTGCCTCTTC"); // 2-4
+
+ size_t s1_s = CNVcaller::convert_to_bit_pattern(s1.c_str());
+ size_t s2_s = CNVcaller::convert_to_bit_pattern(s2.c_str());
+ size_t s3_s = CNVcaller::convert_to_bit_pattern(s3.c_str());
+ size_t s4_s = CNVcaller::convert_to_bit_pattern(s4.c_str());
+ size_t s5_s = CNVcaller::convert_to_bit_pattern(s5.c_str());
+ size_t s6_s = CNVcaller::convert_to_bit_pattern(s6.c_str());
+ size_t s7_s = CNVcaller::convert_to_bit_pattern(s7.c_str());
+ size_t s8_s = CNVcaller::convert_to_bit_pattern(s8.c_str());
BOOST_REQUIRE(map.find(s1_s) != map.end());
BOOST_REQUIRE(map.find(s2_s) != map.end());
@@ -143,7 +187,125 @@ BOOST_AUTO_TEST_CASE( CNV_caller_test_add_32bases_pattern_with_variability_range
BOOST_REQUIRE(map.find(s7_s) != map.end());
BOOST_REQUIRE(map.find(s8_s) != map.end());
- // print all seeds in map
+ BOOST_REQUIRE(map.find(s1_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s2_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s3_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s4_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s5_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s6_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s7_s)->second.get(0).padding_left == 8);
+ BOOST_REQUIRE(map.find(s8_s)->second.get(0).padding_left == 8);
+ }
+
+ // from 46
+ // check that these seed are presents
+ {
+ std::string s1("GGCCTCAGAGAGCAGCGCCTCTTCCGACGGCC"); // reference
+ std::string s2("GGTCTCAGAGAGCAGCGCCTCTTCCGACGGCC");
+ std::string s3("GGCCTCAGAGCAGCGCCTCTTCCGACGGCCCC"); // <--- indel
+ std::string s4("GGCCTCAGAGAGCAGTGCCTCTTCCGACGGCC");
+ std::string s5("GGCCTCAGAGAGCAGCGCCTCTTCCGACAGCC");
+
+ std::string s6("GGTCTCAGAGAGCAGCGCCTCTTCCGACAGCC"); // 1-3
+ std::string s7("GGTCTCAGAGAGCAGTGCCTCTTCCGACGGCC"); // 1-4
+ std::string s8("GGCCTCAGAGAGCAGCGCCTCTTCCGACAGCC"); // 3-4
+
+ size_t s1_s = CNVcaller::convert_to_bit_pattern(s1.c_str());
+ size_t s2_s = CNVcaller::convert_to_bit_pattern(s2.c_str());
+ size_t s3_s = CNVcaller::convert_to_bit_pattern(s3.c_str());
+ size_t s4_s = CNVcaller::convert_to_bit_pattern(s4.c_str());
+ size_t s5_s = CNVcaller::convert_to_bit_pattern(s5.c_str());
+ size_t s6_s = CNVcaller::convert_to_bit_pattern(s6.c_str());
+ size_t s7_s = CNVcaller::convert_to_bit_pattern(s7.c_str());
+ size_t s8_s = CNVcaller::convert_to_bit_pattern(s8.c_str());
+
+ BOOST_REQUIRE(map.find(s1_s) != map.end());
+ BOOST_REQUIRE(map.find(s2_s) != map.end());
+ BOOST_REQUIRE(map.find(s3_s) != map.end());
+ BOOST_REQUIRE(map.find(s4_s) != map.end());
+ BOOST_REQUIRE(map.find(s5_s) != map.end());
+ BOOST_REQUIRE(map.find(s6_s) != map.end());
+ BOOST_REQUIRE(map.find(s7_s) != map.end());
+ BOOST_REQUIRE(map.find(s8_s) != map.end());
+
+ BOOST_REQUIRE(map.find(s1_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s2_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s3_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s4_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s5_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s6_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s7_s)->second.get(0).padding_left == 16);
+ BOOST_REQUIRE(map.find(s8_s)->second.get(0).padding_left == 16);
+ }
+
+ // from 54
+ // check that these seed are presents
+ {
+ std::string s1("AGAGCAGCGCCTCTTCCGACGGCCCCCATCCA"); // reference
+ std::string s3("AGAGCAGTGCCTCTTCCGACGGCCCCCATCCA");
+ std::string s4("AGAGCAGCGCCTCTTCCGACAGCCCCCATCCA");
+ std::string s5("AGAGCAGCGCCTCTTCCGACGGCCTCCATCCA");
+
+ std::string s6("AGAGCAGTGCCTCTTCCGACAGCCCCCATCCA"); // 1-2
+ std::string s7("AGAGCAGTGCCTCTTCCGACGGCCTCCATCCA"); // 1-3
+ std::string s8("AGAGCAGCGCCTCTTCCGACAGCCTCCATCCA"); // 2-3
+
+ std::string s9("AGAGCAGTGCCTCTTCCGACAGCCTCCATCCA"); // 1-2-3
+
+ size_t s1_s = CNVcaller::convert_to_bit_pattern(s1.c_str());
+ size_t s3_s = CNVcaller::convert_to_bit_pattern(s3.c_str());
+ size_t s4_s = CNVcaller::convert_to_bit_pattern(s4.c_str());
+ size_t s5_s = CNVcaller::convert_to_bit_pattern(s5.c_str());
+ size_t s6_s = CNVcaller::convert_to_bit_pattern(s6.c_str());
+ size_t s7_s = CNVcaller::convert_to_bit_pattern(s7.c_str());
+ size_t s8_s = CNVcaller::convert_to_bit_pattern(s8.c_str());
+ size_t s9_s = CNVcaller::convert_to_bit_pattern(s9.c_str());
+
+ BOOST_REQUIRE(map.find(s1_s) != map.end());
+ BOOST_REQUIRE(map.find(s3_s) != map.end());
+ BOOST_REQUIRE(map.find(s4_s) != map.end());
+ BOOST_REQUIRE(map.find(s5_s) != map.end());
+ BOOST_REQUIRE(map.find(s6_s) != map.end());
+ BOOST_REQUIRE(map.find(s7_s) != map.end());
+ BOOST_REQUIRE(map.find(s8_s) != map.end());
+ BOOST_REQUIRE(map.find(s9_s) != map.end());
+
+ BOOST_REQUIRE(map.find(s1_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s3_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s4_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s5_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s6_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s7_s)->second.get(0).padding_left == 24);
+ BOOST_REQUIRE(map.find(s8_s)->second.get(0).padding_left == 24);
+ }
+
+ // from 62
+ // check that these seed are presents
+ {
+ std::string s1("GCCTCTTCCGACGGCCCCCATCCAGTCATCAC"); // reference
+ std::string s2("GCCTCTTCCGACAGCCCCCATCCAGTCATCAC");
+ std::string s3("GCCTCTTCCGACGGCCTCCATCCAGTCATCAC");
+
+ std::string s4("GCCTCTTCCGACAGCCTCCATCCAGTCATCAC"); // 1-2
+
+ size_t s1_s = CNVcaller::convert_to_bit_pattern(s1.c_str());
+ size_t s2_s = CNVcaller::convert_to_bit_pattern(s2.c_str());
+ size_t s3_s = CNVcaller::convert_to_bit_pattern(s3.c_str());
+ size_t s4_s = CNVcaller::convert_to_bit_pattern(s4.c_str());
+
+ BOOST_REQUIRE(map.find(s1_s) != map.end());
+ BOOST_REQUIRE(map.find(s2_s) != map.end());
+ BOOST_REQUIRE(map.find(s3_s) != map.end());
+ BOOST_REQUIRE(map.find(s4_s) != map.end());
+
+ BOOST_REQUIRE(map.find(s1_s)->second.get(0).padding_left == 32);
+ BOOST_REQUIRE(map.find(s2_s)->second.get(0).padding_left == 32);
+ BOOST_REQUIRE(map.find(s3_s)->second.get(0).padding_left == 32);
+ BOOST_REQUIRE(map.find(s4_s)->second.get(0).padding_left == 32);
+
+ }
+
+ // print all seeds in map and check that align with the expected
/* auto i = map.begin();
@@ -847,7 +1009,7 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
seqs.get<0>(0).get<0>(j).get<0>(i),
un_seed.get(j),
calls_to_test.get(j),
- -6,
+ 6,
seqs.get<0>(0).get<0>(j).get<1>(i),
seqs.get<0>(0).get<0>(j).get<2>(i),
seed_pos,
@@ -880,7 +1042,7 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
int offset = 0;
// call id 1
- add_32bases_pattern_with_variability(2,195452561,reference,map,"test/gnomad_variants_variability_test.vcf.gz",1,195452561-195451991+32,NO_PADDING);
+ add_32bases_pattern_with_variability_range_shift(2,195452561,195452561+32,8,reference,map,"test/gnomad_variants_variability_test.vcf.gz",76,195452561-195451992+32,NO_PADDING);
auto i = map.begin();
@@ -900,6 +1062,14 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
for (int j = 0 ; j < seqs.get<0>(i).size() ; j++) {
for (int k = 0 ; k < seqs.get<0>(i).get<0>(j).size() ; k++) {
auto & pair = seqs.get<0>(i).get<0>(j).get<0>(k);
+
+ if (pair.get<1>(3) != 'A' && pair.get<1>(3) != 'C' && pair.get<1>(3) != 'T' && pair.get<1>(3) != 'G')
+ {
+ seqs.get<0>(i).get<0>(j).remove(k);
+ k--;
+ std::cout << "Error: " << pair.get<1>(3) << std::endl;
+ }
+
BOOST_REQUIRE((pair.get<1>(3) == 'A' || pair.get<1>(3) == 'C' || pair.get<1>(3) == 'T' || pair.get<1>(3) == 'G'));
}
}
@@ -910,9 +1080,10 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
// Find all reference points
openfpm::vector<openfpm::vector<aggregate<short int, int, int>>> un_seed;
-// un_seed.resize(calls_to_test.size());
+ un_seed.resize(calls_to_test.size());
// get_reference_point(seed,reference,un_seed);
- un_seed.load("test/unseed_3calls_cnv_variability");
+// un_seed.save("test/MUC20_cnv_variability");
+ un_seed.load("test/MUC20_cnv_variability");
// Calculate evidence for each of the read pair
@@ -941,7 +1112,7 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
seqs.get<0>(0).get<0>(j).get<0>(i),
un_seed.get(j),
calls_to_test.get(j),
- -6,
+ 6,
seqs.get<0>(0).get<0>(j).get<1>(i),
seqs.get<0>(0).get<0>(j).get<2>(i),
seed_pos,
@@ -949,9 +1120,6 @@ BOOST_AUTO_TEST_CASE( CNV_caller_evidence_probability_for_each_event_with_variab
bp_insert_prob);
std::cout << "Probability: " << eml3 << std::endl;
-
- BOOST_REQUIRE(eml3(59) > -20.0 );
- BOOST_REQUIRE(eml3(121) < -20.0);
}
}
--
GitLab