373 |
|
|
374 |
|
int anchor_mismatch = 0; |
375 |
|
|
376 |
+ |
|
377 |
+ |
void parseSegReadName(char* name, char*& name_tags, bool strip_slash, |
378 |
+ |
bool &end, unsigned int &seg_offset, unsigned int& seg_num, |
379 |
+ |
unsigned int & num_segs) { |
380 |
+ |
char* pipe = strrchr(name, '|'); |
381 |
+ |
if (pipe) |
382 |
+ |
{ |
383 |
+ |
if (name_tags) |
384 |
+ |
strcpy(name_tags, pipe); |
385 |
+ |
|
386 |
+ |
char* tag_buf = pipe + 1; |
387 |
+ |
if (strchr(tag_buf, ':')) |
388 |
+ |
{ |
389 |
+ |
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
390 |
+ |
if (seg_num + 1 == num_segs) |
391 |
+ |
end = true; |
392 |
+ |
else |
393 |
+ |
end = false; |
394 |
+ |
} |
395 |
+ |
|
396 |
+ |
*pipe = 0; |
397 |
+ |
} |
398 |
+ |
|
399 |
+ |
// Stripping the slash and number following it gives the insert name |
400 |
+ |
char* slash = strrchr(name, '/'); |
401 |
+ |
if (strip_slash && slash) |
402 |
+ |
*slash = 0; |
403 |
+ |
} |
404 |
+ |
|
405 |
+ |
|
406 |
|
bool SplicedBowtieHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
407 |
|
BowtieHit& bh, |
408 |
|
bool strip_slash, |
456 |
|
unsigned int num_segs = 0; |
457 |
|
|
458 |
|
// Copy the tag out of the name field before we might wipe it out |
459 |
< |
char* pipe = strrchr(name, '|'); |
430 |
< |
if (pipe) |
431 |
< |
{ |
432 |
< |
if (name_tags) |
433 |
< |
strcpy(name_tags, pipe); |
459 |
> |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
460 |
|
|
435 |
– |
char* tag_buf = pipe + 1; |
436 |
– |
if (strchr(tag_buf, ':')) |
437 |
– |
{ |
438 |
– |
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
439 |
– |
if (seg_num + 1 == num_segs) |
440 |
– |
end = true; |
441 |
– |
else |
442 |
– |
end = false; |
443 |
– |
} |
444 |
– |
|
445 |
– |
*pipe = 0; |
446 |
– |
} |
447 |
– |
// Stripping the slash and number following it gives the insert name |
448 |
– |
char* slash = strrchr(name, '/'); |
449 |
– |
if (strip_slash && slash) |
450 |
– |
*slash = 0; |
451 |
– |
|
452 |
– |
//int read_len = strlen(sequence); |
453 |
– |
|
461 |
|
// Add this alignment to the table of hits for this half of the |
462 |
|
// Bowtie map |
463 |
|
|
719 |
|
return false; |
720 |
|
} |
721 |
|
|
722 |
+ |
|
723 |
+ |
char parseCigar(vector<CigarOp>& cigar, const char* cigar_str, |
724 |
+ |
bool &spliced_alignment) { |
725 |
+ |
// Mostly pilfered direct from the SAM tools: |
726 |
+ |
const char* p_cig = cigar_str; |
727 |
+ |
while (*p_cig) |
728 |
+ |
{ |
729 |
+ |
char* t; |
730 |
+ |
int length = (int)strtol(p_cig, &t, 10); |
731 |
+ |
if (length <= 0) |
732 |
+ |
{ |
733 |
+ |
//fprintf (stderr, "CIGAR op has zero length\n"); |
734 |
+ |
return false; |
735 |
+ |
} |
736 |
+ |
char op_char = toupper(*t); |
737 |
+ |
CigarOpCode opcode; |
738 |
+ |
if (op_char == 'M') opcode = MATCH; |
739 |
+ |
else if (op_char == 'I') opcode = INS; |
740 |
+ |
else if (op_char == 'D') opcode = DEL; |
741 |
+ |
else if (op_char == 'N') |
742 |
+ |
{ |
743 |
+ |
if (length > max_report_intron_length) |
744 |
+ |
return false; |
745 |
+ |
opcode = REF_SKIP; |
746 |
+ |
spliced_alignment = true; |
747 |
+ |
} |
748 |
+ |
else if (op_char == 'S') opcode = SOFT_CLIP; |
749 |
+ |
else if (op_char == 'H') opcode = HARD_CLIP; |
750 |
+ |
else if (op_char == 'P') opcode = PAD; |
751 |
+ |
else |
752 |
+ |
{ |
753 |
+ |
fprintf (stderr, "invalid CIGAR operation\n"); |
754 |
+ |
return false; |
755 |
+ |
} |
756 |
+ |
p_cig = t + 1; |
757 |
+ |
//i += length; |
758 |
+ |
cigar.push_back(CigarOp(opcode, length)); |
759 |
+ |
} |
760 |
+ |
if (*p_cig) { |
761 |
+ |
fprintf (stderr, "unmatched CIGAR operation\n"); |
762 |
+ |
return *p_cig; |
763 |
+ |
} |
764 |
+ |
return 0; |
765 |
+ |
} |
766 |
+ |
|
767 |
+ |
void getSAMmismatches(char* &buf, vector<CigarOp>& cigar, |
768 |
+ |
int* mismatches, int& num_mismatches, int& sam_nm, bool& antisense_splice) { |
769 |
+ |
const char* tag_buf = buf; |
770 |
+ |
sam_nm=0; |
771 |
+ |
num_mismatches=0; |
772 |
+ |
while((tag_buf = get_token((char**)&buf,"\t"))) |
773 |
+ |
{ |
774 |
+ |
vector<string> tuple_fields; |
775 |
+ |
tokenize(tag_buf,":", tuple_fields); |
776 |
+ |
if (tuple_fields.size() == 3) |
777 |
+ |
{ |
778 |
+ |
if (tuple_fields[0] == "XS") |
779 |
+ |
{ |
780 |
+ |
if (tuple_fields[2] == "-") |
781 |
+ |
antisense_splice = true; |
782 |
+ |
} |
783 |
+ |
else if (tuple_fields[0] == "NM") |
784 |
+ |
{ |
785 |
+ |
sam_nm = atoi(tuple_fields[2].c_str()); |
786 |
+ |
} |
787 |
+ |
else if (tuple_fields[0] == "NS") |
788 |
+ |
{ |
789 |
+ |
//ignored for now |
790 |
+ |
} |
791 |
+ |
else if (tuple_fields[0] == "MD") |
792 |
+ |
{ |
793 |
+ |
const char* p=tuple_fields[2].c_str(); |
794 |
+ |
int bi=0; //base offset position in the read |
795 |
+ |
while (*p != 0) { |
796 |
+ |
if (isdigit(*p)) { |
797 |
+ |
int v=atoi(p); |
798 |
+ |
do { p++; } while (isdigit(*p)); |
799 |
+ |
bi+=v; |
800 |
+ |
} |
801 |
+ |
while (isupper(*p)) { |
802 |
+ |
p++; |
803 |
+ |
mismatches[num_mismatches]=bi; |
804 |
+ |
num_mismatches++; |
805 |
+ |
bi++; |
806 |
+ |
} |
807 |
+ |
if (*p=='^') { //reference deletion |
808 |
+ |
p++; |
809 |
+ |
while (isupper(*p)) { //insert read bases |
810 |
+ |
p++; bi++; |
811 |
+ |
} |
812 |
+ |
} |
813 |
+ |
} |
814 |
+ |
} |
815 |
+ |
else |
816 |
+ |
{ |
817 |
+ |
//fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str()); |
818 |
+ |
//return false; |
819 |
+ |
} |
820 |
+ |
} |
821 |
+ |
} |
822 |
+ |
/* By convention,the NM field of the SAM record |
823 |
+ |
* counts an insertion or deletion. I dont' think |
824 |
+ |
* we want the mismatch count in the BowtieHit |
825 |
+ |
* record to reflect this. Therefore, subtract out |
826 |
+ |
* the mismatches due to in/dels |
827 |
+ |
*/ |
828 |
+ |
for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){ |
829 |
+ |
if(itr->opcode == INS){ |
830 |
+ |
sam_nm -= itr->length; |
831 |
+ |
} |
832 |
+ |
if(itr->opcode == DEL){ |
833 |
+ |
sam_nm -= itr->length; |
834 |
+ |
} |
835 |
+ |
} |
836 |
+ |
} |
837 |
+ |
|
838 |
|
bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
839 |
|
BowtieHit& bh, |
840 |
|
bool strip_slash, |
896 |
|
unsigned int num_segs = 0; |
897 |
|
|
898 |
|
// Copy the tag out of the name field before we might wipe it out |
899 |
< |
char* pipe = strrchr(name, '|'); |
777 |
< |
if (pipe) |
778 |
< |
{ |
779 |
< |
if (name_tags) |
780 |
< |
strcpy(name_tags, pipe); |
781 |
< |
|
782 |
< |
char* tag_buf = pipe + 1; |
783 |
< |
if (strchr(tag_buf, ':')) |
784 |
< |
{ |
785 |
< |
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
786 |
< |
if (seg_num + 1 == num_segs) |
787 |
< |
end = true; |
788 |
< |
else |
789 |
< |
end = false; |
790 |
< |
} |
791 |
< |
|
792 |
< |
*pipe = 0; |
793 |
< |
} |
794 |
< |
|
795 |
< |
// Stripping the slash and number following it gives the insert name |
796 |
< |
char* slash = strrchr(name, '/'); |
797 |
< |
if (strip_slash && slash) |
798 |
< |
*slash = 0; |
899 |
> |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
900 |
|
|
800 |
– |
char* p_cig = cigar_str; |
801 |
– |
//int len = strlen(sequence); |
901 |
|
vector<CigarOp> cigar; |
902 |
|
bool spliced_alignment = false; |
903 |
< |
// Mostly pilfered direct from the SAM tools: |
904 |
< |
while (*p_cig) |
905 |
< |
{ |
807 |
< |
char* t; |
808 |
< |
int length = (int)strtol(p_cig, &t, 10); |
809 |
< |
if (length <= 0) |
810 |
< |
{ |
811 |
< |
//fprintf (stderr, "CIGAR op has zero length\n"); |
812 |
< |
return false; |
813 |
< |
} |
814 |
< |
char op_char = toupper(*t); |
815 |
< |
CigarOpCode opcode; |
816 |
< |
if (op_char == 'M') opcode = MATCH; |
817 |
< |
else if (op_char == 'I') opcode = INS; |
818 |
< |
else if (op_char == 'D') opcode = DEL; |
819 |
< |
else if (op_char == 'N') |
820 |
< |
{ |
821 |
< |
if (length > max_report_intron_length) |
822 |
< |
return false; |
823 |
< |
opcode = REF_SKIP; |
824 |
< |
spliced_alignment = true; |
825 |
< |
} |
826 |
< |
else if (op_char == 'S') opcode = SOFT_CLIP; |
827 |
< |
else if (op_char == 'H') opcode = HARD_CLIP; |
828 |
< |
else if (op_char == 'P') opcode = PAD; |
829 |
< |
else |
830 |
< |
{ |
831 |
< |
fprintf (stderr, "invalid CIGAR operation\n"); |
832 |
< |
return false; |
833 |
< |
} |
834 |
< |
p_cig = t + 1; |
835 |
< |
//i += length; |
836 |
< |
cigar.push_back(CigarOp(opcode, length)); |
837 |
< |
} |
838 |
< |
if (*p_cig) |
839 |
< |
{ |
840 |
< |
fprintf (stderr, "unmatched CIGAR operation\n"); |
841 |
< |
return false; |
842 |
< |
} |
903 |
> |
if (parseCigar(cigar, cigar_str, spliced_alignment)) { |
904 |
> |
return false; |
905 |
> |
} |
906 |
|
|
907 |
|
//vector<string> attributes; |
908 |
|
//tokenize(tag_buf, " \t",attributes); |
909 |
< |
|
909 |
> |
|
910 |
|
bool antisense_splice = false; |
911 |
< |
unsigned char num_mismatches = 0; |
912 |
< |
unsigned char num_splice_anchor_mismatches = 0; |
913 |
< |
const char* tag_buf = buf; |
914 |
< |
|
852 |
< |
// while((tag_buf = get_token((char**)&buf,"\t"))) |
853 |
< |
// { |
854 |
< |
// vector<string> tuple_fields; |
855 |
< |
// tokenize(tag_buf,":", tuple_fields); |
856 |
< |
// if (tuple_fields.size() == 3) |
857 |
< |
// { |
858 |
< |
// if (tuple_fields[0] == "XS") |
859 |
< |
// { |
860 |
< |
// if (tuple_fields[2] == "-") |
861 |
< |
// antisense_splice = true; |
862 |
< |
// } |
863 |
< |
// else if (tuple_fields[0] == "NM") |
864 |
< |
// { |
865 |
< |
// num_mismatches = atoi(tuple_fields[2].c_str()); |
866 |
< |
// } |
867 |
< |
// else if (tuple_fields[0] == "NS") |
868 |
< |
// { |
869 |
< |
// num_splice_anchor_mismatches = atoi(tuple_fields[2].c_str()); |
870 |
< |
// } |
871 |
< |
// else |
872 |
< |
// { |
873 |
< |
// fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str()); |
874 |
< |
// return false; |
875 |
< |
// } |
876 |
< |
// } |
877 |
< |
// } |
878 |
< |
|
879 |
< |
while((tag_buf = get_token((char**)&buf,"\t"))) |
880 |
< |
{ |
881 |
< |
vector<string> tuple_fields; |
882 |
< |
tokenize(tag_buf,":", tuple_fields); |
883 |
< |
if (tuple_fields.size() == 3) |
884 |
< |
{ |
885 |
< |
if (tuple_fields[0] == "XS") |
886 |
< |
{ |
887 |
< |
if (tuple_fields[2] == "-") |
888 |
< |
antisense_splice = true; |
889 |
< |
} |
890 |
< |
else if (tuple_fields[0] == "NM") |
891 |
< |
{ |
892 |
< |
num_mismatches = atoi(tuple_fields[2].c_str()); |
893 |
< |
} |
894 |
< |
else if (tuple_fields[0] == "NS") |
895 |
< |
{ |
896 |
< |
//ignored for now |
897 |
< |
} |
898 |
< |
else |
899 |
< |
{ |
900 |
< |
//fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str()); |
901 |
< |
//return false; |
902 |
< |
} |
903 |
< |
} |
904 |
< |
} |
911 |
> |
int num_mismatches = 0; |
912 |
> |
int sam_nm = 0; //the value of the NM tag (edit distance) |
913 |
> |
int mismatches[1024]; |
914 |
> |
getSAMmismatches(buf, cigar, mismatches, num_mismatches, sam_nm, antisense_splice); |
915 |
|
|
906 |
– |
|
907 |
– |
/* |
908 |
– |
* By convention,the NM field of the SAM record |
909 |
– |
* counts an insertion or deletion. I dont' think |
910 |
– |
* we want the mismatch count in the BowtieHit |
911 |
– |
* record to reflect this. Therefore, subtract out |
912 |
– |
* the mismatches due to in/dels |
913 |
– |
*/ |
914 |
– |
for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){ |
915 |
– |
if(itr->opcode == INS){ |
916 |
– |
num_mismatches -= itr->length; |
917 |
– |
} |
918 |
– |
if(itr->opcode == DEL){ |
919 |
– |
num_mismatches -= itr->length; |
920 |
– |
} |
921 |
– |
} |
922 |
– |
// vector<string> toks; |
923 |
– |
// tokenize(tag_buf, ":", toks); |
916 |
|
if (spliced_alignment) |
917 |
|
{ |
918 |
|
bh = create_hit(name, |
922 |
|
sam_flag & 0x0010, |
923 |
|
antisense_splice, |
924 |
|
num_mismatches, |
925 |
< |
num_splice_anchor_mismatches, |
925 |
> |
0, |
926 |
|
end); |
935 |
– |
return true; |
936 |
– |
|
927 |
|
} |
928 |
|
else |
929 |
|
{ |
930 |
|
//assert(cigar.size() == 1 && cigar[0].opcode == MATCH); |
941 |
– |
|
931 |
|
bh = create_hit(name, |
932 |
|
text_name, |
933 |
|
text_offset - 1, // SAM files are 1-indexed |
937 |
|
num_mismatches, |
938 |
|
0, |
939 |
|
end); |
951 |
– |
return true; |
940 |
|
} |
941 |
< |
|
942 |
< |
return false; |
941 |
> |
return true; |
942 |
> |
} |
943 |
> |
|
944 |
> |
|
945 |
> |
void spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, |
946 |
> |
int8_t l_len, int mlen, int8_t r_len, CigarOpCode opcode) { |
947 |
> |
//merge the original 'cigar' with the new insert/gap operation |
948 |
> |
// at position l_len into splcigar; |
949 |
> |
//find the original opcode that croses l_len location |
950 |
> |
int ofs=0; |
951 |
> |
bool found=false; |
952 |
> |
for (size_t c = 0 ; c < cigar.size(); ++c) |
953 |
> |
{ |
954 |
> |
ofs+=cigar[c].length; |
955 |
> |
if (!found) { |
956 |
> |
if (ofs>=l_len) { |
957 |
> |
found=true; |
958 |
> |
//inject new code here |
959 |
> |
} |
960 |
> |
else { |
961 |
> |
//not found yet, just copy these codes |
962 |
> |
splcigar.push_back(cigar[c]); |
963 |
> |
} |
964 |
> |
} |
965 |
> |
else { //found yet |
966 |
> |
|
967 |
> |
} |
968 |
> |
} |
969 |
|
} |
970 |
|
|
971 |
|
bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
985 |
|
if (bwt_buf[0] == '@') |
986 |
|
return false; |
987 |
|
|
974 |
– |
//char orientation; |
975 |
– |
/* |
976 |
– |
char mismatches[buf_size]; |
977 |
– |
|
978 |
– |
//char* orientation_str = get_token((char**)&buf,"\t"); |
979 |
– |
|
980 |
– |
mismatches[0] = 0; |
981 |
– |
char* mismatches_str = get_token((char**)&buf, "\t"); |
982 |
– |
if (mismatches_str) |
983 |
– |
strcpy(mismatches, mismatches_str); |
984 |
– |
|
985 |
– |
orientation = orientation_str[0]; |
986 |
– |
text_offset = atoi(text_offset_str); |
987 |
– |
*/ |
988 |
– |
|
988 |
|
char* buf = bwt_buf; |
989 |
|
char* name = get_token((char**)&buf,"\t"); |
990 |
|
char* sam_flag_str = get_token((char**)&buf,"\t"); |
995 |
|
const char* mate_ref_str = get_token((char**)&buf,"\t"); |
996 |
|
const char* mate_pos_str = get_token((char**)&buf,"\t"); |
997 |
|
const char* inferred_insert_sz_str = get_token((char**)&buf,"\t"); |
998 |
< |
|
998 |
> |
int num_mismatches=0; |
999 |
> |
int mismatches[1024]; //list of 0-based mismatch positions in this read |
1000 |
> |
//parsed from SAM's MD:Z: tag |
1001 |
|
const char* seq_str = get_token((char**)&buf,"\t"); |
1002 |
|
if (seq) |
1003 |
|
strcpy(seq, seq_str); |
1025 |
|
|
1026 |
|
int sam_flag = atoi(sam_flag_str); |
1027 |
|
int text_offset = atoi(text_offset_str); |
1028 |
< |
|
1028 |
< |
//############################################## |
1028 |
> |
text_offset--; //SAM is 1-based, Bowtie is 0-based |
1029 |
|
bool end = true; |
1030 |
|
unsigned int seg_offset = 0; |
1031 |
|
unsigned int seg_num = 0; |
1032 |
|
unsigned int num_segs = 0; |
1033 |
|
|
1034 |
|
// Copy the tag out of the name field before we might wipe it out |
1035 |
< |
char* pipe = strrchr(name, '|'); |
1036 |
< |
if (pipe) |
1037 |
< |
{ |
1038 |
< |
if (name_tags) |
1039 |
< |
strcpy(name_tags, pipe); |
1035 |
> |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
1036 |
|
|
1037 |
< |
char* tag_buf = pipe + 1; |
1038 |
< |
if (strchr(tag_buf, ':')) |
1039 |
< |
{ |
1040 |
< |
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
1041 |
< |
if (seg_num + 1 == num_segs) |
1046 |
< |
end = true; |
1047 |
< |
else |
1048 |
< |
end = false; |
1049 |
< |
} |
1037 |
> |
vector<CigarOp> samcigar; |
1038 |
> |
bool spliced_alignment = false; |
1039 |
> |
if (parseCigar(samcigar, cigar_str, spliced_alignment)) { |
1040 |
> |
return false; |
1041 |
> |
} |
1042 |
|
|
1043 |
< |
*pipe = 0; |
1044 |
< |
} |
1045 |
< |
// Stripping the slash and number following it gives the insert name |
1054 |
< |
char* slash = strrchr(name, '/'); |
1055 |
< |
if (strip_slash && slash) |
1056 |
< |
*slash = 0; |
1043 |
> |
bool antisense_splice = false; |
1044 |
> |
int sam_nm = 0; |
1045 |
> |
getSAMmismatches(buf, samcigar, mismatches, num_mismatches, sam_nm, antisense_splice); |
1046 |
|
|
1047 |
< |
//int read_len = strlen(sequence); |
1047 |
> |
//############################################## |
1048 |
|
|
1049 |
|
// Add this alignment to the table of hits for this half of the |
1050 |
|
// Bowtie map |
1076 |
|
if (splice_toks.size() != 2) |
1077 |
|
{ |
1078 |
|
fprintf(stderr, "Warning: found malformed splice record, skipping:\n"); |
1079 |
< |
//fprintf(stderr, "%s (token: %s)\n", text_name, |
1079 |
> |
//fprintf(stderr, "\t%s (token: %s)\n", text_name, |
1080 |
|
// toks[num_extra_toks + splice_field].c_str()); |
1081 |
|
return false; |
1082 |
|
} |
1083 |
< |
|
1083 |
> |
int spl_num_mismatches=0; |
1084 |
|
// |
1085 |
|
// check for an insertion hit |
1086 |
|
// |
1087 |
|
if(toks[num_extra_toks + junction_type_field] == "ins") |
1088 |
|
{ |
1089 |
|
int8_t spliced_read_len = strlen(seq_str); |
1090 |
< |
/* |
1091 |
< |
* The 0-based position of the left edge of the alignment. Note that this |
1092 |
< |
* value may need to be futher corrected to account for the presence of |
1104 |
< |
* of the insertion. |
1105 |
< |
*/ |
1090 |
> |
// The 0-based position of the left edge of the alignment. Note that this |
1091 |
> |
// value may need to be futher corrected to account for the presence of |
1092 |
> |
// of the insertion. |
1093 |
|
uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
1094 |
|
uint32_t right = left + spliced_read_len - 1; |
1095 |
|
|
1096 |
< |
/* |
1110 |
< |
* The 0-based position of the last genomic sequence before the insertion |
1111 |
< |
*/ |
1096 |
> |
// The 0-based position of the last genomic sequence before the insertion |
1097 |
|
uint32_t left_splice_pos = atoi(splice_toks[0].c_str()); |
1098 |
|
|
1099 |
|
string insertedSequence = splice_toks[1]; |
1100 |
< |
/* |
1101 |
< |
* The 0-based position of the first genomic sequence after teh insertion |
1117 |
< |
*/ |
1100 |
> |
// The 0-based position of the first genomic sequence after the insertion |
1101 |
> |
|
1102 |
|
uint32_t right_splice_pos = left_splice_pos + 1; |
1103 |
|
if(left > left_splice_pos){ |
1104 |
|
/* |
1141 |
|
return false; |
1142 |
|
} |
1143 |
|
|
1144 |
< |
char* pch = strtok( mismatches, ","); |
1145 |
< |
unsigned char num_mismatches = 0; |
1144 |
> |
//char* pch = strtok( mismatches, ","); |
1145 |
> |
//unsigned char num_mismatches = 0; |
1146 |
> |
//text_offset holds the left end of the alignment, |
1147 |
> |
//RELATIVE TO the start of the contig |
1148 |
|
|
1149 |
< |
/* |
1150 |
< |
* remember that text_offset holds the left end of the |
1165 |
< |
* alignment, relative to the start of the contig |
1166 |
< |
*/ |
1167 |
< |
|
1168 |
< |
/* |
1169 |
< |
* The 0-based relative position of the left-most character |
1170 |
< |
* before the insertion in the contig |
1171 |
< |
*/ |
1149 |
> |
//The 0-based relative position of the left-most character |
1150 |
> |
//before the insertion in the contig |
1151 |
|
int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str()); |
1152 |
< |
while (pch != NULL) |
1153 |
< |
{ |
1154 |
< |
char* colon = strchr(pch, ':'); |
1155 |
< |
if (colon) |
1156 |
< |
{ |
1157 |
< |
*colon = 0; |
1158 |
< |
int mismatch_pos = atoi(pch); |
1159 |
< |
|
1160 |
< |
/* |
1161 |
< |
* for reversely mapped reads, |
1162 |
< |
* find the correct mismatched position. |
1163 |
< |
*/ |
1164 |
< |
if (sam_flag & 0x0010){ |
1165 |
< |
mismatch_pos = spliced_read_len - mismatch_pos - 1; |
1166 |
< |
} |
1188 |
< |
|
1189 |
< |
/* |
1190 |
< |
* Only count mismatches outside of the insertion region |
1191 |
< |
* If there is a mismatch within the insertion, |
1192 |
< |
* disallow this hit |
1193 |
< |
*/ |
1194 |
< |
if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){ |
1195 |
< |
num_mismatches++; |
1196 |
< |
}else{ |
1197 |
< |
return false; |
1198 |
< |
} |
1152 |
> |
for (int i=0;i<num_mismatches;i++) { |
1153 |
> |
int mismatch_pos = mismatches[i]; |
1154 |
> |
// for reversely mapped reads, |
1155 |
> |
//find the correct mismatched position. |
1156 |
> |
if (sam_flag & 0x0010){ |
1157 |
> |
mismatch_pos = spliced_read_len - mismatch_pos - 1; |
1158 |
> |
} |
1159 |
> |
|
1160 |
> |
//Only count mismatches outside of the insertion region |
1161 |
> |
// If there is a mismatch within the insertion, |
1162 |
> |
// disallow this hit |
1163 |
> |
if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){ |
1164 |
> |
spl_num_mismatches++; |
1165 |
> |
}else{ |
1166 |
> |
return false; |
1167 |
|
} |
1200 |
– |
pch = strtok (NULL, ","); |
1168 |
|
} |
1169 |
|
|
1170 |
< |
|
1171 |
< |
vector<CigarOp> cigar; |
1172 |
< |
cigar.push_back(CigarOp(MATCH, left_match_length)); |
1173 |
< |
cigar.push_back(CigarOp(INS, insertion_match_length)); |
1174 |
< |
cigar.push_back(CigarOp(MATCH, right_match_length)); |
1175 |
< |
|
1170 |
> |
//TODO: merge with the original cigar string |
1171 |
> |
vector<CigarOp> splcigar; |
1172 |
> |
spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS); |
1173 |
> |
/* |
1174 |
> |
splcigar.push_back(CigarOp(MATCH, left_match_length)); |
1175 |
> |
splcigar.push_back(CigarOp(INS, insertion_match_length)); |
1176 |
> |
splcigar.push_back(CigarOp(MATCH, right_match_length)); |
1177 |
> |
*/ |
1178 |
|
/* |
1179 |
|
* For now, disallow hits that don't span |
1180 |
|
* the insertion. If we allow these types of hits, |
1188 |
|
bh = create_hit(name, |
1189 |
|
contig, |
1190 |
|
left, |
1191 |
< |
cigar, |
1192 |
< |
sam_flag & 0x0010, //####### |
1191 |
> |
//splcigar, |
1192 |
> |
splcigar, |
1193 |
> |
sam_flag & 0x0010, |
1194 |
|
junction_strand == "rev", |
1195 |
< |
num_mismatches, |
1195 |
> |
spl_num_mismatches, |
1196 |
|
0, |
1197 |
|
end); |
1198 |
|
return true; |
1199 |
< |
} |
1230 |
< |
|
1199 |
> |
} //"ins" |
1200 |
|
else |
1201 |
|
{ |
1202 |
|
uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
1210 |
|
|
1211 |
|
if ((sam_flag & 0x0010) == 0) //###### |
1212 |
|
{ |
1213 |
< |
if (left_splice_pos + seg_offset < _anchor_length){ |
1214 |
< |
return false; |
1246 |
< |
} |
1213 |
> |
if (left_splice_pos + seg_offset < _anchor_length) |
1214 |
> |
return false; |
1215 |
|
} |
1216 |
|
else |
1217 |
|
{ |
1218 |
< |
if (right_splice_pos + seg_offset < _anchor_length) |
1219 |
< |
return false; |
1218 |
> |
if (right_splice_pos + seg_offset < _anchor_length) |
1219 |
> |
return false; |
1220 |
|
} |
1221 |
|
//uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos; |
1222 |
|
//atoi(toks[right_window_edge_field].c_str()); |
1223 |
|
int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1; |
1224 |
|
|
1225 |
|
string junction_strand = toks[num_extra_toks + strand_field]; |
1226 |
< |
if (!(junction_strand == "rev" || junction_strand == "fwd")|| |
1227 |
< |
!(orientation == '-' || orientation == '+')) |
1226 |
> |
if (!(junction_strand == "rev" || junction_strand == "fwd")) |
1227 |
> |
// || !(orientation == '-' || orientation == '+')) |
1228 |
|
{ |
1229 |
< |
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |
1230 |
< |
//fprintf(stderr, "junction_strand=%s, orientation='%c'\n", |
1231 |
< |
// junction_strand.c_str(), orientation); |
1232 |
< |
return false; |
1229 |
> |
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |
1230 |
> |
//fprintf(stderr, "junction_strand=%s, orientation='%c'\n", |
1231 |
> |
// junction_strand.c_str(), orientation); |
1232 |
> |
return false; |
1233 |
|
} |
1234 |
|
|
1267 |
– |
//vector<string> mismatch_toks; |
1268 |
– |
char* pch = strtok (mismatches,","); |
1235 |
|
int mismatches_in_anchor = 0; |
1236 |
< |
unsigned char num_mismatches = 0; |
1237 |
< |
while (pch != NULL) |
1238 |
< |
{ |
1239 |
< |
char* colon = strchr(pch, ':'); |
1240 |
< |
if (colon) |
1241 |
< |
{ |
1242 |
< |
*colon = 0; |
1243 |
< |
num_mismatches++; |
1278 |
< |
int mismatch_pos = atoi(pch); |
1279 |
< |
if ((orientation == '+' && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) || |
1280 |
< |
(orientation == '-' && abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len) |
1281 |
< |
mismatches_in_anchor++; |
1282 |
< |
} |
1283 |
< |
//mismatch_toks.push_back(pch); |
1284 |
< |
pch = strtok (NULL, ","); |
1285 |
< |
} |
1236 |
> |
for (int i=0;i<num_mismatches;i++) { |
1237 |
> |
spl_num_mismatches++; |
1238 |
> |
int mismatch_pos = mismatches[i]; |
1239 |
> |
if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) || |
1240 |
> |
((sam_flag & 0x0010) != 0 && |
1241 |
> |
abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len) |
1242 |
> |
mismatches_in_anchor++; |
1243 |
> |
} |
1244 |
|
|
1245 |
|
// FIXME: we probably should exclude these hits somewhere, but this |
1246 |
|
// isn't the right place |
1247 |
< |
vector<CigarOp> cigar; |
1248 |
< |
cigar.push_back(CigarOp(MATCH, left_splice_pos)); |
1247 |
> |
vector<CigarOp> splcigar; |
1248 |
> |
CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ; |
1249 |
> |
spliceCigar(splcigar, samcigar, left_splice_pos, gap_len, right_splice_pos, INS); |
1250 |
> |
/* |
1251 |
> |
splcigar.push_back(CigarOp(MATCH, left_splice_pos)); |
1252 |
|
if(toks[num_extra_toks + junction_type_field] == "del"){ |
1253 |
< |
cigar.push_back(CigarOp(DEL, gap_len)); |
1253 |
> |
splcigar.push_back(CigarOp(DEL, gap_len)); |
1254 |
|
}else{ |
1255 |
< |
cigar.push_back(CigarOp(REF_SKIP, gap_len)); |
1255 |
> |
splcigar.push_back(CigarOp(REF_SKIP, gap_len)); |
1256 |
|
} |
1257 |
< |
cigar.push_back(CigarOp(MATCH, right_splice_pos)); |
1258 |
< |
|
1257 |
> |
splcigar.push_back(CigarOp(MATCH, right_splice_pos)); |
1258 |
> |
*/ |
1259 |
|
bh = create_hit(name, |
1260 |
|
contig, |
1261 |
|
left, |
1262 |
< |
cigar, |
1263 |
< |
orientation == '-', |
1262 |
> |
splcigar, |
1263 |
> |
(sam_flag & 0x0010), |
1264 |
|
junction_strand == "rev", |
1265 |
< |
num_mismatches, |
1265 |
> |
spl_num_mismatches, |
1266 |
|
mismatches_in_anchor, |
1267 |
|
end); |
1268 |
|
return true; |
1269 |
|
} |
1270 |
< |
} |
1270 |
> |
} //parse splice data |
1271 |
|
else |
1272 |
|
{ |
1273 |
|
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |