Java Code Examples for htsjdk.samtools.util.CloseableIterator#close()

@Test(dataProvider = "goodInputVcfs")
public void testJavaScript(final File input) throws Exception {
    final File out = VcfTestUtils.createTemporaryIndexedFile("filterVcfTestJS.", ".vcf");
    final FilterVcf filterer = new FilterVcf();
    filterer.INPUT = input;
    filterer.OUTPUT = out;
    filterer.JAVASCRIPT_FILE = quickJavascriptFilter("variant.getStart()%5 != 0");

    final int retval = filterer.doWork();
    Assert.assertEquals(retval, 0);

    //count the number of reads
    final int expectedNumber = 4;
    int count = 0;
    VCFFileReader in = new VCFFileReader(filterer.OUTPUT, false);
    CloseableIterator<VariantContext> iter = in.iterator();
    while (iter.hasNext()) {
        final VariantContext ctx = iter.next();
        count += (ctx.isFiltered() ? 1 : 0);
    }
    iter.close();
    in.close();
    Assert.assertEquals(count, expectedNumber);
}
 

/**
 * Make sure that the order of the output file is identical to the order
 * of the input files by iterating through the output, making sure that,
 * if the context is an indel (snp), the next genomic position in the indel
 * (snp) queue is the same. Also make sure that the context is in the order
 * specified by the input files.
 */
private void validateSnpAndIndelResults(final File output, final Queue<String> indelContigPositions, final Queue<String> snpContigPositions) {
    final VCFFileReader outputReader = new VCFFileReader(output, false);
    final VariantContextComparator outputComparator = outputReader.getFileHeader().getVCFRecordComparator();
    VariantContext last = null;
    final CloseableIterator<VariantContext> iterator = outputReader.iterator();
    while (iterator.hasNext()) {
        final VariantContext outputContext = iterator.next();
        if (outputContext.isIndel()) Assert.assertEquals(getContigPosition(outputContext), indelContigPositions.poll());
        if (outputContext.isSNP()) Assert.assertEquals(getContigPosition(outputContext), snpContigPositions.poll());
        if (last != null) Assert.assertTrue(outputComparator.compare(last, outputContext) <= 0);
        last = outputContext;
    }
    iterator.close();

    // We should have polled everything off the indel (snp) queues
    Assert.assertEquals(indelContigPositions.size(), 0);
    Assert.assertEquals(snpContigPositions.size(), 0);
}
 

/**
 * Checks the ordering and total number of variant context entries in the specified output VCF file.
 * Does NOT check explicitly that the VC genomic positions match exactly those from the inputs. We assume this behavior from other tests.
 *
 * @param output VCF file representing the output of SortVCF
 * @param expectedVariantContextCount the total number of variant context entries from all input files that were merged/sorted
 */
private void validateSortingResults(final File output, final int expectedVariantContextCount) {
    final VCFFileReader outputReader = new VCFFileReader(output, false);
    final VariantContextComparator outputComparator = outputReader.getFileHeader().getVCFRecordComparator();
    VariantContext last = null;
    int variantContextCount = 0;
    final CloseableIterator<VariantContext> iterator = outputReader.iterator();
    while (iterator.hasNext()) {
        final VariantContext outputContext = iterator.next();
        if (last != null) Assert.assertTrue(outputComparator.compare(last, outputContext) <= 0);
        last = outputContext;
        variantContextCount++;
    }
    iterator.close();
    Assert.assertEquals(variantContextCount, expectedVariantContextCount);
}
 

@Override
public void run() {
    try {
        Optional<CloseableIterator<VariantContext>> readerMaybe;
        while ((readerMaybe = vcIterators.next()).isPresent()) {
            final CloseableIterator<VariantContext> reader = readerMaybe.get();
            while (reader.hasNext()) processor.accumulate(reader.next());
            reader.close();

            if (!childrenErrors.isEmpty()) {
                LOG.error(Thread.currentThread() + " aborting: observed error in another child thread.");
                break;
            }
        }
    } catch (final Throwable e) {
        childrenErrors.add(e);
        LOG.error(e, "Unexpected exception encountered in child thread.");
    } finally {
        LOG.debug(String.format("Thread %s is finishing.", Thread.currentThread()));
    }
}
 

private void addVariantsFromFileToBuilder(final PonBuilder ponBuilder, final SAMSequenceRecord samSequenceRecord, final Path file) {
    try (VCFFileReader fileReader = new VCFFileReader(file.toFile(), true)) {
        CloseableIterator<VariantContext> iter =
                fileReader.query(samSequenceRecord.getSequenceName(), 1, samSequenceRecord.getSequenceLength());
        while (iter.hasNext()) {
            ponBuilder.add(iter.next());
        }
        iter.close();
    }
}
 

/**
 * Close any previously-opened iterations over our readers (htsjdk allows only one open iteration per reader).
 */
private void closePreviousIterationsIfNecessary() {
    for ( Map.Entry<SamReader, CloseableIterator<SAMRecord>> readerEntry : readers.entrySet() ) {
        CloseableIterator<SAMRecord> readerIterator = readerEntry.getValue();
        if ( readerIterator != null ) {
            readerIterator.close();
            readerEntry.setValue(null);
        }
    }
}
 

static Queue<String> loadContigPositions(final File inputFile) {
	final VCFFileReader reader = new VCFFileReader(inputFile, false);
	final Queue<String> contigPositions = new LinkedList<String>();
	final CloseableIterator<VariantContext> iterator = reader.iterator();
	while (iterator.hasNext()) contigPositions.add(getContigPosition(iterator.next()));
	iterator.close();
	reader.close();
	return contigPositions;
}
 

@Test
public void testPathologicalOrderingAtTheSamePosition() {
    final AbstractMarkDuplicatesCommandLineProgramTester tester = getTester();
    tester.getSamRecordSetBuilder().setReadLength(68);

    tester.setExpectedOpticalDuplicate(1);

    tester.addMatePair("RUNID:3:1:15013:113051", 0, 129384554, 129384554, false, false, false, false, "68M", "68M", false, false, false, false, false, DEFAULT_BASE_QUALITY);
    tester.addMatePair("RUNID:3:1:15029:113060", 0, 129384554, 129384554, false, false, true, true, "68M", "68M", false, false, false, false, false, DEFAULT_BASE_QUALITY);

    // Create the pathology
    final CloseableIterator<SAMRecord> iterator = tester.getRecordIterator();
    final int[] qualityOffset = {20, 30, 10, 40}; // creates an interesting pathological ordering
    int index = 0;
    while (iterator.hasNext()) {
        final SAMRecord record = iterator.next();
        final byte[] quals = new byte[record.getReadLength()];
        for (int i = 0; i < record.getReadLength(); i++) {
            quals[i] = (byte) (qualityOffset[index] + 10);
        }
        record.setBaseQualities(quals);
        index++;
    }
    iterator.close();

    // Run the test
    tester.runTest();
}
 

/**
 * Do the work after command line has been parsed.
 * RuntimeException may be thrown by this method, and are reported appropriately.
 *
 * @return program exit status.
 */
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    final SamReaderFactory factory = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE);
    if (VALIDATION_STRINGENCY == ValidationStringency.STRICT) {
        factory.validationStringency(ValidationStringency.LENIENT);
    }
    final SamReader reader = factory.open(INPUT);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(reader.getFileHeader(), true, OUTPUT);
    final CloseableIterator<SAMRecord> it = reader.iterator();
    final ProgressLogger progress = new ProgressLogger(Log.getInstance(CleanSam.class));

    // If the read (or its mate) maps off the end of the alignment, clip it
    while (it.hasNext()) {
        final SAMRecord rec = it.next();

        // If the read (or its mate) maps off the end of the alignment, clip it
        AbstractAlignmentMerger.createNewCigarsIfMapsOffEndOfReference(rec);

        // check the read's mapping quality
        if (rec.getReadUnmappedFlag() && 0 != rec.getMappingQuality()) {
            rec.setMappingQuality(0);
        }

        writer.addAlignment(rec);
        progress.record(rec);
    }

    writer.close();
    it.close();
    CloserUtil.close(reader);
    return 0;
}
 

private void closeCurrent() {
  if (mCurrentIt != null) {
    final CloseableIterator<SAMRecord> currentIt = mCurrentIt;
    mCurrentIt = null;
    currentIt.close();
  }
}
 

public void testTabixIndexCreationChr1() throws IOException {

		String configFolder = "test-data/configs/hapmap-chr1";
		// input folder contains no vcf or vcf.gz files
		String inputFolder = "test-data/data/single";

		// create workflow context
		WorkflowTestContext context = buildContext(inputFolder, "hapmap2");
		context.setInput("phasing", "eagle");

		// create step instance
		InputValidation inputValidation = new InputValidationMock(configFolder);

		// run and test
		boolean result = run(context, inputValidation);

		// check if step is failed
		assertEquals(true, result);
		assertTrue(context.hasInMemory("[OK] 1 valid VCF file(s) found."));

		
		// test tabix index and count snps
		String vcfFilename = inputFolder + "/minimac_test.50.vcf.gz";
		VCFFileReader vcfReader = new VCFFileReader(new File(vcfFilename),
				new File(vcfFilename + TabixUtils.STANDARD_INDEX_EXTENSION), true);
		CloseableIterator<VariantContext> snps = vcfReader.query("1", 1, 1000000000);
		int count = 0;
		while (snps.hasNext()) {
			snps.next();
			count++;
		}
		snps.close();
		vcfReader.close();
		
		//check snps
		assertEquals(905, count);

	}
 

public void testTabixIndexCreationChr20() throws IOException {

		String configFolder = "test-data/configs/hapmap-chr1";
		// input folder contains no vcf or vcf.gz files
		String inputFolder = "test-data/data/chr20-phased";

		// create workflow context
		WorkflowTestContext context = buildContext(inputFolder, "hapmap2");

		// create step instance
		InputValidation inputValidation = new InputValidationMock(configFolder);

		// run and test
		boolean result = run(context, inputValidation);

		// check if step is failed
		assertEquals(true, result);
		assertTrue(context.hasInMemory("[OK] 1 valid VCF file(s) found."));

		
		// test tabix index and count snps
		String vcfFilename = inputFolder + "/chr20.R50.merged.1.330k.recode.small.vcf.gz";
		VCFFileReader vcfReader = new VCFFileReader(new File(vcfFilename),
				new File(vcfFilename + TabixUtils.STANDARD_INDEX_EXTENSION), true);
		CloseableIterator<VariantContext> snps = vcfReader.query("20", 1, 1000000000);
		int count = 0;
		while (snps.hasNext()) {
			snps.next();
			count++;
		}
		snps.close();
		vcfReader.close();
		
		//check snps
		assertEquals(7824, count);

	}
 

private static void writeToSlice(@NotNull SAMFileWriter writer, @NotNull CloseableIterator<SAMRecord> iterator) {
    String contig = "";
    while (iterator.hasNext()) {
        SAMRecord record = iterator.next();
        if (record.getContig() != null && !contig.equals(record.getContig())) {
            contig = record.getContig();
            LOGGER.info("Reading contig: {}", contig);
        }
        writer.addAlignment(record);
    }
    iterator.close();
}
 

private Counts getCounts(SamReader reader, LocusInfo locus) {
		
		int depth = 0;
		int altCount = 0;
		int refCount = 0;
		
		CloseableIterator<SAMRecord> iter = reader.queryOverlapping(locus.chromosome, locus.posStart, locus.posStop);
		while (iter.hasNext()) {
			SAMRecord read = iter.next();
			if (!read.getDuplicateReadFlag()) {
				depth += 1;
				
				Object[] baseAndQual = getBaseAndQualAtPosition(read, locus.posStart);
				Character base = (Character) baseAndQual[0];
				int baseQual = (Integer) baseAndQual[1];
				Character refBase = c2r.getSequence(locus.chromosome, locus.posStart, 1).charAt(0);
				
				// Override input with actual reference
//				locus.ref = new String(new char[] { refBase });
				locus.actualRef = new String(new char[] { refBase });
				
				if (locus.isIndel()) {
					if (hasIndel(read, locus)) {
						altCount += 1;
					} else if (!base.equals('N') && base.equals(refBase)) {
						refCount += 1;
					}
				} else if (baseQual >= minBaseQual) {
					
					if (!base.equals('N') && base.equals(refBase)) {
						refCount += 1;
					} else if (base.equals(locus.alt.charAt(0))) {
						altCount += 1;
					}
				}
			}
		}
		
		iter.close();
		
		return new Counts(refCount, altCount, depth);
	}
 

/**
 * Read the data and create a list of lengths
 */
private void read(){
	int counter = 0;
	SAMFileReader reader = new SAMFileReader(bam,index);
	ArrayList<Double> temp = new ArrayList<Double>();
	for (int i = 0; i < genome.size();i++){
		String chr = genome.get(i).getChrom();
		int start = genome.get(i).getStart();
		int stop = genome.get(i).getStop();
		CloseableIterator<SAMRecord> iter = reader.query(chr,start,stop,false);
		while (iter.hasNext()){
			SAMRecord record = null;
			try{
				record = iter.next();
			}
			catch(SAMFormatException ex){
				System.out.println("SAM Record is problematic. Has mapQ != 0 for unmapped read. Will continue anyway");
			}
			if(record != null){
			if(!record.getReadUnmappedFlag() && !record.getMateUnmappedFlag() && record.getFirstOfPairFlag()) {
				if (record.getMappingQuality() >= minQ){
					
					if (Math.abs(record.getInferredInsertSize()) > 100 && Math.abs(record.getInferredInsertSize())
							< 1000){
						counter+=1;
						temp.add((double)Math.abs(record.getInferredInsertSize()));
					}
				}
			}
			}
		}
		iter.close();
	}
	reader.close();
	lengths = new double[counter];
	for (int i = 0;i < temp.size();i++){
		if (temp.get(i) > 100){
			lengths[i] = temp.get(i);
		}
	}
	
}
 

/**
 * Find reads in a given range.
 * 
 * <p>
 * TODO add cigar to the list of returned values
 * <p>
 * TODO add pair information to the list of returned values
 * 
 * @param request
 * @return
 * @throws InterruptedException 
 */
@Override
protected void processDataRequest(DataRequest request) throws InterruptedException {
	
	// Read the given region
	CloseableIterator<SAMRecord> iterator = dataSource.query(request.start.chr, request.start.bp.intValue(), request.end.bp.intValue());
	
	// Produce results
	while (iterator.hasNext()) {

		List<Feature> responseList = new LinkedList<Feature>();

		// Split results into chunks
		for (int c = 0; c < RESULT_CHUNK_SIZE && iterator.hasNext(); c++) {
			SAMRecord record = iterator.next();

			// Region for this read
			Region recordRegion = new Region((long) record.getAlignmentStart(), (long) record.getAlignmentEnd(), request.start.chr);

			// Values for this read
			LinkedHashMap<DataType, Object> values = new LinkedHashMap<DataType, Object>();

			Feature read = new Feature(recordRegion, values);

			if (request.getRequestedContents().contains(DataType.ID)) {
				values.put(DataType.ID, record.getReadName());
			}

			if (request.getRequestedContents().contains(DataType.STRAND)) {
				values.put(DataType.STRAND, BamUtils.getStrand(record, coverageType));					
			}

			if (request.getRequestedContents().contains(DataType.QUALITY)) {
				values.put(DataType.QUALITY, record.getBaseQualityString());
			}

			if (request.getRequestedContents().contains(DataType.CIGAR)) {
				Cigar cigar = new Cigar(read, record.getCigar());
				values.put(DataType.CIGAR, cigar);
			}

			// TODO Deal with "=" and "N" in read string
			if (request.getRequestedContents().contains(DataType.SEQUENCE)) {
				String seq = record.getReadString();
				values.put(DataType.SEQUENCE, seq);
			}

			if (request.getRequestedContents().contains(DataType.MATE_POSITION)) {
				
				BpCoord mate = new BpCoord((Long)(long)record.getMateAlignmentStart(),
						new Chromosome(record.getMateReferenceName()));
				
				values.put(DataType.MATE_POSITION, mate);
			}
			
			if (request.getRequestedContents().contains(DataType.BAM_TAG_NH)) {
				Object ng = record.getAttribute("NH");
				if (ng != null) {
					values.put(DataType.BAM_TAG_NH, (Integer)record.getAttribute("NH"));
				}
			}
			
			/*
			 * NOTE! RegionContents created from the same read area has to be equal in methods equals, hash and compareTo. Primary types
			 * should be ok, but objects (including tables) has to be handled in those methods separately. Otherwise tracks keep adding
			 * the same reads to their read sets again and again.
			 */
			responseList.add(read);
		}

		// Send result			
		super.createDataResult(new DataResult(request.getStatus(), responseList));			
	}

	// We are done
	iterator.close();
}
 

private void calculateCoverage(DataRequest request, BpCoord from, BpCoord to) throws InterruptedException {	
			
	//query data for full average bins, because merging them later would be difficult
	long start = CoverageTool.getBin(from.bp);		
	long end = CoverageTool.getBin(to.bp) + CoverageTool.BIN_SIZE - 1;

	//if end coordinate is smaller than 1 Picard returns a full chromosome and we'll run out of memory
	if (end < 1) {
		end = 1;
	}
			
	CloseableIterator<SAMRecord> iterator = dataSource.query(from.chr, (int)start, (int)end);		
	
	BaseStorage forwardBaseStorage = new BaseStorage();
	BaseStorage reverseBaseStorage = new BaseStorage();
	
	while (iterator.hasNext()) {
		
		SAMRecord record = iterator.next();
		
		LinkedHashMap<DataType, Object> values = new LinkedHashMap<DataType, Object>();

		Region recordRegion = new Region((long) record.getAlignmentStart(), (long) record.getAlignmentEnd(), request.start.chr);
		
		Feature read = new Feature(recordRegion, values);

		values.put(DataType.ID, record.getReadName());
		
		values.put(DataType.STRAND, BamUtils.getStrand(record, coverageType));
		
		Cigar cigar = new Cigar(read, record.getCigar());
		values.put(DataType.CIGAR, cigar);
		
		String seq = record.getReadString();
		values.put(DataType.SEQUENCE, seq);
		
		
		// Split read into continuous blocks (elements) by using the cigar
		List<ReadPart> parts = Cigar.splitElements(read);
		
		for (ReadPart part : parts) {				 
			
			if (read.values.get(DataType.STRAND) == Strand.FORWARD) {
				forwardBaseStorage.addNucleotideCounts(part);
			} else if (read.values.get(DataType.STRAND) == Strand.REVERSE) {
				reverseBaseStorage.addNucleotideCounts(part);			
			}
		}						
	}
			
	// We are done
	iterator.close();
			
	/* Reads that overlap query regions create nucleotide counts outside the query region.
	 * Remove those extra nucleotide counts, because they don't contain all reads of those regions and would show
	 * wrong information. 
	 */
	Region filterRegion = new Region(start, end, from.chr);
	forwardBaseStorage.filter(filterRegion);
	reverseBaseStorage.filter(filterRegion);

	// Send result		
	LinkedList<Feature> resultList = new LinkedList<Feature>();					
	
	createResultList(from, forwardBaseStorage, resultList, Strand.FORWARD);
	createResultList(from, reverseBaseStorage, resultList, Strand.REVERSE);		
	
	LinkedList<Feature> averageCoverage = CoverageTool.average(resultList, from.chr);
	
	if (request.getRequestedContents().contains(DataType.COVERAGE)) {		
		super.createDataResult(new DataResult(request, resultList));
	}
	
	if (request.getRequestedContents().contains(DataType.COVERAGE_AVERAGE)) {
		super.createDataResult(new DataResult(request, averageCoverage));
	}
}