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

The following examples show how to use htsjdk.samtools.SAMRecord#getSAMString() . 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: SomaticLocusCaller.java    From abra2 with MIT License 5 votes vote down vote up 
private boolean hasIndel(SAMRecord read, LocusInfo locus) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= locus.posStop && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case I:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case S:
				readPos += elem.getLength();
				break;
			case D:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				readPos += elem.getLength();
				refPosInRead += elem.getLength();
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return false;
}
Example 2
Source File: SomaticLocusCaller.java    From abra2 with MIT License 5 votes vote down vote up 
private Object[] getBaseAndQualAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return new Object[] { read.getReadString().charAt(readPos) , (int) read.getBaseQualities()[readPos] };
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return new Object[] { 'N', (int) 0 };
}
Example 3
Source File: CadabraProcessor.java    From abra2 with MIT License 5 votes vote down vote up 
private Character getBaseAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return read.getReadString().charAt(readPos);
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return null;
}
Example 4
Source File: SimpleAlleleCounter.java    From abra2 with MIT License 5 votes vote down vote up 
private Pair<Character,Character> getBaseAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return new Pair<Character, Character>(read.getReadString().charAt(readPos), read.getBaseQualityString().charAt(readPos));
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return null;
}
Example 5
Source File: ReorderSam.java    From picard with MIT License 5 votes vote down vote up 
/**
 * Low-level helper function that returns the new reference index for oldIndex according to the
 * ordering map newOrder.  Read is provided in case an error occurs, so that an informative message
 * can be made.
 */
private int newOrderIndex(SAMRecord read, int oldIndex, Map<Integer, Integer> newOrder) {
    if (oldIndex == NO_ALIGNMENT_REFERENCE_INDEX) {
        return NO_ALIGNMENT_REFERENCE_INDEX; // unmapped read
    } else {
        final Integer n = newOrder.get(oldIndex);

        if (n == null) {
            throw new PicardException("BUG: no mapping found for read " + read.getSAMString());
        } else {
            return n;
        }
    }
}
Example 6
Source File: PreprocessingTools.java    From halvade with GNU General Public License v3.0 5 votes vote down vote up 
public int streamElPrep(Reducer.Context context, String output, String rg, 
            int threads, SAMRecordIterator SAMit, 
            SAMFileHeader header, String dictFile, boolean updateRG, boolean keepDups, String RGID) throws InterruptedException, IOException, QualityException {
        long startTime = System.currentTimeMillis();
        String customArgs = HalvadeConf.getCustomArgs(context.getConfiguration(), "elprep", "");  
        String[] command = CommandGenerator.elPrep(bin, "/dev/stdin", output, threads, true, rg, null, !keepDups, customArgs);
//        runProcessAndWait(command);
        ProcessBuilderWrapper builder = new ProcessBuilderWrapper(command, null);
        builder.startProcess(true);        
        BufferedWriter localWriter = builder.getSTDINWriter();
        
        // write header
        final StringWriter headerTextBuffer = new StringWriter();
        new SAMTextHeaderCodec().encode(headerTextBuffer, header);
        final String headerText = headerTextBuffer.toString();
        localWriter.write(headerText, 0, headerText.length());
        
        
        SAMRecord sam;
        int reads = 0;
        while(SAMit.hasNext()) {
            sam = SAMit.next();
            if(updateRG)
                sam.setAttribute(SAMTag.RG.name(), RGID);
            String samString = sam.getSAMString();
            localWriter.write(samString, 0, samString.length());
            reads++;
        }
        localWriter.flush();
        localWriter.close();
                
        int error = builder.waitForCompletion();
        if(error != 0)
            throw new ProcessException("elPrep", error);
        long estimatedTime = System.currentTimeMillis() - startTime;
        Logger.DEBUG("estimated time: " + estimatedTime / 1000);
        if(context != null)
            context.getCounter(HalvadeCounters.TIME_ELPREP).increment(estimatedTime);
        return reads;
    }
Example 7
Source File: CgSamBamSequenceDataSource.java    From rtg-tools with BSD 2-Clause "Simplified" License 4 votes vote down vote up 
/**
 * Unroll both the read bases and quality data for a CG alignment
 * @param rec the alignment
 * @return the unrolled read
 */
public static SamSequence unrollCgRead(SAMRecord rec) {
  final int projectedSplit = rec.getAlignmentStart() * ((rec.getFlags() & SamBamConstants.SAM_MATE_IS_REVERSE) != 0 ? 1 : -1);
  final byte flags = SamSequence.getFlags(rec);
  final byte[] expandedRead;
  final byte[] expandedQual;
  if (rec.getReadUnmappedFlag()) {
    expandedRead = rec.getReadBases();
    expandedQual = rec.getBaseQualities();
  } else {
    final String gc = rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_RAW_READ_INSTRUCTIONS);
    if (gc == null) {
      throw new NoTalkbackSlimException("SAM Record does not contain CGI read reconstruction attribute: " + rec.getSAMString());
    }
    final byte[] gq = FastaUtils.asciiToRawQuality(SamUtils.allowEmpty(rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_OVERLAP_QUALITY)));
    final byte[] gs = SamUtils.allowEmpty(rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_OVERLAP_BASES)).getBytes();
    final boolean legacyLegacy = gq.length == gs.length / 2;
    expandedRead = unrollLegacyRead(rec.getReadBases(), gs, gc);
    if (expandedRead == null) {
      throw new NoTalkbackSlimException("Could not reconstruct read bases for record: " + rec.getSAMString());
    }
    if (rec.getBaseQualities().length == 0) {
      expandedQual = rec.getBaseQualities();
    } else {
      if (!legacyLegacy && gq.length != gs.length) {
        throw new NoTalkbackSlimException("Unexpected length of CG quality information: " + rec.getSAMString());
      }
      final byte[] samQualities = rec.getBaseQualities();
      if (legacyLegacy) {
        expandedQual = unrollLegacyLegacyQualities(samQualities, gq, gc);
      } else {
        final byte[] bytes = unrollLegacyRead(samQualities, gq, gc);
        expandedQual = bytes;
      }
    }
  }

  final byte[] readBytes = CgUtils.unPad(expandedRead, !rec.getReadNegativeStrandFlag());
  final byte[] readQual = CgUtils.unPad(expandedQual, !rec.getReadNegativeStrandFlag());

  return new SamSequence(rec.getReadName(), readBytes, readQual, flags, projectedSplit);
}