Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

SamUtils PairOrientation is confusingly ambiguous about paired reads needing to map to the same contig #1709

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

SamUtils PairOrientation is confusingly ambiguous about paired reads needing to map to the same contig #1709

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 23 additions & 11 deletions src/main/java/htsjdk/samtools/SamPairUtil.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ public class SamPairUtil {
* FR means the read that's mapped to the forward strand comes before the
* read mapped to the reverse strand when their 5'-end coordinates are
* compared.
*
* PairOrientation only makes sense for a pair of reads that are both mapped to the same contig/chromosome
*/
public static enum PairOrientation
{
Expand All @@ -57,31 +59,41 @@ public static enum PairOrientation

/**
* Computes the pair orientation of the given SAMRecord.
* @param r
* @param record the record for which the pair orientation is requested
* @return PairOrientation of the given SAMRecord.
* @throws IllegalArgumentException If the record is not a paired read, or
* one or both reads are unmapped.
* one or both reads are unmapped, or is the two records do are not mapped to the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo here: "or is the two records do are not"

* same reference.
*
* NOTA BENE: This is NOT the orientation byte as used in Picard's MarkDuplicates. For that please look
* in ReadEnds (in Picard).
*/
public static PairOrientation getPairOrientation(final SAMRecord r)
public static PairOrientation getPairOrientation(final SAMRecord record)
{
final boolean readIsOnReverseStrand = r.getReadNegativeStrandFlag();
final boolean readIsOnReverseStrand = record.getReadNegativeStrandFlag();

if(r.getReadUnmappedFlag() || !r.getReadPairedFlag() || r.getMateUnmappedFlag()) {
throw new IllegalArgumentException("Invalid SAMRecord: " + r.getReadName() + ". This method only works for SAMRecords " +
if(record.getReadUnmappedFlag() || !record.getReadPairedFlag() || record.getMateUnmappedFlag()) {
throw new IllegalArgumentException("Invalid SAMRecord: " + record.getReadName() + ". This method only works for SAMRecords " +
"that are paired reads with both reads aligned.");
}

if(readIsOnReverseStrand == r.getMateNegativeStrandFlag() ) {
if (!record.getReferenceIndex().equals(record.getMateReferenceIndex())) {
throw new IllegalArgumentException("Invalid SAMRecord: " + record.getReadName() + ". This method only works for SAMRecords " +
"that are paired reads with both reads aligned to the same reference. Found difference references:" +
record.getReferenceName() + " and " + record.getMateReferenceName() + ".");
}

if(readIsOnReverseStrand == record.getMateNegativeStrandFlag() ) {
return PairOrientation.TANDEM;
}

final long positiveStrandFivePrimePos = ( readIsOnReverseStrand
? r.getMateAlignmentStart() //mate's 5' position ( x---> )
: r.getAlignmentStart() ); //read's 5' position ( x---> )
? record.getMateAlignmentStart() //mate's 5' position ( x---> )
: record.getAlignmentStart() ); //read's 5' position ( x---> )

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

return ( positiveStrandFivePrimePos < negativeStrandFivePrimePos
? PairOrientation.FR
Expand Down
55 changes: 55 additions & 0 deletions src/test/java/htsjdk/samtools/SamPairUtilTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
import org.testng.annotations.Test;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;


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

@Test(dataProvider = "testGetPairOrientationFail", expectedExceptions = IllegalArgumentException.class)
public void testGetPairOrientationFail(final SAMRecord record) {
SamPairUtil.getPairOrientation(record);
}


@Test(dataProvider = "testSetMateInfoMateCigar")
public void testSetMateInfoMateCigar(final String testName,
final int read1Start, final boolean read1Reverse, final String read1Cigar,
Expand Down Expand Up @@ -203,6 +210,54 @@ public Object[][] testGetPairOrientationDataProvider() {
};
}

@DataProvider(name = "testGetPairOrientationFail")
public Iterator<Object[]> testGetPairOrientationFailDataProvider() {
/**
* @param testName
* @param read1Start
* @param read1Length
* @param read1Reverse
* @param read2Start
* @param read2Length
* @param read2Reverse
*/

List<Object[]> tests=new ArrayList<>();
final SAMFileHeader header = new SAMFileHeader();
header.addSequence(new SAMSequenceRecord("chr1", 100000000));
header.addSequence(new SAMSequenceRecord("chr2", 100000000));
{
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
rec.setReadName("Unpaired");
rec.setReadPairedFlag(false);
tests.add(new Object[]{rec});
}
{
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
rec.setReadName("Unmapped");
rec.setReadPairedFlag(true);
rec.setReadUnmappedFlag(true);
tests.add(new Object[]{rec});
}
{
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
rec.setReadName("Unmapped mate");
rec.setReadPairedFlag(true);
rec.setReferenceIndex(0);
rec.setMateUnmappedFlag(true);
tests.add(new Object[]{rec});
}
{
final SAMRecord rec = makeSamRecord(header, 50, 50, false, true);
rec.setReadName("mate on difference reference");
rec.setReadPairedFlag(true);
rec.setReferenceIndex(0);
rec.setMateReferenceIndex(1);
tests.add(new Object[]{rec});
}
return tests.iterator();
}

@DataProvider(name = "testSetMateInfoMateCigar")
public Object[][] testSetMateInfoMateCigarDataProvider() {
/**
Expand Down
Loading