Java Code Examples for htsjdk.samtools.SAMRecord#getReadNegativeStrandFlag()

/**
 * Checks the first ADAPTER_MATCH_LENGTH bases of the read against known adapter sequences and returns
 * true if the read matches an adapter sequence with MAX_ADAPTER_ERRORS mismsatches or fewer.
 * <p>
 * Only unmapped reads and reads with MQ=0 are considers eligible for being adapter
 */

 public boolean isAdapter(final SAMRecord record) {

    // If the read mapped with mapping quality > 0 we do not consider it an adapter read.
    if (!record.getReadUnmappedFlag() && record.getMappingQuality() != 0) {
        return false;
    }

    // if the read is mapped and is aligned to the reverse strand, it needs to be
    // rev-comp'ed before calling isAdapterSequence.
    final boolean revComp = !record.getReadUnmappedFlag() &&
            record.getReadNegativeStrandFlag();

    final byte[] readBases = record.getReadBases();
    return isAdapterSequence(readBases, revComp);
}
 

@Test
// One CODING read (wrong strand, no other gene models overlapping)
public void testCodingReadWrongStrand () {
	SAMRecord r = getFakeRecord(testBAMFile, 2, 8, true);
	boolean negStrandFlag = r.getReadNegativeStrandFlag();
	// gene with 2 exons, 1 coding from 1-10, one UTR from 91-100.  Positive strand gene.
	GeneFromGTF gene = new GeneFromGTF(r.getContig(), 1, 100, false, "A", "coding", "A", "coding", 1);
	final GeneFromGTF.TranscriptFromGTF tx = gene.addTranscript("trans1", 1, 100, 1, 90, 2, "trans1", "trans1", "coding");
	tx.addExon(1, 10);
	tx.addExon(91, 100);
	OverlapDetector<Gene> geneOverlapDetector = new OverlapDetector<>(0, 0);
	geneOverlapDetector.addLhs(gene, gene);
	TagReadWithGeneFunction tagger = new TagReadWithGeneFunction();
	r=tagger.setAnnotations(r, geneOverlapDetector, false);
	Assert.assertEquals(r.getStringAttribute("gn"), gene.getName());
	Assert.assertEquals(r.getStringAttribute("gs"), "+");
	Assert.assertEquals(r.getStringAttribute("gf"), LocusFunction.CODING.name());
	Assert.assertEquals(r.getStringAttribute("XF"), LocusFunction.INTERGENIC.name());
}
 

private void addRead(final GcObject gcObj, final SAMRecord rec, final String group, final byte[] gc, final byte[] refBases) {
    if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) ++gcObj.totalClusters;
    final int pos = rec.getReadNegativeStrandFlag() ? rec.getAlignmentEnd() - scanWindowSize : rec.getAlignmentStart();
    ++gcObj.totalAlignedReads;
    if (pos > 0) {
        final int windowGc = gc[pos];
        if (windowGc >= 0) {
            ++gcObj.readsByGc[windowGc];
            gcObj.basesByGc[windowGc] += rec.getReadLength();
            gcObj.errorsByGc[windowGc] +=
                    SequenceUtil.countMismatches(rec, refBases, bisulfite) +
                            SequenceUtil.countInsertedBases(rec) + SequenceUtil.countDeletedBases(rec);
        }
    }
    if (gcObj.group == null) {
        gcObj.group = group;
    }
}
 

public void incrementCount(SAMRecord read) {
	if (!readIds.contains(read.getReadName())) {
		// Don't allow multiple ends of fragment to be double counted
		count += 1;
	}
	
	totalCount += 1;
	
	if (read.getReadNegativeStrandFlag()) {
		incrementRev();
	} else {
		incrementFwd();
	}
	
	readIds.add(read.getReadName());
	if (!alignmentEnds.containsKey(read.getAlignmentEnd())) {
		alignmentEnds.put(read.getAlignmentEnd(), new ArrayList<String>());
	}
	
	alignmentEnds.get(read.getAlignmentEnd()).add(read.getReadName());
}
 

public Map<Gene, LocusFunction> getLocusFunctionForReadByGene (final SAMRecord rec, final OverlapDetector<Gene> geneOverlapDetector) {
	Map<Gene, LocusFunction> result = new HashMap<>();
	final Interval readInterval = new Interval(rec.getReferenceName(), rec.getAlignmentStart(), rec.getAlignmentEnd(), rec.getReadNegativeStrandFlag(), rec.getReadName());
	final Collection<Gene> overlappingGenes = geneOverlapDetector.getOverlaps(readInterval);

       for (Gene g: overlappingGenes) {
       	LocusFunction f = getLocusFunctionForRead(rec, g);
           result.put(g, f);
       }
       return result;
}
 

private void countRecord(final SAMRecord rec) {
    final boolean isNegativeStrand = rec.getReadNegativeStrandFlag();
    final boolean isReadTwo = rec.getReadPairedFlag() && rec.getSecondOfPairFlag();
    if (isReadTwo) {
        if (isNegativeStrand) this.R2_NEG++;
        else this.R2_POS++;
    } else {
        if (isNegativeStrand) this.R1_NEG++;
        else this.R1_POS++;
    }
}
 

/**
 * Returns 1-based index of first base in read that corresponds to M in CIGAR string.
 * Note that first is relative to 5' end, so that for reverse-strand alignment, the index of
 * the last base aligned is computed relative to the end of the read.
 */
int getIndexOfFirstAlignedBase(final SAMRecord rec) {
    final List<AlignmentBlock> alignmentBlocks = rec.getAlignmentBlocks();
    if (rec.getReadNegativeStrandFlag()) {
        final AlignmentBlock alignmentBlock = alignmentBlocks.get(alignmentBlocks.size() - 1);
        return rec.getReadLength() - CoordMath.getEnd(alignmentBlock.getReadStart(), alignmentBlock.getLength()) + 1;
    } else {
        return alignmentBlocks.get(0).getReadStart();
    }
}
 

public MateKey getOriginalReadInfo(SAMRecord read) {
		int pos = read.getAlignmentStart();
		boolean isUnmapped = read.getReadUnmappedFlag();
		boolean isRc = read.getReadNegativeStrandFlag();
		
		String yo = read.getStringAttribute("YO");
		
		if (yo != null) {
			if (yo.startsWith("N/A")) {
				// Original alignment was unmapped
				isUnmapped = true;
//				isRc = false;
				// Orientation is forced to be opposite of mate during realignment
				// regardless of the original alignment.
				isRc = !read.getMateNegativeStrandFlag();
			} else {
				String[] fields = yo.split(":");
				pos = Integer.parseInt(fields[1]);
				isUnmapped = false;
				isRc = fields[2].equals("-") ? true : false;
			}
		}
		
		int readNum = read.getFirstOfPairFlag() ? 1 : 2;
		
		return new MateKey(read.getReadName(), pos, isUnmapped, isRc, readNum, read.getAlignmentStart());
	}
 

@Override
public int compare( final SAMRecord samRecord1, final SAMRecord samRecord2 ) {
    int cmp = compareCoordinates(samRecord1, samRecord2);
    if ( cmp != 0 ) {
        return cmp;
    }

    // Test of negative strand flag is not really necessary, because it is tested
    // via getFlags(), but it is left here because that is the way it is done
    // in {@link htsjdk.samtools.SAMRecordCoordinateComparator}
    if ( samRecord1.getReadNegativeStrandFlag() == samRecord2.getReadNegativeStrandFlag() ) {
        cmp = samRecord1.getReadName().compareTo(samRecord2.getReadName());
        if ( cmp != 0 ) {
            return cmp;
        }
        cmp = Integer.compare(samRecord1.getFlags(), samRecord2.getFlags());
        if ( cmp != 0 ) {
            return cmp;
        }
        cmp = Integer.compare(samRecord1.getMappingQuality(), samRecord2.getMappingQuality());
        if ( cmp != 0 ) {
            return cmp;
        }
        cmp = Integer.compare(header.getSequenceIndex(samRecord1.getMateReferenceName()), header.getSequenceIndex(samRecord2.getMateReferenceName()));
        if ( cmp != 0 ) {
            return cmp;
        }
        cmp = Integer.compare(samRecord1.getMateAlignmentStart(), samRecord2.getMateAlignmentStart());
        if ( cmp != 0 ) {
            return cmp;
        }
        cmp = Integer.compare(samRecord1.getInferredInsertSize(), samRecord2.getInferredInsertSize());
        return cmp;
    }
    else {
        return samRecord1.getReadNegativeStrandFlag() ? 1: -1;
    }
}
 

static int getCoordinate (final SAMRecord rec) {
	int coordinate=-1;
	if (rec.getReadNegativeStrandFlag())
		coordinate=rec.getUnclippedEnd();
	else
		coordinate=rec.getUnclippedStart();
	if (coordinate <1) coordinate=1;  // thanks for the negative coordinates, makes sense.
	return coordinate;
}
 

/**
 * Minimal test of merging data from separate read 1 and read 2 alignments
 */
@Test(dataProvider="separateTrimmed")
public void testMergingFromSeparatedReadTrimmedAlignments(final File unmapped, final List<File> r1Align, final List<File> r2Align, final int r1Trim, final int r2Trim, final String testName) throws Exception {
     final File output = File.createTempFile("mergeMultipleAlignmentsTest", ".sam");
     output.deleteOnExit();

    doMergeAlignment(unmapped, null, r1Align, r2Align, r1Trim, r2Trim,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReaderFactory factory = SamReaderFactory.makeDefault();
    final SamReader result = factory.open(output);
     for (final SAMRecord sam : result) {
        // Get the alignment record
        final List<File> rFiles = sam.getFirstOfPairFlag() ? r1Align : r2Align;
        SAMRecord alignment = null;
        for (final File f : rFiles) {
            for (final SAMRecord tmp : factory.open(f)) {
                if (tmp.getReadName().equals(sam.getReadName())) {
                    alignment = tmp;
                    break;
                }
            }
            if (alignment != null) break;
        }

        final int trim = sam.getFirstOfPairFlag() ? r1Trim : r2Trim;
        final int notWrittenToFastq = sam.getReadLength() - (trim + alignment.getReadLength());
        final int beginning = sam.getReadNegativeStrandFlag() ? notWrittenToFastq : trim;
        final int end = sam.getReadNegativeStrandFlag() ? trim : notWrittenToFastq;

        if (!sam.getReadUnmappedFlag()) {
            final CigarElement firstMergedCigarElement = sam.getCigar().getCigarElement(0);
            final CigarElement lastMergedCigarElement = sam.getCigar().getCigarElement(sam.getCigar().getCigarElements().size() - 1);
            final CigarElement firstAlignedCigarElement = alignment.getCigar().getCigarElement(0);
            final CigarElement lastAlignedCigarElement = alignment.getCigar().getCigarElement(alignment.getCigar().getCigarElements().size() - 1);

            if (beginning > 0) {
                Assert.assertEquals(firstMergedCigarElement.getOperator(), CigarOperator.S, "First element is not a soft clip");
                Assert.assertEquals(firstMergedCigarElement.getLength(), beginning + ((firstAlignedCigarElement.getOperator() == CigarOperator.S) ? firstAlignedCigarElement.getLength() : 0));
            }
            if (end > 0) {
                Assert.assertEquals(lastMergedCigarElement.getOperator(), CigarOperator.S, "Last element is not a soft clip");
                Assert.assertEquals(lastMergedCigarElement.getLength(), end + ((lastAlignedCigarElement.getOperator() == CigarOperator.S) ? lastAlignedCigarElement.getLength() : 0));
            }
        }
    }

}
 

public boolean checkDuplicates(final SAMRecord record)
{
    if(record.getDuplicateReadFlag())
        return true;

    if(!mMarkDuplicates)
        return false;

    if(mDuplicateReadIds.contains(record.getReadName()))
    {
        mDuplicateReadIds.remove(record.getReadName());
        return true;
    }

    if(!record.getReferenceName().equals(record.getMateReferenceName()) || record.getReadNegativeStrandFlag() == record.getMateNegativeStrandFlag())
        return false;

    int firstStartPos = record.getFirstOfPairFlag() ? record.getStart() : record.getMateAlignmentStart();
    int secondStartPos = record.getFirstOfPairFlag() ? record.getMateAlignmentStart() : record.getStart();
    int readLength = record.getReadLength();
    int insertSize = record.getInferredInsertSize();

    List<int[]> dupDataList = mDuplicateCache.get(firstStartPos);

    if(dupDataList == null)
    {
        dupDataList = Lists.newArrayList();
        mDuplicateCache.put(firstStartPos, dupDataList);
    }
    else
    {
        // search for a match
        if(dupDataList.stream().anyMatch(x -> x[DUP_DATA_SECOND_START] == secondStartPos
                && x[DUP_DATA_READ_LEN] == readLength && insertSize == x[DUP_DATA_INSERT_SIZE]))
        {
            ISF_LOGGER.trace("duplicate fragment: id({}) chr({}) pos({}->{}) otherReadStart({}) insertSize({})",
                    record.getReadName(), record.getReferenceName(), firstStartPos, record.getEnd(), secondStartPos, insertSize);

            // cache so the second read can be identified immediately
            mDuplicateReadIds.add(record.getReadName());
            return true;
        }
    }

    int[] dupData = {secondStartPos, readLength, insertSize};
    dupDataList.add(dupData);

    return false;
}
 

public static int mergeReadPair(SAMRecordWrapper readWrapper, Map<String, SAMRecordWrapper> firstReads, Map<String, SAMRecordWrapper> secondReads) {
	
	int alignmentStart = -1;
	SAMRecord read = readWrapper.getSamRecord();
	
	if (read.getReadPairedFlag() && !read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getReadNegativeStrandFlag() != read.getMateNegativeStrandFlag()) {
		
		SAMRecordWrapper pair = null;
		
		if (read.getFirstOfPairFlag()) {
			pair = secondReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		} else {
			pair = firstReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		}
		
		if (pair != null) {
			SAMRecordWrapper first = null;
			SAMRecordWrapper second = null;
			if (read.getReadNegativeStrandFlag()) {
				first = pair;
				second = readWrapper;
			} else {
				first = readWrapper;
				second = pair;
			}

			
			if (first.getAdjustedAlignmentStart() < second.getAdjustedAlignmentStart() &&
					first.getAdjustedAlignmentEnd() > second.getAdjustedAlignmentStart() &&
					first.getSamRecord().getReadLength() > MIN_OVERLAP && 
					second.getSamRecord().getReadLength() > MIN_OVERLAP) {

				Pair<String, String> merged = mergeSequences(first.getSamRecord().getReadString(), second.getSamRecord().getReadString(),
						first.getSamRecord().getBaseQualityString(), second.getSamRecord().getBaseQualityString());
				
				if (merged != null) {
					readWrapper.setMerged(merged.getFirst(), merged.getSecond(), first.getAdjustedAlignmentStart(), second.getAdjustedAlignmentEnd());
					pair.setMerged(merged.getFirst(), merged.getSecond(), first.getAdjustedAlignmentStart(), second.getAdjustedAlignmentEnd());
					
					alignmentStart = first.getAdjustedAlignmentStart();
				}
			}
		}
	}
	
	return alignmentStart;
}
 

public static boolean hasPossibleAdapterReadThrough(SAMRecord read, Map<String, SAMRecordWrapper> firstReads, Map<String, SAMRecordWrapper> secondReads) {
	
	boolean hasReadThrough = false;
	
	// Check for fragment read through in paired end
	if (read.getReadPairedFlag() && !read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getAlignmentStart() == read.getMateAlignmentStart() &&
			read.getReadNegativeStrandFlag() != read.getMateNegativeStrandFlag()) {
		
		SAMRecordWrapper pair = null;
		
		if (read.getFirstOfPairFlag()) {
			pair = secondReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		} else {
			pair = firstReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		}
		
		if (pair != null && read.getCigar().getCigarElements().size() > 0 && pair.getSamRecord().getCigar().getCigarElements().size() > 0) {
			
			// Looking for something like:
			//     -------->
			//  <--------
			SAMRecord first = null;
			SAMRecord second = null;
			if (read.getReadNegativeStrandFlag()) {
				first = read;
				second = pair.getSamRecord();
			} else {
				first = pair.getSamRecord();
				second = read;
			}
			
			CigarElement firstElement = first.getCigar().getFirstCigarElement();
			CigarElement lastElement = second.getCigar().getLastCigarElement();
			
			if (firstElement.getOperator() == CigarOperator.S && lastElement.getOperator() == CigarOperator.S) {
				// We likely have read through into adapter here.
				hasReadThrough = true;
			}
		}
	}

	return hasReadThrough;
}
 

@Test
public void testMappedToMultipleStrands() throws Exception {
    final File outputMappedToMultipleStands = File.createTempFile("mappedToMultipleStrands", ".sam");
    outputMappedToMultipleStands.deleteOnExit();

    doMergeAlignment(mergingUnmappedBam,
            Collections.singletonList(multipleStrandsAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, outputMappedToMultipleStands,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(outputMappedToMultipleStands);

    for (final SAMRecord sam : result) {
        if (sam.getReadName().equals("test:1") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag() && sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReadString(), "TTTACTGATGTTATGACCATTACTCCGAAAGTGCCAAGATCATGAAGGGCAAGGAGAGAGTGGGATCCCCGGGTAC", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "GTACCCGGGGATCCCACTCTCTCCTTGCCCTTCATGATCTTGGCACTTTCGGAGTAATGGTCATAACATCAGTAAA", "Read aligned to positive strand has unexpected bases.");
            }
        }

        if (sam.getReadName().equals("test:2") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag() && sam.getSecondOfPairFlag()) {
                Assert.assertEquals(sam.getReadString(), "TTATTCACTTAGTGTGTTTTTCCTGAGAACTTGCTATGTGTTAGGTCCTAGGCTGGGTGGGATCCTCTAGAGTCGA", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "TCGACTCTAGAGGATCCCACCCAGCCTAGGACCTAACACATAGCAAGTTCTCAGGAAAAACACACTAAGTGAATAA", "Read aligned to positive strand has unexpected bases.");
            }
        }

        if (sam.getReadName().equals("test:5") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getReadString(), "AGTTTTGGTTTGTCAGACCCAGCCCTGGGCACAGATGAGGAATTCTGGCTTCTCCTCCTGTGGGATCCCCGGGTAC", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "GTACCCGGGGATCCCACAGGAGGAGAAGCCAGAATTCCTCATCTGTGCCCAGGGCTGGGTCTGACAAACCAAAACT", "Read aligned to positive strand has unexpected bases.");
            }
        }
    }
}
 

/**
 * Takes an individual SAMRecord and applies the set of changes/reversions to it that
 * have been requested by program level options.
 */
public void revertSamRecord(final SAMRecord rec) {
    if (RESTORE_ORIGINAL_QUALITIES) {
        final byte[] oq = rec.getOriginalBaseQualities();
        if (oq != null) {
            rec.setBaseQualities(oq);
            rec.setOriginalBaseQualities(null);
        }
    }

    if (REMOVE_DUPLICATE_INFORMATION) {
        rec.setDuplicateReadFlag(false);
    }

    if (REMOVE_ALIGNMENT_INFORMATION) {
        if (rec.getReadNegativeStrandFlag()) {
            rec.reverseComplement(true);
            rec.setReadNegativeStrandFlag(false);
        }

        // Remove all alignment based information about the read itself
        rec.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
        rec.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);

        rec.setInferredInsertSize(0);
        rec.setNotPrimaryAlignmentFlag(false);
        rec.setProperPairFlag(false);
        rec.setReadUnmappedFlag(true);

        // Then remove any mate flags and info related to alignment
        rec.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setMateNegativeStrandFlag(false);
        rec.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setMateUnmappedFlag(rec.getReadPairedFlag());

        if (RESTORE_HARDCLIPS) {
            String hardClippedBases = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG);
            String hardClippedQualities = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG);
            if (hardClippedBases != null && hardClippedQualities != null) {
                // Record has already been reverse complemented if this was on the negative strand
                rec.setReadString(rec.getReadString() + hardClippedBases);
                rec.setBaseQualities(SAMUtils.fastqToPhred(SAMUtils.phredToFastq(rec.getBaseQualities()) + hardClippedQualities));

                // Remove hard clipping storage tags
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG, null);
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG, null);
            }
        }

        // And then remove any tags that are calculated from the alignment
        ATTRIBUTE_TO_CLEAR.forEach(tag -> rec.setAttribute(tag, null));
    }
}
 

private void writeRecord(final SAMRecord read, final Integer mateNumber, final FastqWriter writer,
                         final int basesToTrim, final Integer maxBasesToWrite) {
    final String seqHeader = mateNumber == null ? read.getReadName() : read.getReadName() + "/" + mateNumber;
    String readString = read.getReadString();
    String baseQualities = read.getBaseQualityString();

    // If we're clipping, do the right thing to the bases or qualities
    if (CLIPPING_ATTRIBUTE != null) {
        Integer clipPoint = (Integer) read.getAttribute(CLIPPING_ATTRIBUTE);
        if (clipPoint != null && clipPoint < CLIPPING_MIN_LENGTH) {
            clipPoint = Math.min(readString.length(), CLIPPING_MIN_LENGTH);
        }

        if (clipPoint != null) {
            if (CLIPPING_ACTION.equalsIgnoreCase(CLIP_TRIM)) {
                readString = clip(readString, clipPoint, null, !read.getReadNegativeStrandFlag());
                baseQualities = clip(baseQualities, clipPoint, null, !read.getReadNegativeStrandFlag());
            } else if (CLIPPING_ACTION.equalsIgnoreCase(CLIP_TO_N)) {
                readString = clip(readString, clipPoint, CLIP_TO_N.charAt(0), !read.getReadNegativeStrandFlag());
            } else {
                final char newQual = SAMUtils.phredToFastq(new byte[]{(byte) Integer.parseInt(CLIPPING_ACTION)}).charAt(0);
                baseQualities = clip(baseQualities, clipPoint, newQual, !read.getReadNegativeStrandFlag());
            }
        }
    }

    if (RE_REVERSE && read.getReadNegativeStrandFlag()) {
        readString = SequenceUtil.reverseComplement(readString);
        baseQualities = StringUtil.reverseString(baseQualities);
    }

    if (basesToTrim > 0) {
        readString = readString.substring(basesToTrim);
        baseQualities = baseQualities.substring(basesToTrim);
    }

    // Perform quality trimming if desired, making sure to leave at least one base!
    if (QUALITY != null) {
        final byte[] quals = SAMUtils.fastqToPhred(baseQualities);
        final int qualityTrimIndex = Math.max(1, TrimmingUtil.findQualityTrimPoint(quals, QUALITY));
        if (qualityTrimIndex < quals.length) {
            readString = readString.substring(0, qualityTrimIndex);
            baseQualities = baseQualities.substring(0, qualityTrimIndex);
        }
    }

    if (maxBasesToWrite != null && maxBasesToWrite < readString.length()) {
        readString = readString.substring(0, maxBasesToWrite);
        baseQualities = baseQualities.substring(0, maxBasesToWrite);
    }

    writer.write(new FastqRecord(seqHeader, readString, "", baseQualities));
}
 

@Test
public void testMerger() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReader result = SamReaderFactory.makeDefault().open(output);
    Assert.assertEquals(result.getFileHeader().getSequenceDictionary().getSequences().size(), 8,
            "Number of sequences did not match");
    SAMProgramRecord pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "0");
    Assert.assertEquals(pg.getProgramVersion(), "1.0");
    Assert.assertEquals(pg.getCommandLine(), "align!");
    Assert.assertEquals(pg.getProgramName(), "myAligner");

    final SAMReadGroupRecord rg = result.getFileHeader().getReadGroups().get(0);
    Assert.assertEquals(rg.getReadGroupId(), "0");
    Assert.assertEquals(rg.getSample(), "Hi,Mom!");

    Assert.assertEquals(result.getFileHeader().getSortOrder(), SAMFileHeader.SortOrder.coordinate);

    for (final SAMRecord sam : result) {
        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_adapter") ||
                sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " +
                        sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " +
                        sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(sam.getCigarString().indexOf("S") == -1,
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }

        final boolean pos = !sam.getReadNegativeStrandFlag();
        Assert.assertEquals(sam.getAttribute("aa"), pos ? "Hello" : "olleH");
        Assert.assertEquals(sam.getAttribute("ab"), pos ? "ATTCGG" : "CCGAAT");
        Assert.assertEquals(sam.getAttribute("ac"), pos ? new byte[] {1,2,3} : new byte[] {3,2,1});
        Assert.assertEquals(sam.getAttribute("as"), pos ? new short[]{1,2,3} : new short[]{3,2,1});
        Assert.assertEquals(sam.getAttribute("ai"), pos ? new int[]  {1,2,3} : new int[]  {3,2,1});
        Assert.assertEquals(sam.getAttribute("af"), pos ? new float[]{1,2,3} : new float[]{3,2,1});
    }

    // Quick test to make sure the program record gets picked up from the file if not specified
    // on the command line.
    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            null, null, null, null,
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    CloserUtil.close(result);

    result = SamReaderFactory.makeDefault().open(output);
    pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "1",
            "Program group ID not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramVersion(), "2.0",
            "Program version not picked up correctly from aligned BAM");
    Assert.assertNull(pg.getCommandLine(),
            "Program command line not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramName(), "Hey!",
            "Program name not picked up correctly from aligned BAM");

    CloserUtil.close(result);
}
 

@Test
public void testShortFragmentHardClipping() throws IOException {
    final File output = File.createTempFile("testShortFragmentClipping", ".sam");
    output.deleteOnExit();
    final File alignedSam = new File(TEST_DATA_DIR, "hardclip.aligned.sam");
    final File unmappedSam = new File(TEST_DATA_DIR, "hardclip.unmapped.sam");
    final File ref = new File(TEST_DATA_DIR, "cliptest.fasta");

    final List<String> args = Arrays.asList(
            "UNMAPPED_BAM=" + unmappedSam.getAbsolutePath(),
            "ALIGNED_BAM=" + alignedSam.getAbsolutePath(),
            "OUTPUT=" + output.getAbsolutePath(),
            "REFERENCE_SEQUENCE=" + ref.getAbsolutePath(),
            "HARD_CLIP_OVERLAPPING_READS=true"
            );

    Assert.assertEquals(runPicardCommandLine(args), 0);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    final Map<String, SAMRecord> firstReadEncountered = new HashMap<String, SAMRecord>();
    for (final SAMRecord rec : result) {
        final SAMRecord otherEnd = firstReadEncountered.get(rec.getReadName());
        if (otherEnd == null) {
            firstReadEncountered.put(rec.getReadName(), rec);
        } else {
            final int fragmentStart = Math.min(rec.getAlignmentStart(), otherEnd.getAlignmentStart());
            final int fragmentEnd = Math.max(rec.getAlignmentEnd(), otherEnd.getAlignmentEnd());
            final String[] readNameFields = rec.getReadName().split(":");
            // Read name of each pair includes the expected fragment start and fragment end positions.
            final int expectedFragmentStart = Integer.parseInt(readNameFields[1]);
            final int expectedFragmentEnd = Integer.parseInt(readNameFields[2]);
            Assert.assertEquals(fragmentStart, expectedFragmentStart, rec.getReadName());
            Assert.assertEquals(fragmentEnd, expectedFragmentEnd, rec.getReadName());
            if (readNameFields[0].equals("FR_clip")) {
                Assert.assertEquals(rec.getCigarString(), rec.getReadNegativeStrandFlag()? "20H56M" : "56M20H");
                Assert.assertEquals(otherEnd.getCigarString(), otherEnd.getReadNegativeStrandFlag()? "20H56M" : "56M20H");

                if (!rec.getReadNegativeStrandFlag()) {
                    Assert.assertEquals(rec.getAttribute("XB"), "AGATCGGAAGAGCACACGTC");
                    Assert.assertEquals(rec.getAttribute("XQ"), "BBBBB?BBB?<?A?<7<<=9");
                } else {
                    Assert.assertEquals(rec.getAttribute("XB"), "AGATCGGAAGAGCGTCGTGT");
                    Assert.assertEquals(rec.getAttribute("XQ"), "BCFD=@CBBADCF=CC:CCD");
                }
            } else {
                Assert.assertEquals(rec.getCigarString(), "76M");
                Assert.assertEquals(otherEnd.getCigarString(), "76M");
            }
        }
    }
}
 

private boolean readAnnotationMatchStrand (final Gene gene, final SAMRecord r) {
	boolean negativeStrandRead = r.getReadNegativeStrandFlag();
	boolean geneNegativeStrand = gene.isNegativeStrand();
	return negativeStrandRead==geneNegativeStrand;
}