} else {
left = read2;
right = read1;
}
Cigar leftCigar = left.getCigar();
Cigar rightCigar = right.getCigar();
SAMRecord topHit = getTopHit(read1, read2);
if ((rightCigar.getCigarElement(0).getOperator() == CigarOperator.S) &&
(leftCigar.getCigarElement(left.getCigarLength()-1).getOperator() == CigarOperator.S) &&
(topHit != null)) {
List<CigarElement> leftElements = new ArrayList<CigarElement>();
List<CigarElement> rightElements = new ArrayList<CigarElement>();
// Drop trailing S on left side
for (int i=0; i<left.getCigar().numCigarElements()-1; i++) {
leftElements.add(leftCigar.getCigarElement(i));
}
// Drop leading S on right side
for (int i=1; i<rightCigar.numCigarElements(); i++) {
rightElements.add(rightCigar.getCigarElement(i));
}
// Encountered in dream test data at: chr20 22137016. Not clear how to interpret.
if (rightElements.get(0).getOperator() == CigarOperator.INSERTION) {
return null;
}
// If total element length is longer than the read, then trim first element
// on the left side of the indel (this is likely a deletion??)
// if (totalLength > read1.getReadLength()) {
// // Trim from the right side of the rightmos block of the left elements
// int trimLength = totalLength - read1.getReadLength();
// trimmedElemLength = trimRightmostElement(rightElements, trimLength);
// }
//
// if (trimmedElemLength < 0) {
// // We don't currently handle the case where we need to trim more
// // than the length of the element
// return null;
// }
List<ReadBlock> leftBlocks = ReadBlock.getReadBlocks(left);
ReadBlock lastLeftBlock = leftBlocks.get(leftBlocks.size()-2);
List<ReadBlock> rightBlocks = ReadBlock.getReadBlocks(right);
ReadBlock firstRightBlock = rightBlocks.get(1);
// Confirm no shared bases in read
// int leftStop = lastLeftBlock.getReadStart() + lastLeftBlock.getLength() - 1;
// int rightStart = firstRightBlock.getReadStart();
int leftStop = lastLeftBlock.getReferenceStop();
int rightStart = firstRightBlock.getReferenceStart();
int leftReadStop = lastLeftBlock.getReadStart() + lastLeftBlock.getLength() - 1;
int rightReadStart = firstRightBlock.getReadStart();
int trimLength = 0;
if ((leftStop >= rightStart) || (leftReadStop >= rightReadStart)) {
trimLength = Math.max(leftStop - rightStart + 1, leftReadStop - rightReadStart + 1);
int trimmedElemLength = trimRightmostElement(leftElements, trimLength);
if (trimmedElemLength < 1) {
// This isn't really an indel. Just alternate mappings.
// System.out.println("Element trimmed too far. read1: [" + read1.getSAMString() +
// "] read2: [" + read2.getSAMString() + "]");
return null;
}
}
// Build indel
int leftAlignmentStop = lastLeftBlock.getReferenceStop();
leftAlignmentStop -= trimLength;
int rightAlignmentStart = firstRightBlock.getReferenceStart();
leftReadStop -= trimLength;
int alignmentGap = rightAlignmentStart - leftAlignmentStop - 1;
int readGap = rightReadStart - leftReadStop - 1;
CigarElement indelElement;
if ((alignmentGap > 0) && (readGap == 0)) {
// Deletion
indelElement = new CigarElement(alignmentGap, CigarOperator.DELETION);
} else {
int totalLength = this.getTotalLength(leftElements, rightElements);
// int insertLen = read1.getReadLength() - totalLength;
int insertLen = read1.getReadLength() - totalLength - alignmentGap;
int insertLen2 = readGap - alignmentGap;
if (insertLen != insertLen2) {
// Not really an insert
return null;
}
//TODO: Take a closer look at this
if (insertLen < 1) {
return null;
}
// Right side may have skipped SNPs, so pad it.
if (alignmentGap > 0) {
this.padLeftmostElement(rightElements, alignmentGap);
}
//TODO: Necessary for mismatches?
/*
if ((readGap > 0) && ((totalLength + insertLen) < read1.getReadLength())) {
if (readGap != read1.getReadLength() - totalLength - insertLen) {
System.out.println("WARNING: Invalid read gap padding: " + read1.getSAMString());
}
this.padLeftmostElement(rightElements, readGap);
}
*/
indelElement = new CigarElement(insertLen, CigarOperator.INSERTION);
}
// Check to see if the indel is surrounded by non-clipped bases
if (minIndelBuffer > 0) {
if ((getNonSoftClippedLength(leftElements) < minIndelBuffer) ||
(getNonSoftClippedLength(rightElements) < minIndelBuffer)) {
return null;
}
}
List<CigarElement> elements = new ArrayList<CigarElement>();
elements.addAll(leftElements);
elements.add(indelElement);
elements.addAll(rightElements);
// Create combined read
SAMRecord combinedRead = cloneRead(topHit);
combinedRead.setAlignmentStart(leftBlocks.get(0).getReferenceStart());
combinedRead.setCigar(new Cigar(elements));
combinedRead.setMappingQuality((read1.getMappingQuality() + read2.getMappingQuality()) / 2);
// combinedRead.setMappingQuality(Math.min(read1.getMappingQuality(), read2.getMappingQuality()));
//TODO: Any other ancilliary info to set?