JAL-2909 separate reverse and forward strand reads. fix subtle bug in handling insert...

[jalview.git] / src / jalview / io / BamFile.java

diff --git a/src/jalview/io/BamFile.java b/src/jalview/io/BamFile.java
index 8a04867..8d55c53 100644 (file)
--- a/src/jalview/io/BamFile.java
+++ b/src/jalview/io/BamFile.java
@@ -125,6 +125,44 @@ public class BamFile extends AlignFile
     return null;
   }
 
+  private StringBuilder insertRefSeq(
+          SortedMap<Integer, Integer> insertions, Range xtent)
+  {
+    StringBuilder refseq = new StringBuilder();
+    int inserted = 0;
+    for (int p = xtent.start; p < xtent.end; p++)
+    {
+      refseq.append('N');
+    }
+    for (Map.Entry<Integer, Integer> insert : insertions.entrySet())
+    {
+      int inspos = insert.getKey() - xtent.start + inserted;
+      if (inspos < 0)
+      {
+        System.out.println("Ignoring -ve insert position " + insert.getKey()
+                + " of " + insert.getValue() +
+                " (alpos: " + inspos + ")");
+        inspos = 0;
+        // continue;
+      }
+      for (int i = 0, j = insert.getValue(); i < j; i++)
+      {
+        inserted++;
+        refseq.insert(inspos, '-');
+      }
+    }
+    return refseq;
+  }
+
+  private void padRef(SequenceI ref, CigarParser cig)
+  { // pad to column
+    int padding = cig.firstAlColumn - ref.findIndex(cig.firstRposCol);
+    System.out.println("Padding " + padding + " to move refseq position "
+            + cig.firstRposCol + " to " + cig.firstAlColumn + "col.");
+
+    ref.insertCharAt(0, padding, '-');
+  }
+
   @Override
   public void parse()
   {
@@ -134,10 +172,30 @@ public class BamFile extends AlignFile
       SAMRecordIterator it = fileReader.query(chromosome, start, end,
               false);
       CigarParser parser = new CigarParser('-');
-      SortedMap<Integer, Integer> insertions = parser.getInsertions(it);
+      Range[] xtent = new Range[] { new Range(start, end) };
+      SortedMap<Integer, Integer> insertions[] = parser.getInsertions(it,
+              xtent);
       it.close();
 
+      SequenceI refSeq = new Sequence("chr:" + chromosome,
+              insertRefSeq(insertions[0], xtent[0])
+                      .toString());
+
+      refSeq.setStart(xtent[0].start);
+      refSeq.setEnd(xtent[0].end);
+      SequenceI revRefSeq = new Sequence("rev:chr:" + chromosome,
+              insertRefSeq(insertions[1], xtent[0])
+                      .toString());
+      revRefSeq.setStart(xtent[0].start);
+      revRefSeq.setEnd(xtent[0].end);
+
+      // Hack to move the seqs along
+      padRef(refSeq, parser);
+      padRef(revRefSeq, parser);
+
       it = fileReader.query(chromosome, start, end, false);
+
+      ArrayList<SequenceI> fwd = new ArrayList(), rev = new ArrayList();
       while (it.hasNext())
       {
         SAMRecord rec = it.next();
@@ -150,7 +208,9 @@ public class BamFile extends AlignFile
         }
 
         // make dataset sequence: start at 1, end at read length
-        SequenceI seq = new Sequence(rec.getReadName(),
+        SequenceI seq = new Sequence(
+                "" + (rec.getReadNegativeStrandFlag() ? "rev:" : "")
+                        + rec.getReadName(),
                 rec.getReadString().toLowerCase());
         seq.setStart(1);
         seq.setEnd(rec.getReadLength());
@@ -163,7 +223,8 @@ public class BamFile extends AlignFile
                   (float) q.next() - ' ', "bamfile"));
           p++;
         }
-        String newRead = parser.parseCigarToSequence(rec, insertions,
+        String newRead = parser.parseCigarToSequence(rec,
+                insertions[rec.getReadNegativeStrandFlag() ? 1 : 0],
                 alignmentStart, seq);
 
         // make alignment sequences
@@ -172,7 +233,28 @@ public class BamFile extends AlignFile
 
         // set start relative to soft clip; assume end is set by Sequence code
         alsq.setStart(rec.getStart() - rec.getUnclippedStart() + 1);
-        seqs.add(alsq);
+        if (rec.getReadNegativeStrandFlag())
+        {
+          rev.add(alsq);
+        }
+        else
+        {
+          fwd.add(alsq);
+        }
+      }
+      // Add forward
+      if (fwd.size() > 0)
+      {
+        seqs.add(refSeq); // FIXME needs to be optional, and properly padded
+        seqs.addAll(fwd);
+        fwd.clear();
+      }
+      // and reverse strand reads.
+      if (rev.size() > 0)
+      {
+        seqs.add(revRefSeq); // FIXME needs to be optional and properly padded
+        seqs.addAll(rev);
+        rev.clear();
       }
     }
   }