We can detecte and open bam file

Author: natir

FastQt.pro

@@ -19,7 +19,7 @@ QMAKE_CXXFLAGS += -std=c++11
 unix {
 # COMPILE HTSLIB
 mytarget.target = $$PWD/htslib/libhts.so
-mytarget.commands = cd $$PWD/htslib; make -j4
+mytarget.commands = cd $$PWD/htslib; make -j 4
 mytarget_clean.commands = cd $$PWD/htslib; make clean
 QMAKE_EXTRA_TARGETS += mytarget
 PRE_TARGETDEPS += $$PWD/htslib/libhts.so

analysis/analysisrunner.cpp

@@ -64,7 +64,13 @@ void AnalysisRunner::run()
     QIODevice * file = Q_NULLPTR;
 
     file = new QFile(mFilename);
-    if (is_gz(file))
+    bool ubam = false;
+    if (is_ubam(file))
+    {
+        file = new QFile(mFilename);
+        ubam = true;
+    }
+    else if (is_gz(file))
     {
         file = new KCompressionDevice(mFilename, KCompressionDevice::GZip);
         if (!is_fastq(file))
@@ -105,26 +111,34 @@ void AnalysisRunner::run()
         mSequenceCount = 0;
         mProgression   = 0;
 
-        FastqReader reader(file);
+        AbstractSequenceReader* reader = nullptr;
+        if (ubam)
+        {
+            reader = new BamReader(static_cast<QFile*>(file));
+        }
+        else
+        {
+            reader = new FastqReader(file);
+        }
         mStartTime.start();
 
         // pre compute total size for sequencial access .
         //emitUpdate(tr("Analysis ..."));
-        reader.computeTotalSize();
+        reader->computeTotalSize();
 
 
 
         for (Analysis * a : mAnalysisHash)
             a->before();
 
 
-        while (reader.next() && !mCancel)
+        while (reader->next() && !mCancel)
         {
 
             // check if first sequence is valid..Means it's probably a good file
             if (mSequenceCount == 0)
             {
-                if (!reader.sequence().isValid())
+                if (!reader->sequence().isValid())
                 {
                     qCritical()<<Q_FUNC_INFO<<"Cannot read sequence. Are you sure it's a Fastq file ?";
                     setStatus(Canceled);
@@ -137,7 +151,7 @@ void AnalysisRunner::run()
             // this is critcal and can decrease the speed. Send message only 1 sequence / 1000
             if (mSequenceCount % 1000 == 0)
             {
-                int percentNow = reader.percentComplete();
+                int percentNow = reader->percentComplete();
                 // if percentNow is still null, return empty percent ...
                 if ( (percentNow >= mProgression + 5) || (percentNow == 0))
                 {
@@ -149,7 +163,7 @@ void AnalysisRunner::run()
 
             for (Analysis * a : mAnalysisHash)
             {
-                a->processSequence(reader.sequence());
+                a->processSequence(reader->sequence());
             }
         }
 
@@ -163,6 +177,7 @@ void AnalysisRunner::run()
 
         mDuration = mStartTime.elapsed();
 
+        delete reader;
     }
 
     else

main.cpp

@@ -1,3 +1,4 @@
+
 /*
 Copyright Copyright 2016-17 Sacha Schutz
 
@@ -29,7 +30,7 @@ Copyright Copyright 2016-17 Sacha Schutz
 #include "cliparser.h"
 #include "maincli.h"
 
-#include "bamreader.h"
+#include "htslib/hfile.h"
 
 int main(int argc, char *argv[])
 {

sequence/bamreader.h

@@ -7,7 +7,7 @@
 class BamReader : public AbstractSequenceReader
 {
 public:
-    BamReader(QFile * device);
+    BamReader(QFile *device);
     bool next() Q_DECL_OVERRIDE;
 
 

utils/format_detection.cpp

@@ -113,3 +113,34 @@ bool is_fastq(QIODevice* file)
         return false;
     }
 }
+
+bool is_ubam(QIODevice* file)
+{
+    /*Check is a gzip file*/
+    if(is_gz(file))
+    {
+        if (file->open(QIODevice::ReadOnly))
+        {
+            /*Check is a bgzp file*/
+            /*We use https://samtools.github.io/hts-specs/SAMv1.pdf page 10*/
+            file->seek(12);
+            QByteArray magic_number = file->read(3);
+            file->close();
+
+            if(magic_number.length() != 3)
+                return false;
+
+            return static_cast<unsigned char>(magic_number.at(0)) == 'B' \
+                   && static_cast<unsigned char>(magic_number.at(1)) == 'C' \
+                   && static_cast<unsigned char>(magic_number.at(2)) == 2;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}

utils/format_detection.h

@@ -32,4 +32,6 @@ bool is_xz(QIODevice* file);
 
 bool is_fastq(QIODevice* file);
 
+bool is_ubam(QIODevice* file);
+
 #endif // FORMAT_DETECTION_H