Examples of SAMFileReader

Examples of net.sf.samtools.SAMFileReader

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
   
    log("Adjusting reads.");
   
    RealignmentWriter writer = getRealignmentWriter(outputSam, isTightAlignment, tempDir);
   
    SAMFileReader contigReader = new SAMFileReader(new File(alignedToContigSam));
    contigReader.setValidationStringency(ValidationStringency.SILENT);
   
    SamStringReader samStringReader = new SamStringReader(samHeader);
   
    for (SAMRecord read : contigReader) {
     
      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());

      SAMRecord readToOutput = null;
     
      // Only adjust reads that align to contig with no indel and shorter edit distance than the original alignment
      String matchingString = read.getReadLength() + "M";
      if ((read.getCigarString().equals(matchingString)) &&
        (read.getReadUnmappedFlag() == false)  &&
        (!orig.getCigarString().contains("N")) &&  // Don't remap introns
        (SAMRecordUtils.getEditDistance(read, null) < SAMRecordUtils.getOrigEditDistance(orig)) &&
        (!isFiltered(orig))) {
       
        SAMRecord origRead = orig;
        String contigReadStr = read.getReferenceName();
       
        int numBestHits = SAMRecordUtils.getIntAttribute(read, "X0");
        int numSubOptimalHits = SAMRecordUtils.getIntAttribute(read, "X1");
       
        int totalHits = numBestHits + numSubOptimalHits;
       
        List<HitInfo> bestHits = getBestHits(contigReadStr, samStringReader, read, matchingString);
       
        Map<String, SAMRecord> outputReadAlignmentInfo = convertBestHitsToAlignmentInfo(bestHits, origRead);
       
        readToOutput = getUpdatedReadInfo(outputReadAlignmentInfo, read,
            orig, origRead, c2r, totalHits, isTightAlignment);
       
      }
     
      adjustForStrand(read.getReadNegativeStrandFlag(), orig);
      writer.addAlignment(readToOutput, orig);
    }

    int realignedCount = writer.flush();
    contigReader.close();
   
    log("Done adjusting reads.  Number of reads realigned: " + realignedCount);
  }
View Full Code Here

Examples of net.sf.samtools.SAMFileReader

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
  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

Examples of net.sf.samtools.SAMFileReader

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
      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

Examples of net.sf.samtools.SAMFileReader

25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
  }
 
  @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

Examples of net.sf.samtools.SAMFileReader

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
    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

Examples of net.sf.samtools.SAMFileReader

19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
  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

Examples of net.sf.samtools.SAMFileReader

86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
  }
 
  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

Examples of net.sf.samtools.SAMFileReader

14
15
16
17
18
19
20
21
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

Examples of net.sf.samtools.SAMFileReader

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    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

Examples of net.sf.samtools.SAMFileReader

238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
        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
TOP
Copyright © 2018 www.massapi.com. All rights reserved.
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.