dnaa patch

gistfile1.diff

Subject: [PATCH] There may be cases in PE/MP data where the second read is not available.
 We have to consider those reads (singletons) as read1 not read2.

---
 dqc/dqc_postalignqc.c |   12 +++++++++++-
 1 files changed, 11 insertions(+), 1 deletions(-)

diff --git a/dqc/dqc_postalignqc.c b/dqc/dqc_postalignqc.c
index d3a4881..250b4d1 100644
--- a/dqc/dqc_postalignqc.c
+++ b/dqc/dqc_postalignqc.c
@@ -100,6 +100,7 @@ void dqc_postalignqc_core(dqc_postalignqc_opt *opt, char *bams[], int32_t n_bams
        dqc_postalignqc_bamstats_t *stats=NULL;
        counts_t *start_site_dist=NULL;
        RGBinary rg;
+  int32_t which; // read1 or read2 ?
 
        // Init
        RGBinaryReadBinary(&rg, NTSpace, opt->fasta_file_name); // reed in the reference
@@ -127,13 +128,22 @@ void dqc_postalignqc_core(dqc_postalignqc_opt *opt, char *bams[], int32_t n_bams
                                        if((b->core.flag & BAM_FREAD2)) {
                                                e->num_ends = 2; // always paired end
                                        }
+
+          // If this is a singleton, default to read1
+          if ((b->core.flag & BAM_FREAD1) == 0 && (b->core.flag & BAM_FREAD2) == 0) {
+            which = 1;
+          }
+          else {
+            which = (b->core.flag & BAM_FREAD1) ? 1 : 2;
+          }
+
                                        dqc_postalignqc_referrdist_update(e,
                                                        b,
                                                        &rg,
                                                        0, // gaps are not errors
                                                        1, // trim end gaps
                                                        (0 == bam_aux_getCSi(b, 0)) ? NTSpace : ColorSpace,
-                                                       (b->core.flag & BAM_FREAD1) ? 1 : 2);
+                                                       which);
                                }
                                // Insert size distribution
                                if(!(b->core.flag & BAM_FUNMAP) && // mapped
-- 
1.6.5.2