Skip to content

Commit 9fe9042

Browse files
authored
SamUtils PairOrientation is confusingly ambiguous about paired reads needing to map to the same contig (#1709)
* - Improved documentation to clarify that PairOrientation (in SamUtils) pertains only to reads on the same contig. - Now getPairOrientation will throw an exception if thrown when the two reads are mapped to different contigs - Added tests to verify that the exception is thrown in the 4 cases specified.
1 parent a9aea4d commit 9fe9042

File tree

2 files changed

+80
-12
lines changed

2 files changed

+80
-12
lines changed

src/main/java/htsjdk/samtools/SamPairUtil.java

Lines changed: 25 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,8 @@ public class SamPairUtil {
4545
* FR means the read that's mapped to the forward strand comes before the
4646
* read mapped to the reverse strand when their 5'-end coordinates are
4747
* compared.
48+
*
49+
* PairOrientation only makes sense for a pair of reads that are both mapped to the same contig/chromosome
4850
*/
4951
public static enum PairOrientation
5052
{
@@ -57,31 +59,42 @@ public static enum PairOrientation
5759

5860
/**
5961
* Computes the pair orientation of the given SAMRecord.
60-
* @param r
62+
* @param record the record for which the pair orientation is requested
6163
* @return PairOrientation of the given SAMRecord.
6264
* @throws IllegalArgumentException If the record is not a paired read, or
63-
* one or both reads are unmapped.
65+
* one or both reads are unmapped, or if the two records are not mapped to the
66+
* same reference.
67+
*
68+
* NOTA BENE: This is NOT the orientation byte as used in Picard's MarkDuplicates. For that please look
69+
* in ReadEnds (in Picard).
6470
*/
65-
public static PairOrientation getPairOrientation(final SAMRecord r)
71+
public static PairOrientation getPairOrientation(final SAMRecord record)
6672
{
67-
final boolean readIsOnReverseStrand = r.getReadNegativeStrandFlag();
73+
final boolean readIsOnReverseStrand = record.getReadNegativeStrandFlag();
74+
75+
if(record.getReadUnmappedFlag() || !record.getReadPairedFlag() || record.getMateUnmappedFlag()) {
76+
throw new IllegalArgumentException("Invalid SAMRecord: " + record.getReadName() +
77+
". This method only works for SAMRecords that are paired reads with both reads aligned.");
78+
}
6879

69-
if(r.getReadUnmappedFlag() || !r.getReadPairedFlag() || r.getMateUnmappedFlag()) {
70-
throw new IllegalArgumentException("Invalid SAMRecord: " + r.getReadName() + ". This method only works for SAMRecords " +
71-
"that are paired reads with both reads aligned.");
80+
if (!record.getReferenceIndex().equals(record.getMateReferenceIndex())) {
81+
throw new IllegalArgumentException("Invalid SAMRecord: " + record.getReadName() +
82+
". This method only works for SAMRecords that are paired reads with both reads aligned to the" +
83+
" same reference. Found difference references:" + record.getReferenceName() + " and " +
84+
record.getMateReferenceName() + ".");
7285
}
7386

74-
if(readIsOnReverseStrand == r.getMateNegativeStrandFlag() ) {
87+
if(readIsOnReverseStrand == record.getMateNegativeStrandFlag() ) {
7588
return PairOrientation.TANDEM;
7689
}
7790

7891
final long positiveStrandFivePrimePos = ( readIsOnReverseStrand
79-
? r.getMateAlignmentStart() //mate's 5' position ( x---> )
80-
: r.getAlignmentStart() ); //read's 5' position ( x---> )
92+
? record.getMateAlignmentStart() //mate's 5' position ( x---> )
93+
: record.getAlignmentStart() ); //read's 5' position ( x---> )
8194

8295
final long negativeStrandFivePrimePos = ( readIsOnReverseStrand
83-
? r.getAlignmentEnd() //read's 5' position ( <---x )
84-
: r.getAlignmentStart() + r.getInferredInsertSize() ); //mate's 5' position ( <---x )
96+
? record.getAlignmentEnd() //read's 5' position ( <---x )
97+
: record.getAlignmentStart() + record.getInferredInsertSize() ); //mate's 5' position ( <---x )
8598

8699
return ( positiveStrandFivePrimePos < negativeStrandFivePrimePos
87100
? PairOrientation.FR

src/test/java/htsjdk/samtools/SamPairUtilTest.java

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030
import org.testng.annotations.Test;
3131

3232
import java.util.ArrayList;
33+
import java.util.Iterator;
3334
import java.util.List;
3435

3536

@@ -49,6 +50,12 @@ public void testGetPairOrientation(final String testName,
4950
Assert.assertEquals(SamPairUtil.getPairOrientation(rec2), expectedOrientation, testName + " second end");
5051
}
5152

53+
@Test(dataProvider = "testGetPairOrientationFail", expectedExceptions = IllegalArgumentException.class)
54+
public void testGetPairOrientationFail(final SAMRecord record) {
55+
SamPairUtil.getPairOrientation(record);
56+
}
57+
58+
5259
@Test(dataProvider = "testSetMateInfoMateCigar")
5360
public void testSetMateInfoMateCigar(final String testName,
5461
final int read1Start, final boolean read1Reverse, final String read1Cigar,
@@ -203,6 +210,54 @@ public Object[][] testGetPairOrientationDataProvider() {
203210
};
204211
}
205212

213+
@DataProvider(name = "testGetPairOrientationFail")
214+
public Iterator<Object[]> testGetPairOrientationFailDataProvider() {
215+
/**
216+
* @param testName
217+
* @param read1Start
218+
* @param read1Length
219+
* @param read1Reverse
220+
* @param read2Start
221+
* @param read2Length
222+
* @param read2Reverse
223+
*/
224+
225+
List<Object[]> tests=new ArrayList<>();
226+
final SAMFileHeader header = new SAMFileHeader();
227+
header.addSequence(new SAMSequenceRecord("chr1", 100000000));
228+
header.addSequence(new SAMSequenceRecord("chr2", 100000000));
229+
{
230+
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
231+
rec.setReadName("Unpaired");
232+
rec.setReadPairedFlag(false);
233+
tests.add(new Object[]{rec});
234+
}
235+
{
236+
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
237+
rec.setReadName("Unmapped");
238+
rec.setReadPairedFlag(true);
239+
rec.setReadUnmappedFlag(true);
240+
tests.add(new Object[]{rec});
241+
}
242+
{
243+
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
244+
rec.setReadName("Unmapped mate");
245+
rec.setReadPairedFlag(true);
246+
rec.setReferenceIndex(0);
247+
rec.setMateUnmappedFlag(true);
248+
tests.add(new Object[]{rec});
249+
}
250+
{
251+
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
252+
rec.setReadName("mate on difference reference");
253+
rec.setReadPairedFlag(true);
254+
rec.setReferenceIndex(0);
255+
rec.setMateReferenceIndex(1);
256+
tests.add(new Object[]{rec});
257+
}
258+
return tests.iterator();
259+
}
260+
206261
@DataProvider(name = "testSetMateInfoMateCigar")
207262
public Object[][] testSetMateInfoMateCigarDataProvider() {
208263
/**

0 commit comments

Comments
 (0)