Examples of net.sf.samtools.SAMRecord

net.sf.samtools.SAMRecord

    // 0x800 = supplementary alignment from SAM spec v1.5
    return ((read.getFlags() & 0x800)  == 0) && (!read.getNotPrimaryAlignmentFlag()); 
  }
  
  private SAMRecord combineChimericReads(SAMRecord read1, SAMRecord read2) {
    SAMRecord left = null;
    SAMRecord right = null;
    
    if (read1.getAlignmentStart() < read2.getAlignmentStart()) {
      left = read1;
      right = read2;
    } 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?
      
      return combinedRead;

View Full Code Here

    }
  }
  
  private void pruneLikelyInserts(List<SAMRecord> sortedReads) {
    if (sortedReads.size() == 3) {
      SAMRecord first  = sortedReads.get(0);
      SAMRecord middle = sortedReads.get(1);
      SAMRecord last   = sortedReads.get(2);
      
//      if ((middle.getAlignmentStart() <= first.getAlignmentEnd()) && (middle.getAlignmentEnd() >= last.getAlignmentStart())) {
//        sortedReads.remove(middle);
//      }

View Full Code Here

      
      if (readCache.isEmpty()) {
        shouldReadFromFile = true;
      }
      else {
        SAMRecord last = readCache.get(readCache.size()-1);
        if (getAlignmentStart(last) <= currentPos && last.getReferenceName().equals(currentChr)) {
          shouldReadFromFile = true;
        }
      }
      
      while (shouldReadFromFile && samIter.hasNext()) {
        SAMRecord read = samIter.next();
        
        if (readCache.isEmpty() && !read.getReferenceName().equals(currentChr)) {
          currentChr = read.getReferenceName();
          currentPos = getAlignmentStart(read);
        }
        
        readCache.add(read);
        
        if (getAlignmentStart(read) > currentPos || !read.getReferenceName().equals(currentChr)) {
          shouldReadFromFile = false;
        }
      }      
    }

View Full Code Here

      
      String nextChr = null;
      int nextPos = -1;
      
      while (cacheIter.hasNext()) {
        SAMRecord read = cacheIter.next();
        
        if (read.getAlignmentEnd() < currentPos && read.getReferenceName().equals(currentChr)) {
          // We've gone past the end of this read, so remove from cache.
          cacheIter.remove();
        } else if (getAlignmentStart(read) <= currentPos && read.getAlignmentEnd() >= currentPos) {
          // This read spans the current locus of interest.
          reads.add(read);
        } else {
          // This read is beyond the current locus.
          if (nextChr == null) {
            nextChr = read.getReferenceName();
            nextPos = getAlignmentStart(read);
          }
        }
      }

View Full Code Here

    Map<String, Set<String>> posCounts = new HashMap<String, Set<String>>();


    CloseableIterator<SAMRecord> iter = reader.queryOverlapping(region.getSeqname(), (int) region.getStart(), (int) region.getEnd());
    
    while (iter.hasNext()) {
      SAMRecord contig = iter.next();
      
      List<SAMRecord> contigChunks = chunkRead(contig);
      for (SAMRecord chunk : contigChunks) {
        chunks.put(chunk.getReadString(), chunk);
      }

View Full Code Here

      endIdx = Math.min(startIdx + readLength*2, contig.getReadLength()-1);
//      int endIdx = Math.min(end - contig.getAlignmentStart(), startIdx+contig.getReadLength()-1); // inclusive
      String bases = contig.getReadString().substring(startIdx, endIdx+1);
      Cigar cigar = SAMRecordUtils.subset(contig.getCigar(), startIdx, endIdx);
      
      SAMRecord chunk = cloneRead(contig);
      // Using a static variable because contigs may overlap multiple regions and cause dupes
      chunk.setReadName(chunk.getReadName() + "_" + chunkIdx++);
      chunk.setAlignmentStart(actualPos + calcIndelOffset(startIdx, contig.getCigar()));
      chunk.setCigar(cigar);
      chunk.setReadString(bases);
      chunk.clearAttributes();
//      chunk.setAttribute("ZZ", contig.getCigarString());
      
      chunk.setAttribute("ZZ", contig.getReferenceName() + ":" + contig.getAlignmentStart() +
          ":" + contig.getCigarString()); 
      
      chunks.add(chunk);
      
      pos += readLength;

View Full Code Here

  }
  
  private SAMRecord getOrigRecord(SAMRecord read, SAMFileHeader samHeader) {
    String origSamStr = read.getReadName();
    origSamStr = origSamStr.replace(Sam2Fastq.FIELD_DELIMITER, "\t");
    SAMRecord orig;
    try {
      orig = samStringReader.getRead(origSamStr);
    } catch (RuntimeException e) {
      System.out.println("Error processing: [" + origSamStr + "]");
      System.out.println("Contig read: [" + read.getSAMString() + "]");
      e.printStackTrace();
      throw e;
    }
    orig.setHeader(samHeader);
    
//    orig.setReadString(read.getReadString());
//    orig.setBaseQualityString(read.getBaseQualityString());


    return orig;

View Full Code Here

      }
    }
    




    SAMRecord lastRead = reads.get(reads.size()-1);
    int regionStart = reads.get(0).getAlignmentStart();
    int regionEnd = lastRead.getAlignmentEnd() > 0 ? lastRead.getAlignmentEnd() : lastRead.getAlignmentStart();
    
    String output = "region_" + reads.get(0).getReferenceName() + "_" + regionStart + "_" + regionEnd;
    String contigs = "";
    
    // Make this set of reads eligible for GC

View Full Code Here

          
          int candidateReadCount = 0;
          
          while (iter.hasNext() && !isMaxReadsForRegionExceeded) {
            
            SAMRecord read = iter.next();
            readCount++;
                        
            // Don't allow same read to be counted twice.
            if ( (!realigner.isFiltered(read)) && 
               (!read.getDuplicateReadFlag()) && 
               (!read.getReadFailsVendorQualityCheckFlag()) &&
               (read.getMappingQuality() >= realigner.getMinMappingQuality() || read.getReadUnmappedFlag()) &&
               //(Sam2Fastq.isPrimary(read)) &&
//               (!isHardClipped(read)) &&
               ((!checkForDupes) || (!readIds.contains(getIdentifier(read))))) {
              
              if (read.getReadString().length() > readLength) {
                reader.close();
                throw new IllegalArgumentException("Maximum read length of: " + readLength +
                    " exceeded for: " + read.getSAMString());
              }
  
              Integer numBestHits = (Integer) read.getIntegerAttribute("X0");
              boolean hasAmbiguousInitialAlignment = numBestHits != null && numBestHits > 1;
              
              if (!hasAmbiguousInitialAlignment) {
                if (!checkForDupes) {
                  readIds.add(getIdentifier(read));
                }
                
                if (!isAssemblyCandidate && isAssemblyTriggerCandidate(read, c2r)) {
                  candidateReadCount++;
                }
                
                reads.add(read);
                
                if (shouldSearchForSv && isSvCandidate(read, region)) {
                  svCandidates.add(new Position(read.getMateReferenceName(), read.getMateAlignmentStart(), input));
                }
                
                if (reads.size() > maxReadsPerSample) {
                  isMaxReadsForRegionExceeded = true;
                  break;
                }
              }
            }
          }
          
          reader.close();
          
          if (candidateReadCount > minCandidateCount(reads.size(), region)) {
            isAssemblyCandidate = true;
          }
        }
        
        if (reads.size() < minReadCount) {
          minReadCount = reads.size();
        }
      }
      
      if (isMaxReadsForRegionExceeded) {
        isAssemblyCandidate = false;
        shouldSearchForSv = false;
        System.out.println("Max reads exceeded for region: " + prefix);
      }
      
      readIds = null;
      
      StringBuffer readBuffer = new StringBuffer();
      
      if (isAssemblyCandidate) {
        
        int downsampleTarget = desiredNumberOfReads(regions);
        
        for (List<SAMRecord> reads : readsList) {
          // Default to always keep
          double keepProbability = 1.1;
          
          if (reads.size() > downsampleTarget) {
            keepProbability = (double) downsampleTarget / (double) reads.size();
          }
          
          Random random = new Random(1);
          
          for (SAMRecord read : reads) {
            if (random.nextDouble() < keepProbability) {
              readBuffer.append(read.getReadNegativeStrandFlag() ? "1" : "0");
              
              if (read.getReadString().length() == readLength) {
                readBuffer.append(read.getReadString());
                readBuffer.append(read.getBaseQualityString());
              } else {
                StringBuffer basePadding = new StringBuffer();
                StringBuffer qualPadding = new StringBuffer();
                
                for (int i=0; i<readLength-read.getReadString().length(); i++) {
                  basePadding.append('N');
                  qualPadding.append('!');
                }
                
                readBuffer.append(read.getReadString() + basePadding.toString());
                readBuffer.append(read.getBaseQualityString() + qualPadding.toString());              
              }
            }
          }
          
          // Make this set of reads eligible for GC

View Full Code Here

            Cigar newCigar = shiftCigarLeft(read.getCigar(), i);
            
            String shiftedReadAltRef = c2r.getAlternateReference(read, newCigar);
            
            if ((shiftedReadAltRef != null) && (origReadAltRef.equals(shiftedReadAltRef))) {          
              SAMRecord newRead = cloneRead(read);
              newRead.setCigar(newCigar);
              newRead.setAttribute("IS", read.getCigarString());
              return newRead;
            }
          }
        }
      }

View Full Code Here

0 1 2 3 4 5 6 7

TOP

Related Classes of net.sf.samtools.SAMRecord

abra.bamsplitter.BamSplitter

abra.bamsplitter.BamSplitterThread

abra.BetaPairValidatingRealignmentWriter

abra.cadabra.ForwardShiftInsertIterator

abra.cadabra.ReadLocusReader$ReadLocusIterator

abra.cadabra.SomaticLocusCaller

abra.CombineChimera3

abra.CombineChimeraTest

abra.ContigChopper

abra.IndelShifter

All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.