diff --git a/uk/ac/sanger/artemis/components/alignment/BamUtils.java b/uk/ac/sanger/artemis/components/alignment/BamUtils.java new file mode 100644 index 0000000000000000000000000000000000000000..b786157bca588d2225757a6e653d53c2e25363da --- /dev/null +++ b/uk/ac/sanger/artemis/components/alignment/BamUtils.java @@ -0,0 +1,217 @@ +/* BamUtils + * + * created: 2011 + * + * This file is part of Artemis + * + * Copyright(C) 2011 Genome Research Limited + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or(at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + */ +package uk.ac.sanger.artemis.components.alignment; + +import java.util.Hashtable; +import java.util.List; +import java.util.Vector; + +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.FeatureVector; +import uk.ac.sanger.artemis.components.FileViewer; + +class BamUtils +{ + /** + * + * @param features + * @param refName + * @param samFileReaderHash + * @param bamList + * @param seqNames + * @param offsetLengths + * @param concatSequences + * @param seqLengths + * @param samRecordFlagPredicate + * @param samRecordMapQPredicate + * @param contained + * @param useIntrons + */ + protected static void countReads(final FeatureVector features, + 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 SAMRecordFlagPredicate samRecordFlagPredicate, + final SAMRecordMapQPredicate samRecordMapQPredicate, + final boolean contained, + final boolean useIntrons) + { + Hashtable<String, List<Integer>> featureReadCount = new Hashtable<String, List<Integer>>(); + 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>(); + + for(String bam : bamList) + { + int cnt = 0; + if(!useIntrons && f.getSegments().size() > 1) + { + for(int j=0; j<f.getSegments().size(); j++) + { + start = f.getSegments().elementAt(j).getStart().getPosition(); + end = f.getSegments().elementAt(j).getEnd().getPosition(); + cnt += getCount(start, end, bam, refName, samFileReaderHash, + seqNames, offsetLengths, concatSequences, seqLengths, + samRecordFlagPredicate, samRecordMapQPredicate, contained); + } + } + else + cnt = getCount(start, end, bam, refName, samFileReaderHash, + seqNames, offsetLengths, concatSequences, seqLengths, + samRecordFlagPredicate, samRecordMapQPredicate, contained); + + sampleCounts.add(cnt); + } + + featureReadCount.put(f.getSystematicName(), sampleCounts); + } + + StringBuffer buff = new StringBuffer(); + buff.append("\t"); + for(String bam : bamList) + buff.append(bam+"\t"); + buff.append("\n"); + + for (String fId : featureReadCount.keySet() ) { + buff.append(fId+"\t"); + List<Integer> cnts = featureReadCount.get(fId); + for(int i=0; i<cnts.size(); i++) + buff.append(cnts.get(i)+(i<cnts.size()-1 ? "\t" : "")); + buff.append("\n"); + } + + FileViewer viewer = new FileViewer ("Read Count", true); + viewer.getTextPane().setText(buff.toString()); + } + + /** + * Count the reads in a range. + * @param start + * @param end + * @param bam + * @param refName + * @param samFileReaderHash + * @param seqNames + * @param offsetLengths + * @param concatSequences + * @param seqLengths + * @param samRecordFlagPredicate + * @param samRecordMapQPredicate + * @param contained + * @return + */ + private static int getCount( + final int start, + final int end, + final String bam, + final String refName, + final Hashtable<String, SAMFileReader> samFileReaderHash, + final Vector<String> seqNames, + final Hashtable<String, Integer> offsetLengths, + final boolean concatSequences, + final Hashtable<String, Integer> seqLengths, + final SAMRecordFlagPredicate samRecordFlagPredicate, + final SAMRecordMapQPredicate samRecordMapQPredicate, + final boolean contained) + { + int cnt = 0; + if(concatSequences) + { + int len = 0; + int lastLen = 1; + for(String name : seqNames) + { + int thisLength = seqLengths.get(name); + len += thisLength; + + if( (lastLen >= start && lastLen < end) || + (len >= start && len < end) || + (start >= lastLen && start < len) || + (end >= lastLen && end < len) ) + { + int offset = offsetLengths.get(name); + int thisStart = start - offset; + if(thisStart < 1) + thisStart = 1; + int thisEnd = end - offset; + if(thisEnd > thisLength) + thisEnd = thisLength; + + cnt = count(bam, samFileReaderHash, name, thisStart, thisEnd, + samRecordFlagPredicate, samRecordMapQPredicate, contained); + + } + lastLen = len; + } + } + else + { + cnt = count(bam, samFileReaderHash, refName, start, end, + samRecordFlagPredicate, samRecordMapQPredicate, contained); + } + return cnt; + } + + private static int count(String bam, + Hashtable<String, SAMFileReader> samFileReaderHash, + String refName, + int start, + int end, + SAMRecordFlagPredicate samRecordFlagPredicate, + SAMRecordMapQPredicate samRecordMapQPredicate, + boolean contained) + { + int cnt = 0; + SAMFileReader inputSam = samFileReaderHash.get(bam); + final CloseableIterator<SAMRecord> it = inputSam.query(refName, start, end, contained); + + while ( it.hasNext() ) + { + SAMRecord samRecord = it.next(); + if( samRecordFlagPredicate == null || + !samRecordFlagPredicate.testPredicate(samRecord)) + { + if(samRecordMapQPredicate == null || + samRecordMapQPredicate.testPredicate(samRecord)) + { + cnt++; + } + } + } + it.close(); + return cnt; + } +} + diff --git a/uk/ac/sanger/artemis/components/alignment/BamView.java b/uk/ac/sanger/artemis/components/alignment/BamView.java index 46557a4573c13db87a1df4893c0877bf10b1bc84..c60eef4f8c882932ae45ea128b78b7fe70ad3e1e 100644 --- a/uk/ac/sanger/artemis/components/alignment/BamView.java +++ b/uk/ac/sanger/artemis/components/alignment/BamView.java @@ -67,6 +67,7 @@ import java.util.List; import java.util.Vector; import javax.swing.BorderFactory; +import javax.swing.Box; import javax.swing.BoxLayout; import javax.swing.ButtonGroup; import javax.swing.JButton; @@ -1977,7 +1978,7 @@ public class BamView extends JPanel private void createMenus(JComponent menu) { - final JMenuItem addBam = new JMenuItem("Add Bam ..."); + final JMenuItem addBam = new JMenuItem("Add BAM ..."); menu.add(addBam); addBam.addActionListener(new ActionListener() { @@ -2000,6 +2001,32 @@ public class BamView extends JPanel bamFilesMenu.setFont(addBam.getFont()); menu.add(bamFilesMenu); + + final JMenu analyse = new JMenu("Analyse"); + menu.add(analyse); + final JMenuItem readCount = new JMenuItem("Read count for selected features ..."); + analyse.add(readCount); + readCount.addActionListener(new ActionListener() + { + public void actionPerformed(ActionEvent e) + { + 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); + + BamUtils.countReads(features, (String)combo.getSelectedItem(), samFileReaderHash, bamList, + seqNames, offsetLengths, concatSequences, seqLengths, + samRecordFlagPredicate, samRecordMapQPredicate, + !overlap.isSelected(), spliced.isSelected()); + } + }); + for(int i=0; i<bamList.size(); i++) addToViewMenu(i);