samValidator.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cstdlib>
   4 #include<cassert>
   5 #include<string>
   6 #include<set>
   7
   8 #include <stdint.h>
   9 #include "sam/bam.h"
  10 #include "sam/sam.h"
  11
  12 #include "my_assert.h"
  13
  14 using namespace std;
  15
  16 samfile_t *in;
  17 bam1_t *b, *b2;
  18
  19 set<string> used;
  20
  21 bool isValid;
  22
  23 int main(int argc, char* argv[]) {
  24         if (argc != 2) {
  25                 printf("Usage: rsem-sam-validator <input.sam/input.bam>\n");
  26                 exit(-1);
  27         }
  28
  29         string input_file(argv[1]);
  30         size_t pos = input_file.find_last_of('.');
  31         general_assert(pos != string::npos, "Input file does not have a suffix!");
  32         ++pos;
  33
  34         string suffix = input_file.substr(pos);
  35         size_t len = suffix.length();
  36         for (size_t i = 0; i < len; i++) suffix[i] = tolower(suffix[i]);
  37
  38         general_assert(suffix == "sam" || suffix == "bam", "Cannot recognize input file's file type! The file suffix is neither sam or bam.");
  39
  40         in = (suffix == "sam" ? samopen(argv[1], "r", NULL) : samopen(argv[1], "rb", NULL));
  41         general_assert(in != 0, "Cannot open input file!");
  42
  43         used.clear();
  44         b = bam_init1(); b2 = bam_init1();
  45
  46         isValid = true;
  47
  48         long cnt = 0;
  49         string cqname(""), qname;
  50         uint64_t creadlen = 0, readlen;
  51         bool cispaired = false, ispaired;
  52
  53         printf("."); fflush(stdout);
  54         do {
  55                 if (samread(in, b) < 0) break;
  56                 assert(b->core.l_qseq > 0);
  57
  58                 qname.assign(bam1_qname(b));
  59
  60                 // if this is a paired-end read
  61                 ispaired = b->core.flag & 0x0001;
  62                 if (ispaired) {
  63
  64                         isValid = (samread(in, b2) >= 0) && (qname == bam1_qname(b2)) && (b2->core.flag & 0x0001);
  65                         if (!isValid) { printf("\nOnly find one mate for paired-end read %s!\n", qname.c_str()); continue; }
  66                         assert(b2->core.l_qseq > 0);
  67                         isValid = !((b->core.flag & 0x0040) && (b->core.flag & 0x0080)) && !((b2->core.flag & 0x0040) & (b2->core.flag & 0x0080));
  68                         if (!isValid) { printf("\nRead %s has more than 2 segments (e.g. tripled or more ended reads)!\n", qname.c_str()); continue; }
  69                         isValid = !(((b->core.flag & 0x0040) && (b2->core.flag & 0x0040)) || ((b->core.flag & 0x0080) && (b2->core.flag & 0x0080)));
  70                         if (!isValid) { printf("\nThe two mates of paired-end read %s are not adjacent!\n", qname.c_str()); continue; }
  71
  72
  73                         // both mates are mapped
  74                         if (!(b->core.flag & 0x0004) && !(b2->core.flag & 0x0004)) {
  75                                 isValid = (b->core.tid == b2->core.tid) && (b->core.pos == b2->core.mpos) && (b2->core.pos == b->core.mpos);
  76                                 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; }
  77                         }
  78
  79                         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);
  80                 }
  81                 else readlen = b->core.l_qseq;
  82
  83                 if (cqname != qname) {
  84                         isValid = used.find(qname) == used.end();
  85                         if (!isValid) { printf("\nThe alignments of read %s are not grouped together!", qname.c_str()); continue; }
  86                         used.insert(cqname);
  87                         cqname = qname;
  88                         creadlen = readlen;
  89                         cispaired = ispaired;
  90                 }
  91                 else {
  92                         assert(cqname != "");
  93                         isValid = (creadlen == readlen);
  94                         if (!isValid) { printf("\nRead %s have different read/mate lengths!\n", qname.c_str()); continue; }
  95                         isValid = (cispaired == ispaired);
  96                         if (!isValid) { printf("\nRead %s is detected as both single-end read and paired-end read!\n", qname.c_str()); continue; }
  97                 }
  98
  99                 ++cnt;
 100                 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
 101
 102         } while(isValid);
 103
 104         bam_destroy1(b); bam_destroy1(b2);
 105         samclose(in);
 106
 107         if (isValid) printf("\nThe input file is valid!\n");
 108         else printf("The input file is not valid!\n");
 109
 110         return 0;
 111 }