Java Code Examples for htsjdk.samtools.SamReader#close()

private void check(final File bam) throws IOException {
  final SamReader normalReader = SamUtils.makeSamReader(bam);
  final Iterator<SAMRecord> norm = normalReader.iterator();
  try {
    try (final RecordIterator<SAMRecord> closed = new SamClosedFileReader(bam, null, null, normalReader.getFileHeader())) {
      while (norm.hasNext() && closed.hasNext()) {
        final SAMRecord a = norm.next();
        final SAMRecord b = closed.next();
        assertEquals(a.getSAMString().trim(), b.getSAMString().trim());
      }
      assertEquals(norm.hasNext(), closed.hasNext());
    }
  } finally {
    normalReader.close();
  }
}
 

@Test
public void testReader() throws Exception {
  int expectedCount = 0;
  SamReader samReader = SamReaderFactory.makeDefault().open(new File(input));
  for (SAMRecord r : samReader) {
    expectedCount++;
  }
  samReader.close();

  AnySAMInputFormat inputFormat = new AnySAMInputFormat();
  List<InputSplit> splits = inputFormat.getSplits(jobContext);
  assertEquals(1, splits.size());
  RecordReader<LongWritable, SAMRecordWritable> reader = inputFormat
      .createRecordReader(splits.get(0), taskAttemptContext);
  reader.initialize(splits.get(0), taskAttemptContext);

  int actualCount = 0;
  while (reader.nextKeyValue()) {
    actualCount++;
  }
  reader.close();

  assertEquals(expectedCount, actualCount);
}
 

private void convertFileAndVerifyRecordCount(final int expectedCount,
                   final String fastqFilename1,
                   final String fastqFilename2,
                   final FastqQualityFormat version,
                   final boolean permissiveFormat,
                   final boolean useSequentialFastqs) throws IOException {
    final File samFile = convertFile(fastqFilename1, fastqFilename2, version, permissiveFormat, useSequentialFastqs);
    final SamReader samReader = SamReaderFactory.makeDefault().open(samFile);
    final SAMRecordIterator iterator = samReader.iterator();
    int actualCount = 0;
    while (iterator.hasNext()) {
        iterator.next();
        actualCount++;
    }
    samReader.close();
    Assert.assertEquals(actualCount, expectedCount);
}
 

protected static Map<String, Map<String, MatePair>> createPUPairsMap(final File samFile) throws IOException {
    IOUtil.assertFileIsReadable(samFile);
    final SamReader reader = SamReaderFactory.makeDefault().open(samFile);
    final Map<String, Map<String, MatePair>> map = new LinkedHashMap<>();

    Map<String,MatePair> curFileMap;
    for (final SAMRecord record : reader ) {
        final String platformUnit = record.getReadGroup().getPlatformUnit();
        curFileMap = map.get(platformUnit);
        if(curFileMap == null) {
            curFileMap = new LinkedHashMap<>();
            map.put(platformUnit, curFileMap);
        }

        MatePair mpair = curFileMap.get(record.getReadName());
        if (mpair == null) {
             mpair = new MatePair();
             curFileMap.put(record.getReadName(), mpair);
        }
        mpair.add(record);
    }
    reader.close();
    return map;
}
 

private Map<String,MatePair> createSamMatePairsMap(final File samFile) throws IOException {
    IOUtil.assertFileIsReadable(samFile);
    final SamReader reader = SamReaderFactory.makeDefault().open(samFile);

    final Map<String,MatePair> map = new LinkedHashMap<String,MatePair>();
    for (final SAMRecord record : reader ) {
        MatePair mpair = map.get(record.getReadName());
        if (mpair == null) {
             mpair = new MatePair();
             map.put(record.getReadName(), mpair);
        }
        mpair.add(record);
    }
    reader.close();
    return map;
}
 

public void call(String normal, String tumor, String vcf, String reference, int minBaseQual, double minMaf) throws IOException {
	loadLoci(vcf);
	c2r = new CompareToReference2();
	c2r.init(reference);
	this.minBaseQual = minBaseQual;
	this.minMaf = minMaf;
	
	System.err.println("Processing positions");

	SamReader normalReader = SAMRecordUtils.getSamReader(normal);
	
	SamReader tumorReader = SAMRecordUtils.getSamReader(tumor);
       
	for (LocusInfo locus : loci) {
		locus.normalCounts = getCounts(normalReader, locus);
		locus.tumorCounts = getCounts(tumorReader, locus);
	}
	
	normalReader.close();
	tumorReader.close();
	
	System.err.println("Writing results");
	
	outputResults();
}
 

private static void sliceFromVCF(@NotNull CommandLine cmd) throws IOException {
    String inputPath = cmd.getOptionValue(INPUT);
    String vcfPath = cmd.getOptionValue(VCF);
    int proximity = Integer.parseInt(cmd.getOptionValue(PROXIMITY, "500"));

    SamReaderFactory readerFactory = createFromCommandLine(cmd);
    SamReader reader = readerFactory.open(new File(inputPath));

    QueryInterval[] intervals = getIntervalsFromVCF(vcfPath, reader.getFileHeader(), proximity);
    CloseableIterator<SAMRecord> iterator = reader.queryOverlapping(intervals);
    SAMFileWriter writer = new SAMFileWriterFactory().setCreateIndex(true)
            .makeBAMWriter(reader.getFileHeader(), true, new File(cmd.getOptionValue(OUTPUT)));

    writeToSlice(writer, iterator);

    writer.close();
    reader.close();
}
 

@Test
public void testShardedWriting() throws Exception {
  final String OUTPUT = helper.getTestOutputGcsFolder() + OUTPUT_FNAME;
  String[] ARGS = {
      "--project=" + helper.getTestProject(),
      "--output=" + OUTPUT,
      "--numWorkers=18",
      "--runner=DataflowRunner",
      "--wait=true",
      "--stagingLocation=" + helper.getTestStagingGcsFolder(),
      "--references=" + TEST_CONTIG,
      "--BAMFilePath=" + TEST_BAM_FNAME,
      "--lociPerWritingShard=1000000"
  };
  SamReader reader = null;
  try {
    helper.touchOutput(OUTPUT);

    ShardedBAMWriting.main(ARGS);

    reader = helper.openBAM(OUTPUT);
    Assert.assertTrue(reader.hasIndex());
    final int sequenceIndex = reader.getFileHeader().getSequenceIndex("11");
    BAMIndexMetaData metaData = reader.indexing().getIndex().getMetaData(sequenceIndex);
    Assert.assertEquals(EXPECTED_ALL_READS - EXPECTED_UNMAPPED_READS,
        metaData.getAlignedRecordCount());
    // Not handling unmapped reads yet
    // Assert.assertEquals(EXPECTED_UNMAPPED_READS,
    // metaData.getUnalignedRecordCount());

  } finally {
    if (reader != null) {
      reader.close();
    }
    helper.deleteOutput(OUTPUT);
  }
}
 

@Test
public void testMapReduceJob() throws Exception {
  Configuration conf = new Configuration();

  FileSystem fileSystem = FileSystem.get(conf);
  Path inputPath = new Path(input);
  Path outputPath = fileSystem.makeQualified(new Path("target/out"));
  fileSystem.delete(outputPath, true);

  Job job = Job.getInstance(conf);
  FileInputFormat.setInputPaths(job, inputPath);
  job.setInputFormatClass(SAMInputFormat.class);
  job.setOutputKeyClass(LongWritable.class);
  job.setOutputValueClass(SAMRecordWritable.class);
  job.setNumReduceTasks(0);
  FileOutputFormat.setOutputPath(job, outputPath);

  boolean success = job.waitForCompletion(true);
  assertTrue(success);

  List<String> samStrings = new ArrayList<String>();
  SamReader samReader = SamReaderFactory.makeDefault().open(new File(input));
  for (SAMRecord r : samReader) {
    samStrings.add(r.getSAMString().trim());
  }
  samReader.close();

  File outputFile = new File(new File(outputPath.toUri()), "part-m-00000");
  BufferedReader br = new BufferedReader(new FileReader(outputFile));
  String line;
  int index = 0;
  while ((line = br.readLine()) != null) {
    String value = line.substring(line.indexOf("\t") + 1); // ignore key
    assertEquals(samStrings.get(index++), value);
  }
  br.close();
}
 

/**
 * Test that clipping of FR reads for fragments shorter than read length happens only when it should.
 */
@Test
public void testShortFragmentClipping() throws Exception {
    final File output = File.createTempFile("testShortFragmentClipping", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "cliptest.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "cliptest.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "cliptest.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    final Map<String, SAMRecord> firstReadEncountered = new HashMap<String, SAMRecord>();

    for (final SAMRecord rec : result) {
        final SAMRecord otherEnd = firstReadEncountered.get(rec.getReadName());
        if (otherEnd == null) {
            firstReadEncountered.put(rec.getReadName(), rec);
        } else {
            final int fragmentStart = Math.min(rec.getAlignmentStart(), otherEnd.getAlignmentStart());
            final int fragmentEnd = Math.max(rec.getAlignmentEnd(), otherEnd.getAlignmentEnd());
            final String[] readNameFields = rec.getReadName().split(":");
            // Read name of each pair includes the expected fragment start and fragment end positions.
            final int expectedFragmentStart = Integer.parseInt(readNameFields[1]);
            final int expectedFragmentEnd = Integer.parseInt(readNameFields[2]);
            Assert.assertEquals(fragmentStart, expectedFragmentStart, rec.getReadName());
            Assert.assertEquals(fragmentEnd, expectedFragmentEnd, rec.getReadName());
        }
    }
    result.close();
}
 

public static HeaderInfo getHeaderFromBAMFile(Storage.Objects storage, String BAMPath, Iterable<Contig> explicitlyRequestedContigs) throws IOException {
  HeaderInfo result = null;

  // Open and read start of BAM
  LOG.info("Reading header from " + BAMPath);
  final SamReader samReader = BAMIO
      .openBAM(storage, BAMPath, ValidationStringency.DEFAULT_STRINGENCY);
  final SAMFileHeader header = samReader.getFileHeader();
  Contig firstContig = getFirstExplicitContigOrNull(header, explicitlyRequestedContigs);
  if (firstContig == null) {
    final SAMSequenceRecord seqRecord = header.getSequence(0);
    firstContig = new Contig(seqRecord.getSequenceName(), -1, -1);
  }

  LOG.info("Reading first chunk of reads from " + BAMPath);
  final SAMRecordIterator recordIterator = samReader.query(
      firstContig.referenceName, (int)firstContig.start + 1, (int)firstContig.end + 1, false);

  Contig firstShard = null;
  while (recordIterator.hasNext() && result == null) {
    SAMRecord record = recordIterator.next();
    final int alignmentStart = record.getAlignmentStart();
    if (firstShard == null && alignmentStart > firstContig.start &&
        (alignmentStart < firstContig.end || firstContig.end == -1)) {
      firstShard = new Contig(firstContig.referenceName, alignmentStart, alignmentStart);
      LOG.info("Determined first shard to be " + firstShard);
      result = new HeaderInfo(header, firstShard);
    }
  }
  recordIterator.close();
  samReader.close();

  if (result == null) {
    throw new IOException("Did not find reads for the first contig " + firstContig.toString());
  }
  LOG.info("Finished header reading from " + BAMPath);
  return result;
}
 

@Test
public void testRemoveNmMdAndUqOnOverlappingReads() throws IOException {
    final File output = File.createTempFile("testRemoveNmMdAndUqOnOverlappingReads", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "removetags.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "removetags.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "removetags.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, SAMFileHeader.SortOrder.queryname);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    for (final SAMRecord rec : result) {
        boolean hasTags = false;
        if (rec.getReadName().startsWith("CLIPPED")) {
            final String[] readNameFields = rec.getReadName().split(":");
            final int index = rec.getFirstOfPairFlag() ? 1 : 2;
            hasTags = Integer.parseInt(readNameFields[index]) == 1;
        }
        if (hasTags) {
            Assert.assertNull(rec.getAttribute("MD"));
            Assert.assertNull(rec.getAttribute("NM"));
        } else {
            Assert.assertNotNull(rec.getAttribute("MD"));
            Assert.assertNotNull(rec.getAttribute("NM"));
        }
    }
    result.close();
}
 

private static void sliceFromURLs(@NotNull URL indexUrl, @NotNull URL bamUrl, @NotNull CommandLine cmd) throws IOException {
    File indexFile = downloadIndex(indexUrl);
    indexFile.deleteOnExit();

    SamReader reader = createFromCommandLine(cmd).open(SamInputResource.of(bamUrl).index(indexFile));

    BAMIndex bamIndex;
    if (indexFile.getPath().contains(".crai")) {
        SeekableStream craiIndex = CRAIIndex.openCraiFileAsBaiStream(indexFile, reader.getFileHeader().getSequenceDictionary());
        bamIndex = new DiskBasedBAMFileIndex(craiIndex, reader.getFileHeader().getSequenceDictionary());
    } else {
        bamIndex = new DiskBasedBAMFileIndex(indexFile, reader.getFileHeader().getSequenceDictionary(), false);
    }

    Optional<Pair<QueryInterval[], BAMFileSpan>> queryIntervalsAndSpan = queryIntervalsAndSpan(reader, bamIndex, cmd);
    Optional<Chunk> unmappedChunk = getUnmappedChunk(bamIndex, HttpUtils.getHeaderField(bamUrl, "Content-Length"), cmd);
    List<Chunk> sliceChunks = sliceChunks(queryIntervalsAndSpan, unmappedChunk);
    SamReader cachingReader = createCachingReader(indexFile, bamUrl, cmd, sliceChunks);

    SAMFileWriter writer = new SAMFileWriterFactory().setCreateIndex(true)
            .makeBAMWriter(reader.getFileHeader(), true, new File(cmd.getOptionValue(OUTPUT)));

    queryIntervalsAndSpan.ifPresent(pair -> {
        LOGGER.info("Slicing bam on bed regions...");
        CloseableIterator<SAMRecord> bedIterator = getIterator(cachingReader, pair.getKey(), pair.getValue().toCoordinateArray());
        writeToSlice(writer, bedIterator);
        LOGGER.info("Done writing bed slices.");
    });

    unmappedChunk.ifPresent(chunk -> {
        LOGGER.info("Slicing unmapped reads...");
        CloseableIterator<SAMRecord> unmappedIterator = cachingReader.queryUnmapped();
        writeToSlice(writer, unmappedIterator);
        LOGGER.info("Done writing unmapped reads.");
    });

    reader.close();
    writer.close();
    cachingReader.close();
}
 

public void runDiff() throws Exception {
  SamReaderFactory readerFactory =  SamReaderFactory
      .makeDefault()
      .validationStringency(ValidationStringency.SILENT)
      .enable(SamReaderFactory.Option.CACHE_FILE_BASED_INDEXES);
  LOG.info("Opening file 1 for diff: " + BAMFile1);
  SamReader reader1 = readerFactory.open(new File(BAMFile1));
  LOG.info("Opening file 2 for diff: " + BAMFile2);
  SamReader reader2 = readerFactory.open(new File(BAMFile2));


  try {
    Iterator<Contig> contigsToProcess = null;
    if (options.contigsToProcess != null && !options.contigsToProcess.isEmpty()) {
      Iterable<Contig> parsedContigs = Contig.parseContigsFromCommandLine(options.contigsToProcess);
      referencesToProcess = Sets.newHashSet();
      for (Contig c : parsedContigs) {
        referencesToProcess.add(c.referenceName);
      }
      contigsToProcess = parsedContigs.iterator();
      if (!contigsToProcess.hasNext()) {
        return;
      }
    }
    LOG.info("Comparing headers");
    if (!compareHeaders(reader1.getFileHeader(), reader2.getFileHeader())) {
      error("Headers are not equal");
      return;
    }
    LOG.info("Headers are equal");
    do {
      SAMRecordIterator it1;
      SAMRecordIterator it2;
      if (contigsToProcess == null) {
        LOG.info("Checking all the reads");
        it1 = reader1.iterator();
        it2 = reader2.iterator();
      } else {
        Contig contig = contigsToProcess.next();
        LOG.info("Checking contig " + contig.toString());
        processedContigs++;
        it1 = reader1.queryOverlapping(contig.referenceName, (int)contig.start,  (int)contig.end);
        it2 = reader2.queryOverlapping(contig.referenceName, (int)contig.start,  (int)contig.end);
      }

      if (!compareRecords(it1, it2)) {
        break;
      }

      it1.close();
      it2.close();
    } while (contigsToProcess != null && contigsToProcess.hasNext());
  } catch (Exception ex) {
    throw ex;
  } finally {
    reader1.close();
    reader2.close();
  }
  LOG.info("Processed " + processedContigs + " contigs, " +
      processedLoci + " loci, " + processedReads + " reads.");
}
 

private void getSamHeaderAndReadLength() throws IOException {
	
	Logger.info("Identifying header and determining read length");
	
	this.samHeaders = new SAMFileHeader[this.inputSams.length];
	
	for (int i=0; i<this.inputSams.length; i++) {
	
		SamReader reader = SAMRecordUtils.getSamReader(inputSams[i]);
		try {
			samHeaders[i] = reader.getFileHeader();
			samHeaders[i].setSortOrder(SAMFileHeader.SortOrder.unsorted);
			
			Iterator<SAMRecord> iter = reader.iterator();
			
			int cnt = 0;
			while ((iter.hasNext()) && (cnt < 1000000)) {
				SAMRecord read = iter.next();
				this.readLength = Math.max(this.readLength, read.getReadLength());
				this.maxMapq = Math.max(this.maxMapq, read.getMappingQuality());
				
				// Assumes aligner sets proper pair flag correctly
				if ((isPairedEnd) && (read.getReadPairedFlag()) && (read.getProperPairFlag())) {
					this.minInsertLength = Math.min(this.minInsertLength, Math.abs(read.getInferredInsertSize()));
					this.maxInsertLength = Math.max(this.maxInsertLength, Math.abs(read.getInferredInsertSize()));
				}
				
				cnt += 1;
			}
			
			// Allow some fudge in insert length
			minInsertLength = Math.max(minInsertLength - 2*readLength, 0);
			maxInsertLength = maxInsertLength + 2*readLength;
			
		} finally {
			reader.close();
		}
	}
	
	Logger.info("Min insert length: " + minInsertLength);
	Logger.info("Max insert length: " + maxInsertLength);
			
	Logger.info("Max read length is: " + readLength);
	if (assemblerSettings.getMinContigLength() < 1) {
		assemblerSettings.setMinContigLength(Math.max(readLength+1, MIN_CONTIG_LENGTH));
	}
	Logger.info("Min contig length: " + assemblerSettings.getMinContigLength());
}
 

public void loadReads(File[] bams) throws IOException{

int count = 0;
int numOp = 0;

for(File bam : bams){
    HLA.log.appendln("Loading reads from:\t" + bam.getName());
    Object2IntOpenHashMap<String> readLoadingSet = new Object2IntOpenHashMap<String>();
    readLoadingSet.defaultReturnValue(0);
    
    final SamReader reader = SamReaderFactory.makeDefault().open(bam);
    
    //Kourami bam checker added
    if(!checkHeader(reader.getFileHeader())){
	HLA.log.appendln("Unexpected BAM :\t"+ bam.getName() 
			+"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			+ "Please use the recommended preprocessing steps explained on the github page:\n"
			+ "https://github.com/Kingsford-Group/kourami");
	System.err.println("Unexpected BAM :\t"+ bam.getName() 
			   +"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			   + "Please use the recommended preprocessing steps explained on the github page:\n"
			   + "https://github.com/Kingsford-Group/kourami");
	HLA.log.outToFile();
	System.exit(1);
    }

    for(final SAMRecord samRecord : reader){
	if(count == 0){
	    HLA.READ_LENGTH = samRecord.getReadLength();
	    HLA.log.appendln("Setting HLA.READ_LEGNTH = " + HLA.READ_LENGTH);
	}
	//added checking to process reads matching to HLA-type sequences
	//discarding decoy hits (DQB2, DQA2)
	boolean qc = false;
	if( (samRecord.getReferenceName().indexOf("*") > -1) 
	    && !samRecord.getReadUnmappedFlag() 
	    && !samRecord.isSecondaryOrSupplementary() 
	    && !this.startWIns(samRecord)){
	    count++;
	    if(samRecord.getReadPairedFlag())
		numOp += processRecord(samRecord, readLoadingSet);
	    else
		numOp += processRecordUnpaired(samRecord);
	}
	
	if(HLA.DEBUG && count%10000 == 0)
	    HLA.log.appendln("Processed 10000 reads...");
    }
    reader.close();
}
HLA.log.appendln("Loaded a total of " + count + " mapped reads.");
HLA.log.appendln("A total of " + numOp + " bases");
   }
 

/**
 * Run IlluminaBasecallsToSam on a few test cases, and verify that results agree with hand-checked expectation.
 */
@Test(dataProvider="data")
public void testBasic(final String LANE, final String readStructure,
                      final String fivePrimerAdapter, final String threePrimerAdapter,
                      final int expectedNumClippedRecords) throws Exception {
    // Create the SAM file from Gerald output
    final File samFile = File.createTempFile("." + LANE + ".illuminaBasecallsToSam", ".sam");
    samFile.deleteOnExit();
    final List<String> illuminaArgv = CollectionUtil.makeList("BASECALLS_DIR=" + TEST_DATA_DIR,
            "LANE=" + LANE,
            "RUN_BARCODE=" + RUN_BARCODE,
            "READ_STRUCTURE=" + readStructure,
            "OUTPUT=" + samFile,
            "SEQUENCING_CENTER=BI",
            "ALIAS=" + ALIAS
    );
    if (fivePrimerAdapter != null) {
        illuminaArgv.addAll(CollectionUtil.makeList(
            "ADAPTERS_TO_CHECK=null",
            "FIVE_PRIME_ADAPTER=" + fivePrimerAdapter,
            "THREE_PRIME_ADAPTER=" + threePrimerAdapter
        ));
    }
    Assert.assertEquals(runPicardCommandLine(illuminaArgv), 0);

    // Read the file and confirm it contains what is expected
    final SamReader samReader = SamReaderFactory.makeDefault().open(samFile);

    // look for clipped adaptor attribute in lane 3 PE (2) and in lane 6 (1) non-PE
    int count = 0;
    for (final SAMRecord record : samReader) {
        if (record.getIntegerAttribute(ReservedTagConstants.XT) != null) {
            count++;
            if ((count == 1 || count == 2) && LANE.equals("2")){
                Assert.assertEquals (114, (int)record.getIntegerAttribute(ReservedTagConstants.XT));
            } else if (count == 1 || count == 2 && LANE.equals("1")) {
                Assert.assertEquals(68, (int) record.getIntegerAttribute(ReservedTagConstants.XT));
            }
        }
    }

    // Check the total number of clipped records
    Assert.assertEquals(count, expectedNumClippedRecords);

    samReader.close();
    samFile.delete();
}
 

public void countSplices(String input, Set<SpliceJunction> annotatedJunctions) throws IOException {
	
	
	SamReader reader = SAMRecordUtils.getSamReader(input);
	
	updateCounts(reader);
	
	List<SpliceJunction> junctions = new ArrayList<SpliceJunction>(uniqueReads.keySet());
	Collections.sort(junctions, new SpliceJunctionComparator(reader.getFileHeader()));

	reader.close();
	
	for (SpliceJunction junction : junctions) {
		
		int annotated = annotatedJunctions.contains(junction) ? 1 : 0;
		
		String rec = String.format("%s\t%d\t%d\t.\t.\t%d\t%d\t%d\t.", junction.chrom, junction.start, junction.stop,
				annotated, uniqueReads.get(junction), multiMapReads.get(junction));
		
		System.out.println(rec);
	}
}
 

private void processUnmapped(SAMFileWriter output, String inputBam) throws IOException {
	
	Logger.debug("Processing unmapped reads...");
	
	SamReader reader = SAMRecordUtils.getSamReader(inputBam);

	// This should give us only read pairs with both ends unmapped
	Iterator<SAMRecord> iter = reader.queryUnmapped();
	
	while (iter.hasNext()) {
		SAMRecord read = iter.next();
		output.addAlignment(read);
	}
	
	reader.close();
}
 

public static void main(String[] args) throws IOException {
	Params params = new Params();
	JCommander jc = new JCommander(params);
	jc.parse(args);

	Log.setGlobalLogLevel(params.logLevel);

	if (args.length == 0 || params.help) {
		usage(jc);
		System.exit(1);
	}

	if (params.reference == null) {
		System.out.println("Reference file not found, will try downloading...");
	}

	ReferenceSource referenceSource = null;
	if (params.reference != null) {
		System.setProperty("reference", params.reference.getAbsolutePath());
		referenceSource = new ReferenceSource(params.reference);
	} else {
		String prop = System.getProperty("reference");
		if (prop != null) {
			referenceSource = new ReferenceSource(new File(prop));
		}
	}

	SamReaderFactory factory = SamReaderFactory.make().validationStringency(params.validationLevel);
	SamInputResource r;
	if ("file".equalsIgnoreCase(params.url.getProtocol()))
		r = SamInputResource.of(params.url.getPath());
	else
		r = SamInputResource.of(params.url);
	SamReader reader = factory.open(r);
	SAMRecordIterator iterator = reader.iterator();

	SAMFileWriterFactory writerFactory = new SAMFileWriterFactory();
	SAMFileWriter writer = null;
	OutputStream os = new BufferedOutputStream(new FileOutputStream(params.outputFile));
	switch (params.outputFormat) {
	case BAM:
		writer = writerFactory.makeBAMWriter(reader.getFileHeader(),
				reader.getFileHeader().getSortOrder() == SortOrder.coordinate, os);
		break;
	case CRAM:
		writer = writerFactory.makeCRAMWriter(reader.getFileHeader(), os, params.reference);
		break;

	default:
		System.out.println("Unknown output format: " + params.outputFormat);
		System.exit(1);
	}

	while (iterator.hasNext()) {
		writer.addAlignment(iterator.next());
	}
	writer.close();
	reader.close();
}