add option to rpkm calculation

git-svn-id: svn+ssh://svn.internal.sanger.ac.uk/repos/svn/pathsoft/artemis/trunk@15896 ee4ac58c-ac51-4696-9907-e4b3aa274f04

uk/ac/sanger/artemis/components/alignment/BamUtils.java

@@ -23,21 +23,30 @@
  */
 package uk.ac.sanger.artemis.components.alignment;
 
+import java.awt.Container;
+import java.awt.Frame;
+import java.text.DecimalFormat;
 import java.util.Hashtable;
 import java.util.List;
 import java.util.Vector;
 
+import javax.swing.JFrame;
+
 import net.sf.samtools.SAMFileReader;
 import net.sf.samtools.SAMRecord;
 import net.sf.samtools.util.CloseableIterator;
 import uk.ac.sanger.artemis.Feature;
+import uk.ac.sanger.artemis.FeatureSegmentVector;
 import uk.ac.sanger.artemis.FeatureVector;
+import uk.ac.sanger.artemis.components.EntryEdit;
 import uk.ac.sanger.artemis.components.FileViewer;
+import uk.ac.sanger.artemis.components.MultiComparator;
+import uk.ac.sanger.artemis.io.Range;
 
 class BamUtils
 {
   /**
-   * 
+   * Read count for selected reads.
    * @param features
    * @param refName
    * @param samFileReaderHash
@@ -50,6 +59,7 @@ class BamUtils
    * @param samRecordMapQPredicate
    * @param contained
    * @param useIntrons
+   * @param mappedReads
    */
   protected static void countReads(final FeatureVector features,
                                    final String refName,
@@ -62,26 +72,30 @@ class BamUtils
                                    final SAMRecordFlagPredicate samRecordFlagPredicate,
                                    final SAMRecordMapQPredicate samRecordMapQPredicate,
                                    final boolean contained,
-                                   final boolean useIntrons)
-  {    
-    Hashtable<String, List<Integer>> featureReadCount = new Hashtable<String, List<Integer>>();
+                                   final boolean useIntrons,
+                                   final int mappedReads[])
+  {
+    Hashtable<String, List<Float>> featureReadCount = new Hashtable<String, List<Float>>();
+
     for(int i=0; i<features.size(); i++)
     {     
       Feature f = features.elementAt(i);
       
-      int start = f.getFirstBase();
-      int end   = f.getLastBase();
-      List<Integer> sampleCounts = new Vector<Integer>();
+      int start  = f.getFirstBase();
+      int end    = f.getLastBase();
+      float fLenKb = getFeatureLength(f);
+      List<Float> sampleCounts = new Vector<Float>();
         
-      for(String bam : bamList)
+      for(int j=0; j<bamList.size(); j++)
       {
-        int cnt = 0;
+        String bam = bamList.get(j);
+        float cnt = 0;
         if(!useIntrons && f.getSegments().size() > 1)
         {
-          for(int j=0; j<f.getSegments().size(); j++)
+          for(int k=0; k<f.getSegments().size(); k++)
           {
-            start = f.getSegments().elementAt(j).getStart().getPosition();
-            end   = f.getSegments().elementAt(j).getEnd().getPosition();
+            start = f.getSegments().elementAt(k).getStart().getPosition();
+            end   = f.getSegments().elementAt(k).getEnd().getPosition();
             cnt += getCount(start, end, bam, refName, samFileReaderHash, 
               seqNames, offsetLengths, concatSequences, seqLengths, 
               samRecordFlagPredicate, samRecordMapQPredicate, contained); 
@@ -91,13 +105,16 @@ class BamUtils
           cnt = getCount(start, end, bam, refName, samFileReaderHash, 
               seqNames, offsetLengths, concatSequences, seqLengths, 
               samRecordFlagPredicate, samRecordMapQPredicate, contained); 
+       
+        if(mappedReads != null)
+          cnt = cnt / ((((float)mappedReads[j]) / 1000000.f)*fLenKb);
         
         sampleCounts.add(cnt);
       }
-
       featureReadCount.put(f.getSystematicName(), sampleCounts);
     }
 
+    DecimalFormat df = new DecimalFormat("0.00##");
     StringBuffer buff = new StringBuffer();
     for(String bam : bamList)
       buff.append("#BAM: "+bam+"\n");
@@ -105,16 +122,32 @@ class BamUtils
     
     for (String fId : featureReadCount.keySet() ) {
       buff.append(fId+"\t");
-      List<Integer> cnts = featureReadCount.get(fId);
+      List<Float> cnts = featureReadCount.get(fId);
       for(int i=0; i<cnts.size(); i++)
-        buff.append(cnts.get(i)+(i<cnts.size()-1 ? "\t" : ""));
+        buff.append(df.format(cnts.get(i)) + (i<cnts.size()-1 ? "\t" : ""));
       buff.append("\n");
     }
 
-    FileViewer viewer = new FileViewer ("Read Count", true, false, true);
+    FileViewer viewer;
+    if(mappedReads != null)
+      viewer = new FileViewer ("RPKM", true, false, true);
+    else
+      viewer = new FileViewer ("Read Count", true, false, true);
     viewer.getTextPane().setText(buff.toString());
   }
 
+  private static float getFeatureLength(Feature f)
+  {
+    FeatureSegmentVector segs = f.getSegments();
+    int len = 0;
+    for(int i=0; i<segs.size(); i++)
+    {
+      Range r = segs.elementAt(i).getRawRange();
+      len += r.getEnd()-r.getStart();
+    }
+    return ((float)len) / 1000.f;
+  }
+  
   /**
    * Count the reads in a range.
    * @param start
@@ -183,6 +216,81 @@ class BamUtils
     return cnt;
   }
 
+  /**
+   * Calculate the total number of mapped reads.
+   * @param refName
+   * @param samFileReaderHash
+   * @param bamList
+   * @param seqNames
+   * @param offsetLengths
+   * @param concatSequences
+   * @param seqLengths
+   * @param sequenceLength
+   * @return
+   */
+  protected static int[] getTotalMappedReads(
+      final String refName,
+      final Hashtable<String, SAMFileReader> samFileReaderHash,
+      final List<String> bamList, 
+      final Vector<String> seqNames,
+      final Hashtable<String, Integer> offsetLengths,
+      final boolean concatSequences, 
+      final Hashtable<String, Integer> seqLengths,
+      final int sequenceLength)
+  {
+    int MAX_BASE_CHUNK = 2000*60;
+    int mapped[] = new int[bamList.size()];
+    boolean contained = false;
+    SAMRecordFlagPredicate samRecordFlagPredicate = new SAMRecordFlagPredicate(SAMRecordFlagPredicate.READ_UNMAPPED_FLAG); 
+    SAMRecordMapQPredicate samRecordMapQPredicate = new SAMRecordMapQPredicate(-1);
+
+    for (int i = 0; i < sequenceLength; i += MAX_BASE_CHUNK)
+    {
+      int sbegc = i;
+      int sendc = i + MAX_BASE_CHUNK - 1;
+
+      for (int j=0; j<bamList.size(); j++)
+      {
+        String bam = bamList.get(j);
+        if (concatSequences)
+        {
+          int len = 0;
+          int lastLen = 1;
+          for (String name : seqNames)
+          {
+            int thisLength = seqLengths.get(name);
+            len += thisLength;
+
+            if ((lastLen >= sbegc && lastLen < sendc)
+                || (len >= sbegc && len < sendc)
+                || (sbegc >= lastLen && sbegc < len)
+                || (sendc >= lastLen && sendc < len))
+            {
+              int offset = offsetLengths.get(name);
+              int thisStart = sbegc - offset;
+              if (thisStart < 1)
+                thisStart = 1;
+              int thisEnd = sendc - offset;
+              if (thisEnd > thisLength)
+                thisEnd = thisLength;
+
+              mapped[j] += count(bam, samFileReaderHash, name, thisStart, thisEnd,
+                  samRecordFlagPredicate, samRecordMapQPredicate, contained);
+
+            }
+            lastLen = len;
+          }
+        }
+        else
+        {
+          mapped[j] += count(bam, samFileReaderHash, refName, sbegc, sendc,
+              samRecordFlagPredicate, samRecordMapQPredicate, contained);
+        }
+      }
+    }
+    return mapped;
+  }
+  
   private static int count(String bam, 
                     Hashtable<String, SAMFileReader> samFileReaderHash, 
                     String refName, 
@@ -212,5 +320,18 @@ class BamUtils
     it.close();
     return cnt;
   }
+  
+  protected static Container getBamContainer(BamView bamView)
+  {
+    Frame fs[] = JFrame.getFrames();
+    for(Frame f: fs)
+    {
+      if( f instanceof JFrame && 
+         ((JFrame)f) instanceof EntryEdit ||
+         ((JFrame)f) instanceof MultiComparator)
+        return ((JFrame)f).getContentPane();
+    }
+    return bamView;
+  }
 }
 

uk/ac/sanger/artemis/components/alignment/BamView.java

@@ -29,6 +29,8 @@ import java.awt.BorderLayout;
 import java.awt.Color;
 import java.awt.Component;
 import java.awt.Composite;
+import java.awt.Container;
+import java.awt.Cursor;
 import java.awt.Dimension;
 import java.awt.FlowLayout;
 import java.awt.FontMetrics;
@@ -2004,12 +2006,14 @@ public class BamView extends JPanel
     
     final JMenu analyse = new JMenu("Analyse");
     menu.add(analyse);
-    final JMenuItem readCount = new JMenuItem("Read count for selected features ...");
+    final JMenuItem readCount = new JMenuItem("Read count of selected features ...");
     analyse.add(readCount);
     readCount.addActionListener(new ActionListener()
     {
       public void actionPerformed(ActionEvent e)
       {
+        Container c = BamUtils.getBamContainer(BamView.this);
+        c.setCursor(new Cursor(Cursor.WAIT_CURSOR));
         FeatureVector features = feature_display.getSelection().getAllFeatures();
         
         JCheckBox overlap = new JCheckBox("Include all overlapping reads", true);
@@ -2023,7 +2027,46 @@ public class BamView extends JPanel
         BamUtils.countReads(features, (String)combo.getSelectedItem(), samFileReaderHash, bamList,
             seqNames, offsetLengths, concatSequences, seqLengths, 
             samRecordFlagPredicate, samRecordMapQPredicate,
-            !overlap.isSelected(), spliced.isSelected());
+            !overlap.isSelected(), spliced.isSelected(), null);
+        c.setCursor(new Cursor(Cursor.DEFAULT_CURSOR));
+      } 
+    });
+    
+    final JMenuItem rpkm = new JMenuItem("RPKM value of selected features ...");
+    analyse.add(rpkm);
+    rpkm.addActionListener(new ActionListener()
+    {
+      public void actionPerformed(ActionEvent e)
+      {
+        Container c = BamUtils.getBamContainer(BamView.this);
+        c.setCursor(new Cursor(Cursor.WAIT_CURSOR));
+        FeatureVector features = feature_display.getSelection().getAllFeatures();
+        
+        JCheckBox overlap = new JCheckBox("Include all overlapping reads", true);
+        overlap.setToolTipText("Include reads that partially overlap the feature");
+        JCheckBox spliced = new JCheckBox("Introns included", true);
+        Box yBox = Box.createVerticalBox();
+        yBox.add(overlap);
+        yBox.add(spliced);
+        JOptionPane.showMessageDialog(null, yBox, "Read Count Option", JOptionPane.INFORMATION_MESSAGE);
+        
+        int seqlen = 0;
+        if(feature_display != null)
+          seqlen = feature_display.getSequenceLength();
+        else if(bases != null)
+          seqlen = bases.getLength();
+
+        int mappedReads[] =
+          BamUtils.getTotalMappedReads((String)combo.getSelectedItem(),
+              samFileReaderHash, bamList, seqNames, offsetLengths, concatSequences, 
+              offsetLengths, seqlen);
+
+        logger4j.debug("TOTAL MAPPED READS "+mappedReads);
+        BamUtils.countReads(features, (String)combo.getSelectedItem(), samFileReaderHash, bamList,
+            seqNames, offsetLengths, concatSequences, seqLengths, 
+            samRecordFlagPredicate, samRecordMapQPredicate,
+            !overlap.isSelected(), spliced.isSelected(), mappedReads);
+        c.setCursor(new Cursor(Cursor.DEFAULT_CURSOR));
       } 
     });