Examples of net.sf.samtools.SAMFileReader

net.sf.samtools.SAMFileReader
Class for reading and querying SAM/BAM files. Delegates to appropriate concrete implementation.

  private BufferedReader refReader;


  public void compare(String sam, String refFileName, int maxDiff) throws IOException, FileNotFoundException {
    refReader = new BufferedReader(new FileReader(refFileName));
    
    SAMFileReader reader = new SAMFileReader(new File(sam));
    reader.setValidationStringency(ValidationStringency.SILENT);
    
    int count = 0;
    int totalMapped = 0;
    
    for (SAMRecord read : reader) {
      if (!read.getReadUnmappedFlag()) {
        String seq = read.getReferenceName();
        if (!seq.equals(currSeqName)) {
          loadSeqRef(seq);
        }
        
//        String refStr = reference.substring(read.getAlignmentStart()-1, read.getAlignmentStart()+99);
//        if (!refStr.equals(read.getReadString())) {
//        StringBuffer diffStr = new StringBuffer();
        if (numDifferences(read) > maxDiff) {
          System.out.println("------------");
          System.out.println("read: " + read.getSAMString());
//          System.out.println("ref: " + refStr);
//          System.out.println("dif: " + diffStr.toString());
          count += 1;
        }
        
        totalMapped += 1;
      }
    }
    
    System.out.println("count: " + count + " out of: " + totalMapped);
    
    reader.close();
  }

View Full Code Here

      int readLength, CompareToReference2 c2r) {
    
    Map<String, List<Integer>> locations = new HashMap<String, List<Integer>>();
    
    for (String file : files) {
          SAMFileReader reader = new SAMFileReader(new File(file));
          reader.setValidationStringency(ValidationStringency.SILENT);


          header = reader.getFileHeader();
          init();
          
          for (SAMRecord read : reader) {
            if (isInRegion(read)) {
              if (read.getCigarString().contains("I") || read.getCigarString().contains("D")) {
                addLocation(read, locations);
              } else if (read.getCigarString().contains("S") && c2r.numHighQualityMismatches(read, minBaseQuality) > 1) {
                addLocation(read, locations);
              }
            }
          }
          
          reader.close();
    }
    
    return locationsToRegions(locations, readLength);
  }

View Full Code Here

  }
  
  @Override
  public void go() throws Exception {
    System.err.println("Starting chromosome: " + chromosome);
    SAMFileReader rdr = new SAMFileReader(new File(inputFile));
    rdr.setValidationStringency(ValidationStringency.SILENT);
    
    Iterator<SAMRecord> iter = rdr.queryContained(chromosome, 0, 0);
      
    while (iter.hasNext()) {
      SAMRecord read = (SAMRecord) iter.next();
      
      outputWriter.addAlignment(read);
    }
    
    rdr.close();
    System.err.println("Chromosome: " + chromosome + " done.");
  }

View Full Code Here

    File dir = new File(outputDirectory);
    if (!dir.exists()) {
      dir.mkdir();
    }
    
    SAMFileReader rdr = new SAMFileReader(new File(filename));
    
    ThreadManager threads = new ThreadManager(numThreads);
    
    Map<String, SAMFileWriter> outputWriterMap = new HashMap<String, SAMFileWriter>();
    
    SAMFileWriterFactory writerFactory = new SAMFileWriterFactory();
    writerFactory.setUseAsyncIo(false);
    
    // Farm each chromosome out to its own thread.
    for (SAMSequenceRecord chr : rdr.getFileHeader().getSequenceDictionary().getSequences()) {    
      SAMFileWriter writer = writerFactory.makeSAMOrBAMWriter(
          rdr.getFileHeader(), false, new File(outputDirectory + "/" + chr.getSequenceName() + ".bam"));
      
      outputWriterMap.put(chr.getSequenceName(), writer);
      
      BamSplitterThread thread = new BamSplitterThread(threads, filename, chr.getSequenceName(), writer);
      threads.spawnThread(thread);
    }
    threads.waitForAllThreadsToComplete();
    
    // Now go back and retrieve the unmapped reads.
    System.err.println("Processing unmapped reads");
    Iterator<SAMRecord> iter = rdr.queryUnmapped();
    while (iter.hasNext()) {
      SAMRecord read = iter.next();
    
      // If this read is not assigned a position, but the mate is, include in the output BAM associated with mate's chromosome.
      if (read.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX && read.getMateReferenceIndex() != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
        SAMFileWriter writer = outputWriterMap.get(read.getMateReferenceName());
        writer.addAlignment(read);
      }
    }
    
    for (SAMFileWriter writer : outputWriterMap.values()) {
      writer.close();
    }
    
    rdr.close();
    
    long e = System.currentTimeMillis();
    
    System.err.println("BAMSplitter done.  Elapsed minutes: " + (double) (e-s)/1000.0/60.0);
  }

View Full Code Here


  public Map<String, Integer> countReadsSupportingBreakpoints(String svSam, int readLength, SAMFileHeader samHeader) {
    
    samStringReader = new SamStringReader(samHeader);
    
    SAMFileReader reader = new SAMFileReader(new File(svSam));
    reader.setValidationStringency(ValidationStringency.SILENT);


    String fullMatch = readLength + "M";
    
    // Require 90% of the read to overlap the breakpoint
    int minStart = (int) (readLength * .10);
    int maxStart = (int) (readLength *.9) + 1;
    
    // TODO: Need way to query mapped reads only
    for (SAMRecord read : reader) {
      if (!read.getReadUnmappedFlag() && read.getCigarString().equals(fullMatch)) {
        if (read.getAlignmentStart() >= minStart && read.getAlignmentStart() <= maxStart) {
          int editDistance = SAMRecordUtils.getIntAttribute(read, "NM");
          
          if (editDistance <= MAX_EDIT_DISTANCE) {
            SAMRecord orig = getOrigRecord(read, samHeader);
            int origEditDistance = SAMRecordUtils.getIntAttribute(orig, "YX");
            if (editDistance < origEditDistance) {
              //TODO: Inspect alternate alignments
              String[] refFields = read.getReferenceName().split("_");
              if (refFields.length >= 6) {
                String breakpointGroupId = refFields[0] + "_" + refFields[1] + "\t" + refFields[2] + ":" + refFields[3] + "\t" +
                    refFields[4] + ":" + refFields[5];
                Integer count = breakpointCounts.get(breakpointGroupId);
                if (count == null) {
                  breakpointCounts.put(breakpointGroupId, 1);
                } else {
                  breakpointCounts.put(breakpointGroupId, count + 1);
                }
              } else {
                System.out.println("Error analyzing breakpoint for: " + read.getSAMString());
              }
            }
          }
        }
      }
    }
    
    reader.close();
    
    return breakpointCounts;
  }

View Full Code Here

  }
  
  public static void main(String[] args) {
    SVReadCounter svc = new SVReadCounter();
    
    SAMFileReader reader = new SAMFileReader(new File("/home/lmose/dev/abra/sv/test_sv.sam"));
    
    Map<String, Integer> counts = svc.countReadsSupportingBreakpoints("/home/lmose/dev/abra/sv/sv.bam", 100, reader.getFileHeader());
    
    for (String bp : counts.keySet()) {
      System.out.println(bp + "\t" + counts.get(bp));
    }
    
    reader.close();
  }

View Full Code Here

public class ReadLocusReader implements Iterable<ReadsAtLocus> {


  private SAMFileReader samReader;
  
  public ReadLocusReader(String samFile) {
        samReader = new SAMFileReader(new File(samFile));
        samReader.setValidationStringency(ValidationStringency.SILENT);
  }

View Full Code Here

    return chunks;
  }
  
  public void chopClopDrop(List<Feature> regions, String input, String output) {
    
    SAMFileReader reader = new SAMFileReader(new File(input));
    
    SAMFileHeader header = reader.getFileHeader();
    header.setSortOrder(SortOrder.unsorted);
    
    SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(
        header, true, new File(output));
    
    for (Feature region : regions) {
      Map<String, SAMRecord> chunks = chop(reader, region);
      chunks = clop(chunks);
      chunks = drop(chunks);
      
      for (SAMRecord chunk : chunks.values()) {
        out.addAlignment(chunk);
      }
    }
    
    reader.close();
    out.close();
  }

View Full Code Here

        readIds = new HashSet<String>();
        List<SAMRecord> reads = new ArrayList<SAMRecord>();
        readsList.add(reads);
        
        for (Feature region : regions) {
          SAMFileReader reader = new SAMFileReader(new File(input));
          reader.setValidationStringency(ValidationStringency.SILENT);
          
          Iterator<SAMRecord> iter;
          if (region != null) {
            iter = reader.queryOverlapping(region.getSeqname(), (int) region.getStart(), (int) region.getEnd());
          } else {
            iter = reader.iterator();
          }
          
          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;
          }
        }

View Full Code Here


//    String in = "/home/lmose/dev/abra/leftalign.sam";
//    String out = "/home/lmose/dev/abra/la.out.sam";
//    String ref = "/home/lmose/reference/chr1/chr1.fa";
    
    SAMFileReader reader = new SAMFileReader(new File(in));
    
    SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(
        reader.getFileHeader(), false, new File(out));
    
    CompareToReference2 c2r = new CompareToReference2();
    c2r.init(ref);
    
    IndelShifter indelShifter = new IndelShifter();


    for (SAMRecord read : reader) {
      SAMRecord shiftedRead = indelShifter.shiftIndelsLeft(read, c2r);
      writer.addAlignment(shiftedRead);
    }
    
    writer.close();
    reader.close();
  }

View Full Code Here

0 1 2 3 4 5

TOP

Related Classes of net.sf.samtools.SAMFileReader

abra.bamsplitter.BamSplitter

abra.bamsplitter.BamSplitterThread

abra.BetaPairValidatingRealignmentWriter

abra.cadabra.ReadLocusReader

abra.cadabra.SomaticLocusCaller

abra.CompareToReference

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.