- Add a function to calculate the value of the OA tag (#1354)

* - Add a function to calculate the value of the OA tag
samtools · Jun 17, 2019 · a61e233 · a61e233
1 parent f13d075
commit a61e233
Show file tree

Hide file tree

Showing 4 changed files with 82 additions and 6 deletions.
diff --git a/src/main/java/htsjdk/samtools/SAMRecordSetBuilder.java b/src/main/java/htsjdk/samtools/SAMRecordSetBuilder.java
@@ -231,16 +231,17 @@ private SAMRecord createReadNoFlag(final String name, final int contig, final in
             rec.setReferenceIndex(contig);
             rec.setAlignmentStart(start);
         }
-        if (!recordUnmapped) {
+        rec.setReadUnmappedFlag(recordUnmapped);
+        if (recordUnmapped) {
+            rec.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
+        } else {
             rec.setReadNegativeStrandFlag(negativeStrand);
             if (null != cigar) {
                 rec.setCigarString(cigar);
             } else if (!rec.getReadUnmappedFlag()) {
                 rec.setCigarString(readLength + "M");
             }
-            rec.setMappingQuality(255);
-        } else {
-            rec.setReadUnmappedFlag(true);
+            rec.setMappingQuality(SAMRecord.UNKNOWN_MAPPING_QUALITY);
         }
         rec.setAttribute(SAMTag.RG.name(), READ_GROUP_ID);
 
@@ -371,7 +372,7 @@ public void addPair(final String name, final int contig, final int start1, final
         if (useNmFlag) {
             end1.setAttribute(ReservedTagConstants.NM, 0);
         }
-        end1.setMappingQuality(255);
+        end1.setMappingQuality(SAMRecord.UNKNOWN_MAPPING_QUALITY);
         end1.setReadPairedFlag(true);
         end1.setProperPairFlag(true);
         end1.setFirstOfPairFlag(end1IsFirstOfPair);
@@ -393,7 +394,7 @@ public void addPair(final String name, final int contig, final int start1, final
         if (useNmFlag) {
             end2.setAttribute(ReservedTagConstants.NM, 0);
         }
-        end2.setMappingQuality(255);
+        end2.setMappingQuality(SAMRecord.UNKNOWN_MAPPING_QUALITY);
         end2.setReadPairedFlag(true);
         end2.setProperPairFlag(true);
         end2.setFirstOfPairFlag(!end1IsFirstOfPair);

diff --git a/src/main/java/htsjdk/samtools/SAMUtils.java b/src/main/java/htsjdk/samtools/SAMUtils.java
@@ -30,6 +30,7 @@
 import htsjdk.samtools.util.CoordMath;
 import htsjdk.samtools.util.RuntimeEOFException;
 import htsjdk.samtools.util.StringUtil;
+import htsjdk.tribble.annotation.Strand;
 
 import java.io.File;
 import java.io.IOException;
@@ -41,6 +42,7 @@
 import java.util.Collections;
 import java.util.List;
 import java.util.Map;
+import java.util.Optional;
 import java.util.TreeMap;
 import java.util.regex.Pattern;
 
@@ -1265,4 +1267,32 @@ public static List<SAMRecord> getOtherCanonicalAlignments(final SAMRecord record
     public static boolean isReferenceSequenceIncompatibleWithBAI(final SAMSequenceRecord sequence) {
         return sequence.getSequenceLength() > GenomicIndexUtil.BIN_GENOMIC_SPAN;
     }
+
+    /**
+     * Function to create the OA tag value from a record. The OA tag contains the mapping information
+     * of a record encoded as a comma-separated string (REF,POS,STRAND,CIGAR,MAPPING_QUALITY,NM_TAG_VALUE)
+     * @param record to use for generating the OA tag
+     * @return the OA tag string value
+     */
+    public static String calculateOATagValue(SAMRecord record) {
+        if (record.getReferenceName().contains(",")) {
+            throw new SAMException(String.format("Reference name for record %s contains a comma character.", record.getReadName()));
+        }
+        final String oaValue;
+        if (record.getReadUnmappedFlag()) {
+            oaValue = String.format("*,0,%s,*,%s,",
+                    record.getReadNegativeStrandFlag() ? Strand.NEGATIVE : Strand.POSITIVE,
+                    record.getMappingQuality());
+        } else {
+            oaValue = String.format("%s,%s,%s,%s,%s,%s",
+                    record.getReferenceName(),
+                    record.getAlignmentStart(),
+                    record.getReadNegativeStrandFlag() ? Strand.NEGATIVE : Strand.POSITIVE,
+                    record.getCigarString(),
+                    record.getMappingQuality(),
+                    Optional.ofNullable(record.getAttribute(SAMTag.NM.name())).orElse(""));
+        }
+        return oaValue;
+
+    }
 }
diff --git a/src/main/java/htsjdk/samtools/util/SequenceUtil.java b/src/main/java/htsjdk/samtools/util/SequenceUtil.java
@@ -25,6 +25,7 @@
 
 import htsjdk.samtools.*;
 import htsjdk.samtools.fastq.FastqConstants;
+import htsjdk.tribble.annotation.Strand;
 import htsjdk.utils.ValidationUtils;
 
 import java.io.File;

diff --git a/src/test/java/htsjdk/samtools/SAMUtilsTest.java b/src/test/java/htsjdk/samtools/SAMUtilsTest.java
@@ -28,7 +28,9 @@
 import org.testng.annotations.DataProvider;
 import org.testng.annotations.Test;
 
+import java.util.ArrayList;
 import java.util.Arrays;
+import java.util.Iterator;
 import java.util.List;
 
 public class SAMUtilsTest extends HtsjdkTest {
@@ -284,4 +286,46 @@ public void testCompressedBasesToBytes() {
                 'R', 'R', 'S', 'S', 'V', 'V', 'W', 'W', 'Y', 'Y', 'H', 'H', 'K', 'K', 'D', 'D', 'B', 'B'};
         Assert.assertEquals(new String(bytes), new String(expectedBases));
     }
+
+
+    @DataProvider()
+    public Iterator<Object[]> getOAValues(){
+        final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
+        final List<Object[]> tests = new ArrayList<>();
+
+        {
+            final SAMRecord record = builder.addFrag("test1", 0, 15, false, false, "36M", null, 45);
+
+            tests.add(new Object[]{record, "chr1,15,+,36M,255,"});
+        }
+        {
+            final SAMRecord record = builder.addFrag("test2", 0, 15, true, true, "36M", null, 45);
+            record.setMappingQuality(60);
+            // builder ignores this request....
+            record.setReadNegativeStrandFlag(true);
+            tests.add(new Object[]{record, "*,0,-,*,60,"});
+        }
+        {
+            final SAMRecord record = builder.addFrag("test3", 2, 15, false, false, "36M", null, 45);
+            record.setAttribute(SAMTag.NM.name(),33);
+
+            tests.add(new Object[]{record, "chr3,15,+,36M,255,33"});
+        }
+        {
+            final SAMRecord record = builder.addFrag("test4", 1, 115, true, false, "12S12M12I", null, 45);
+            record.setAttribute(SAMTag.NM.name(),0);
+            record.setMappingQuality(60);
+
+            tests.add(new Object[]{record, "chr2,115,-,12S12M12I,60,0"});
+        }
+
+        return tests.iterator();
+
+    }
+
+
+    @Test(dataProvider = "getOAValues")
+    public void testOAValues(final SAMRecord record, final String expectedOA){
+        Assert.assertEquals(SAMUtils.calculateOATagValue(record), expectedOA);
+    }
 }