13 #include "my_assert.h"
24 int main(int argc, char* argv[]) {
26 printf("Usage: rsem-sam-validator <input.sam/input.bam>\n");
30 string input_file(argv[1]);
31 size_t pos = input_file.find_last_of('.');
32 general_assert(pos != string::npos, "Input file does not have a suffix!");
35 string suffix = input_file.substr(pos);
36 size_t len = suffix.length();
37 for (size_t i = 0; i < len; i++) suffix[i] = tolower(suffix[i]);
39 general_assert(suffix == "sam" || suffix == "bam", "Cannot recognize input file's file type! The file suffix is neither sam or bam.");
41 in = (suffix == "sam" ? samopen(argv[1], "r", NULL) : samopen(argv[1], "rb", NULL));
42 general_assert(in != 0, "Cannot open input file!");
45 b = bam_init1(); b2 = bam_init1();
50 string cqname(""), qname;
51 uint64_t creadlen = 0, readlen;
52 bool cispaired = false, ispaired;
54 printf("."); fflush(stdout);
56 if (samread(in, b) < 0) break;
57 assert(b->core.l_qseq > 0);
59 qname.assign(bam1_qname(b));
61 // if this is a paired-end read
62 ispaired = b->core.flag & 0x0001;
65 isValid = (samread(in, b2) >= 0) && (qname == bam1_qname(b2)) && (b2->core.flag & 0x0001);
66 if (!isValid) { printf("\nOnly find one mate for paired-end read %s!\n", qname.c_str()); continue; }
67 assert(b2->core.l_qseq > 0);
68 isValid = !((b->core.flag & 0x0040) && (b->core.flag & 0x0080)) && !((b2->core.flag & 0x0040) & (b2->core.flag & 0x0080));
69 if (!isValid) { printf("\nRead %s has more than 2 segments (e.g. tripled or more ended reads)!\n", qname.c_str()); continue; }
70 isValid = !(((b->core.flag & 0x0040) && (b2->core.flag & 0x0040)) || ((b->core.flag & 0x0080) && (b2->core.flag & 0x0080)));
71 if (!isValid) { printf("\nThe two mates of paired-end read %s are not adjacent!\n", qname.c_str()); continue; }
74 // both mates are mapped
75 if (!(b->core.flag & 0x0004) && !(b2->core.flag & 0x0004)) {
76 isValid = (b->core.tid == b2->core.tid) && (b->core.pos == b2->core.mpos) && (b2->core.pos == b->core.mpos);
77 if (!isValid) { printf("\nOne of paired-end read %s's alignment does not have two mates adjacent to each other!\n", qname.c_str()); continue; }
80 readlen = ((b->core.flag & 0x0040) ? (uint64_t(b->core.l_qseq) << 32) + b2->core.l_qseq : (uint64_t(b2->core.l_qseq) << 32) + b->core.l_qseq);
82 else readlen = b->core.l_qseq;
84 if (cqname != qname) {
85 isValid = used.find(qname) == used.end();
86 if (!isValid) { printf("\nThe alignments of read %s are not grouped together!", qname.c_str()); continue; }
94 isValid = (creadlen == readlen);
95 if (!isValid) { printf("\nRead %s have different read/mate lengths!\n", qname.c_str()); continue; }
96 isValid = (cispaired == ispaired);
97 if (!isValid) { printf("\nRead %s is detected as both single-end read and paired-end read!\n", qname.c_str()); continue; }
101 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
105 bam_destroy1(b); bam_destroy1(b2);
108 if (isValid) printf("\nThe input file is valid!\n");
109 else printf("The input file is not valid!\n");