Java Code Examples for htsjdk.samtools.CigarOperator#MATCH_OR_MISMATCH
The following examples show how to use
htsjdk.samtools.CigarOperator#MATCH_OR_MISMATCH .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: IndelShifter.java From abra2 with MIT License | 4 votes |
|
public Cigar shiftAllIndelsLeft(int refStart, int refEnd, String chromosome, Cigar cigar, String seq, CompareToReference2 c2r) {
try {
List<CigarElement> elems = new ArrayList<CigarElement>(cigar.getCigarElements());
int elemSize = elems.size();
for (int i=1; i<elemSize-1; i++) {
int refOffset = 0;
int readOffset = 0;
for (int j=0; j<i-1; j++) {
refOffset += getRefOffset(elems.get(j));
readOffset += getReadOffset(elems.get(j));
}
CigarElement prev = elems.get(i-1);
CigarElement elem = elems.get(i);
CigarElement next = elems.get(i+1);
if ((elem.getOperator() == CigarOperator.DELETION || elem.getOperator() == CigarOperator.INSERTION) &&
prev.getOperator() == CigarOperator.MATCH_OR_MISMATCH && next.getOperator() == CigarOperator.MATCH_OR_MISMATCH) {
// subset cigar here and attempt to shift
Cigar subCigar = new Cigar(Arrays.asList(prev, elem, next));
String subSeq = seq.substring(readOffset, readOffset+subCigar.getReadLength());
Cigar newCigar = shiftIndelsLeft(refStart + refOffset, refStart + refOffset + subCigar.getReferenceLength(),
chromosome, subCigar, subSeq, c2r);
//TODO: Merge abutting indels if applicable
elems.set(i-1, newCigar.getCigarElement(0));
elems.set(i, newCigar.getCigarElement(1));
if (newCigar.getCigarElements().size() == 3) {
elems.set(i+1, newCigar.getCigarElement(2));
} else {
elems.remove(i+1);
elemSize -= 1;
}
}
}
mergeAbuttingElements(elems);
return new Cigar(elems);
} catch (RuntimeException e) {
e.printStackTrace();
Logger.error("Error on: " + refStart + ", " + refEnd + "," + chromosome + ", " + cigar + ", " + seq);
throw e;
}
}
Example 2
| Source File: Utils.java From cramtools with Apache License 2.0 | 4 votes |
|
/**
* A rip off samtools bam_md.c
*
* @param record
* @param ref
* @param calcMD
* @param calcNM
*/
public static void calculateMdAndNmTags(SAMRecord record, byte[] ref, boolean calcMD, boolean calcNM) {
if (!calcMD && !calcNM)
return;
Cigar cigar = record.getCigar();
List<CigarElement> cigarElements = cigar.getCigarElements();
byte[] seq = record.getReadBases();
int start = record.getAlignmentStart() - 1;
int i, x, y, u = 0;
int nm = 0;
StringBuffer str = new StringBuffer();
int size = cigarElements.size();
for (i = y = 0, x = start; i < size; ++i) {
CigarElement ce = cigarElements.get(i);
int j, l = ce.getLength();
CigarOperator op = ce.getOperator();
if (op == CigarOperator.MATCH_OR_MISMATCH || op == CigarOperator.EQ || op == CigarOperator.X) {
for (j = 0; j < l; ++j) {
int z = y + j;
if (ref.length <= x + j)
break; // out of boundary
int c1 = 0;
int c2 = 0;
// try {
c1 = seq[z];
c2 = ref[x + j];
if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) {
// a match
++u;
} else {
str.append(u);
str.appendCodePoint(ref[x + j]);
u = 0;
++nm;
}
}
if (j < l)
break;
x += l;
y += l;
} else if (op == CigarOperator.DELETION) {
str.append(u);
str.append('^');
for (j = 0; j < l; ++j) {
if (ref[x + j] == 0)
break;
str.appendCodePoint(ref[x + j]);
}
u = 0;
if (j < l)
break;
x += l;
nm += l;
} else if (op == CigarOperator.INSERTION || op == CigarOperator.SOFT_CLIP) {
y += l;
if (op == CigarOperator.INSERTION)
nm += l;
} else if (op == CigarOperator.SKIPPED_REGION) {
x += l;
}
}
str.append(u);
if (calcMD)
record.setAttribute(SAMTag.MD.name(), str.toString());
if (calcNM)
record.setAttribute(SAMTag.NM.name(), nm);
}
Example 3
| Source File: AlignmentsTags.java From cramtools with Apache License 2.0 | 4 votes |
|
/**
* Same as {@link htsjdk.samtools.util.SequenceUtil#calculateMdAndNmTags}
* but allows for reference excerpt instead of a full reference sequence.
*
* @param record
* SAMRecord to inject NM and MD tags into
* @param ref
* reference sequence bases, may be a subsequence starting at
* refOffest
* @param refOffset
* array offset of the reference bases, 0 is the same as the
* whole sequence. This value will be subtracted from the
* reference array index in calculations.
* @param calcMD
* calculate MD tag if true
* @param calcNM
* calculate NM tag if true
*/
public static void calculateMdAndNmTags(final SAMRecord record, final byte[] ref, final int refOffset,
final boolean calcMD, final boolean calcNM) {
if (!calcMD && !calcNM)
return;
final Cigar cigar = record.getCigar();
final List<CigarElement> cigarElements = cigar.getCigarElements();
final byte[] seq = record.getReadBases();
final int start = record.getAlignmentStart() - 1;
int i, x, y, u = 0;
int nm = 0;
final StringBuilder str = new StringBuilder();
final int size = cigarElements.size();
for (i = y = 0, x = start; i < size; ++i) {
final CigarElement ce = cigarElements.get(i);
int j;
final int length = ce.getLength();
final CigarOperator op = ce.getOperator();
if (op == CigarOperator.MATCH_OR_MISMATCH || op == CigarOperator.EQ || op == CigarOperator.X) {
for (j = 0; j < length; ++j) {
final int z = y + j;
if (refOffset + ref.length <= x + j)
break; // out of boundary
int c1 = 0;
int c2 = 0;
// try {
c1 = seq[z];
c2 = ref[x + j - refOffset];
if ((c1 == c2) || c1 == 0) {
// a match
++u;
} else {
str.append(u);
str.appendCodePoint(ref[x + j - refOffset]);
u = 0;
++nm;
}
}
if (j < length)
break;
x += length;
y += length;
} else if (op == CigarOperator.DELETION) {
str.append(u);
str.append('^');
for (j = 0; j < length; ++j) {
if (ref[x + j - refOffset] == 0)
break;
str.appendCodePoint(ref[x + j - refOffset]);
}
u = 0;
if (j < length)
break;
x += length;
nm += length;
} else if (op == CigarOperator.INSERTION || op == CigarOperator.SOFT_CLIP) {
y += length;
if (op == CigarOperator.INSERTION)
nm += length;
} else if (op == CigarOperator.SKIPPED_REGION) {
x += length;
}
}
str.append(u);
if (calcMD)
record.setAttribute(SAMTag.MD.name(), str.toString());
if (calcNM)
record.setAttribute(SAMTag.NM.name(), nm);
}
