fixing off by one bug in reporting alleles in the VCF file

mgymrek · Jul 28, 2014 · 5fd8768 · 5fd8768
1 parent 61e937a
commit 5fd8768
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 4 deletions.
diff --git a/src/Genotyper.cpp b/src/Genotyper.cpp
@@ -119,7 +119,8 @@ void Genotyper::GetAlleles(const list<AlignedRead>& aligned_reads,
   for (list<AlignedRead>::const_iterator it = aligned_reads.begin();
        it != aligned_reads.end(); it++) {
     if (it->mate) continue;
-    if (it->diffFromRef != 0 && std::find(alleles->begin(), alleles->end(), it->diffFromRef) == alleles->end()) {
+    if (it->diffFromRef != 0 && std::find(alleles->begin(), alleles->end(), it->diffFromRef) == alleles->end() &&
+	it->diffFromRef+it->refCopyNum*it->period > 0) {
       alleles->push_back(it->diffFromRef);
     }
   }
@@ -378,7 +379,7 @@ void Genotyper::Genotype(const list<AlignedRead>& read_list) {
   str_record.stop = read_list.front().msEnd;
   str_record.repseq = read_list.front().repseq;
   str_record.refcopy = static_cast<float>((read_list.front().msEnd-
-					   read_list.front().msStart))/
+					   read_list.front().msStart+1))/
     static_cast<float>(read_list.front().period);
   if (ref_nucleotides->find
       (pair<string,int>(str_record.chrom, str_record.start))

diff --git a/src/VCFWriter.cpp b/src/VCFWriter.cpp
@@ -174,7 +174,7 @@ void VCFWriter::WriteRecord(const STRRecord& str_record) {
   output_stream << "END=" << str_record.stop << ";"
 		<< "MOTIF=" << str_record.repseq << ";"
 		<< "REF=" << str_record.refcopy << ";"
-		<< "RL=" << str_record.stop - str_record.start<< ";"
+		<< "RL=" << str_record.stop - str_record.start + 1 << ";"
 		<< "RPA=" << repeats_per_allele.str() << ";"
 		<< "RU=" << str_record.repseq_in_ref << ";"
 		<< "VT=STR" << "\t";

diff --git a/src/main_allelotype.cpp b/src/main_allelotype.cpp
@@ -430,7 +430,7 @@ int main(int argc, char* argv[]) {
 	ref_str.stop         = str_end-extend;
 	ref_str.chrom        = chrom;
 	reference_strs.push_back(ref_str);
-	string refnuc = ref_record.nucleotides.substr(extend, ref_record.nucleotides.length()-2*extend);
+	string refnuc = ref_record.nucleotides.substr(extend-1, ref_record.nucleotides.length()-2*extend+1);
 	string repseq_in_ref = ref_record.nucleotides.substr(extend, repseq.size());
 	pair<string, int> locus = pair<string,int>(chrom, start);
 	ref_nucleotides.insert(pair< pair<string, int>, string>(locus, refnuc));