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 "utils.h"
  13 #include "my_assert.h"
  14
  15 using namespace std;
  16
  17 samfile_t *in;
  18 bam1_t *b, *b2;
  19
  20 set<string> used;
  21
  22 bool isValid;
  23
  24 int main(int argc, char* argv[]) {
  25         if (argc != 2) {
  26                 printf("Usage: rsem-sam-validator <input.sam/input.bam>\n");
  27                 exit(-1);
  28         }
  29
  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!");
  33         ++pos;
  34
  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]);
  38
  39         general_assert(suffix == "sam" || suffix == "bam", "Cannot recognize input file's file type! The file suffix is neither sam or bam.");
  40
  41         in = (suffix == "sam" ? samopen(argv[1], "r", NULL) : samopen(argv[1], "rb", NULL));
  42         general_assert(in != 0, "Cannot open input file!");
  43
  44         used.clear();
  45         b = bam_init1(); b2 = bam_init1();
  46
  47         isValid = true;
  48
  49         HIT_INT_TYPE cnt = 0;
  50         string cqname(""), qname;
  51         uint64_t creadlen = 0, readlen;
  52         bool cispaired = false, ispaired;
  53
  54         printf("."); fflush(stdout);
  55         do {
  56                 if (samread(in, b) < 0) break;
  57                 assert(b->core.l_qseq > 0);
  58
  59                 qname.assign(bam1_qname(b));
  60
  61                 // if this is a paired-end read
  62                 ispaired = b->core.flag & 0x0001;
  63                 if (ispaired) {
  64
  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; }
  72
  73
  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! If you're running covert-sam-for-rsem now, this might mean the read contains duplicate alignments.\n", qname.c_str()); continue; }
  78                         }
  79
  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);
  81                 }
  82                 else readlen = b->core.l_qseq;
  83
  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; }
  87                         used.insert(cqname);
  88                         cqname = qname;
  89                         creadlen = readlen;
  90                         cispaired = ispaired;
  91                 }
  92                 else {
  93                         assert(cqname != "");
  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; }
  98                 }
  99
 100                 ++cnt;
 101                 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
 102
 103         } while(isValid);
 104
 105         bam_destroy1(b); bam_destroy1(b2);
 106         samclose(in);
 107
 108         if (isValid) printf("\nThe input file is valid!\n");
 109         else printf("The input file is not valid!\n");
 110
 111         return 0;
 112 }