Java Code Examples for htsjdk.samtools.util.IOUtil#openFileForWriting()

@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(this.ANNOTATIONS_FILE);
	IOUtil.assertFileIsWritable(this.OUTPUT);
	
	OverlapDetector<Gene> od = GeneAnnotationReader.loadAnnotationsFile(this.ANNOTATIONS_FILE, this.SEQUENCE_DICTIONARY);
	
	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(OUTPUT));
	
	for (Gene g: od.getAll()) {
		for (Gene.Transcript t: g) {
			RefFlatRecord r = convertTranscript(t);
			out.println(r.toString());
		}
	}
	out.close();
	return 0;
}
 

@Override
protected int doWork() {

	IOUtil.assertFileIsReadable(INPUT);
	IOUtil.assertFileIsWritable(OUTPUT);
	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(OUTPUT));

	writeHeader(out);

	TagOfTagResults<String,String> results= getResults(this.INPUT, this.PRIMARY_TAG, this.SECONDARY_TAG, this.FILTER_PCR_DUPLICATES, this.MINIMUM_MAPPING_QUALITY);

	for (String k: results.getKeys()) {
		Set<String> values = results.getValues(k);
		writeStats(k, values, out);
	}
	out.close();

	return(0);
}
 

private void writeContigReport (final File outFile, final ObjectCounter<ContigResult> contigResults) {
	PrintStream writer = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(outFile));
	List<String> header = new ArrayList<>();
	header.add("INPUT_1="+this.INPUT_1.toString());
	header.add("INPUT_2="+this.INPUT_2.toString());
	header.add("READ_QUALITY="+this.READ_QUALITY);
	header.add("TRIM_CONTIG_STRING="+this.TRIM_CONTIG_STRING);
	String h = StringUtils.join(header, "\t");
	writer.print("#");
	writer.println(h);

	String [] colNames = {"ORIGINAL_CONTIG", "NEW_CONTIGS", "NUM_READS", "MAPPED_UNIQUELY"};
	writer.println(StringUtil.join("\t", colNames));
	List<ContigResult> crList = contigResults.getKeysOrderedByCount(true);
	for (ContigResult cr: crList) {
		String [] body = {cr.getOldContig(), StringUtil.join(",", cr.getNewContigs()), Integer.toString(contigResults.getCountForKey(cr)), Boolean.toString(cr.isNewReadMapsUniquely())};
		writer.println(StringUtil.join("\t", body));
	}
	writer.close();
}
 

private void writeCommonSubstitutionReport (final BottomUpCollapseResult result, final BarcodeSubstitutionCollection commonChanges, final File outFile) {
	if (outFile==null) return;
	PrintStream outReport = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(outFile));
	String [] header = {"position", "intended_base", "substitution_base", "count", "repaired"};
	outReport.println(StringUtil.join("\t", header));

	String [] bases= {"A", "C", "G", "T"};

	BarcodeSubstitutionCollection allChanges = result.gatherAllPatterns();
	List<Integer> positions = allChanges.getPositions();
	for (int position : positions)
		for (String intendedBase: bases)
			for (String neighborBase: bases)
				if (!intendedBase.equals(neighborBase)) {
					int count = allChanges.getCount(intendedBase, neighborBase, position);
					int commonCount = commonChanges.getCount(intendedBase, neighborBase, position);
					boolean used = commonCount >0;
					String [] line = {Integer.toString(position+1), intendedBase, neighborBase, Integer.toString(count), Boolean.toString(used).toUpperCase()};
					outReport.println(StringUtil.join("\t", line));
				}
	outReport.close();
}
 

void writeReport (final ObjectCounter<String> umisPerCell, final Collection<BarcodeNeighborGroup> neighborGroups, final Map<String, String> intendedSequenceMap, final Map<String, Double> umiBias, final File outFile) {
	if (outFile==null) return;  // if the output file is null, this is a no-op.
	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(outFile));

	IntendedSequenceBuilder b = new IntendedSequenceBuilder(umisPerCell, umiBias);

	// write header.
	String [] header = {"intended_sequence", "related_sequences", "num_related", "deleted_base", "deleted_base_pos", "non_incorporated_rate", "intended_UMIs", "related_median_UMIs", "intended_TBias", "related_median_TBias"};
	out.println(StringUtils.join(header, "\t"));

	// for each neighbor group, write out a report line.
	for (BarcodeNeighborGroup ng: neighborGroups) {
		// this is slightly janky, but the intended sequence map holds raw neighbor barcodes to the intended sequence, so you only need the first barcode in the list.
		String intendedSequence = intendedSequenceMap.get(ng.getNeighborCellBarcodes().iterator().next());
		IntendedSequence is = b.build(intendedSequence, ng);
		writeLine(is, out);
	}
	CloserUtil.close(out);
}
 

@Test
public void testReadBamList() throws IOException {
    // Create a bam_list with one relative path and one absolute, and confirm that reader
    // resolves them correctly.
    final File bamList = TestUtils.getTempReportFile("testReadBamList.", ".bam_list");
    final PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(bamList));
    final String relative = "foo.bam";
    out.println(relative);
    final String absolute = "/an/absolute/path.bam";
    out.println(absolute);
    out.close();
    final List<File> expected = Arrays.asList(
      new File(bamList.getCanonicalFile().getParent(), relative),
      new File(absolute)
    );
    Assert.assertEquals(expected, FileListParsingUtils.readFileList(bamList));
    // Confirm that when reading a symlink to a bam_list, relative paths are resolved relative to the directory
    // of the actually bam_list file, not the directory containing the symlink.
    final File otherDir = Files.createTempDirectory("testReadBamList").toFile();
    otherDir.deleteOnExit();
    final File symlink = new File(otherDir, "testReadBamList.bam_list");
    symlink.deleteOnExit();
    Files.createSymbolicLink(symlink.toPath(), bamList.toPath());
    Assert.assertEquals(expected, FileListParsingUtils.readFileList(symlink));
}
 

private void writeOutput (ObjectCounter<String> data, ObjectCounter<String> aggregateCounts, Map<String, List<String>> mapping, MapBarcodesByEditDistance m) {				
	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(new File ("/downloads/arpy/mutational_collapse_testdata.result_jim.txt")));			
	String [] header = {"sequence",  "counts", "parent",  "edist",  "fam_seqs", "fam_counts"};
	out.println(StringUtils.join(header, "\t"));
	
	for (String parentSeq: mapping.keySet()) {		
		if (parentSeq.equals("ACCTAGATTCAGTTAGCCAG"))
				System.out.println();
		List<String> sequences = mapping.get(parentSeq);														
		int famSeqs=sequences.size()+1;
		String [] line = {parentSeq, Integer.toString(data.getCountForKey(parentSeq)), parentSeq, "0", Integer.toString(famSeqs) ,Integer.toString(aggregateCounts.getCountForKey(parentSeq))};
		out.println(StringUtils.join(line, "\t"));
														
		for (String v: sequences) {
			int ed = HammingDistance.getHammingDistance(parentSeq, v);
			// for merged results, the family seqs size is always 1 and the fam counts is always 0.
			String [] line2 = {v, Integer.toString(data.getCountForKey(v)), parentSeq, Integer.toString(ed), "1", "0"};
			out.println(StringUtils.join(line2, "\t"));																				
		}			
	}
			
	out.close();
}
 

@Override
protected int doWork() {
	// validation
	IOUtil.assertFileIsReadable(INPUT);
	IOUtil.assertFileIsWritable(OUTPUT);
	IOUtil.assertFileIsReadable(this.CELL_BC_FILE);
	
	List<String> cellBarcodes = ParseBarcodeFile.readCellBarcodeFile(this.CELL_BC_FILE);
	log.info("Using " + cellBarcodes.size() + " cells in analysis");
	IntervalList snpIntervals = IntervalList.fromFile(INTERVALS);
	log.info("Using " + snpIntervals.getIntervals().size() + " SNP intervals in analysis");
	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(OUTPUT));
	writeHeader(out);

	SamReader reader = SamReaderFactory.makeDefault().enable(SamReaderFactory.Option.EAGERLY_DECODE).open(this.INPUT);			
	
	SNPUMIBasePileupIterator sbpi = getIter(reader, snpIntervals, cellBarcodes);
	
	MultiCellDigitalAlleleCountsIterator multiIter = new MultiCellDigitalAlleleCountsIterator(new DigitalAlleleCountsIterator(sbpi, BASE_QUALITY));

	// sort cell barcodes alphabetically for output.
	Collections.sort(cellBarcodes);
	@SuppressWarnings("unused")
	int counter=0;
	while (multiIter.hasNext()) {
		MultiCellDigitalAlleleCounts mcdac = multiIter.next();
		processMCDAC(cellBarcodes, mcdac, out, AUTO_FLUSH_OUTPUTS);
		counter++;
		if (counter%PROGRESS_RATE==0) log.info("Processed " + counter + " SNPs");
	}
	log.info("Processed " + counter +" total SNPs");
	out.close();
	multiIter.close();
	return 0;
}
 

@Override
 protected int doWork() {
     IOUtil.assertFileIsReadable(SEQUENCE_DICTIONARY);
     IOUtil.assertFileIsReadable(GTF);
     IOUtil.assertFileIsWritable(this.OUTPUT);
     initialize();

     FilteringGTFParser parser = parseGTF();

     try {
         PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(OUTPUT));
         writeHeader(out);

         // if no enhancement needed, just write out the results.
         if (!ENHANCE_GTF)
	writeRecords (out, parser);
else {

             EnhanceGTFRecords e = new EnhanceGTFRecords();
             List<GTFRecord> records = e.enhanceGTFRecords(parser);
             Collections.sort(records, new GenomicOrderComparator(dict));
             writeRecords(out, records);
         }
         out.close();
         return 0;
     } finally {
         CloserUtil.close(parser);
     }
 }
 

@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(INPUT);
	IOUtil.assertFileIsWritable(OUTPUT);

	ObjectCounter<String> counter = getBamTagCounts(INPUT, this.TAG, this.MINIMUM_MAPPING_QUALITY, this.FILTER_PCR_DUPLICATES);

       PrintStream writer = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.OUTPUT));
       writeHeader(writer);
       writeHistogram(counter, writer);

	return 0;
}
 

public void writeDenseDgeFile(final File output, final boolean formatAsInteger, final DgeHeader dgeHeader) {
       try {
		IOUtil.assertFileIsWritable(output);
		List<String> cellBarcodes = this.getCellBarcodes();
		List<String> genes = this.getGenes();
		BufferedWriter out = new BufferedWriter(new OutputStreamWriter(IOUtil.openFileForWriting(output)));

           if (dgeHeader != null)
			new DgeHeaderCodec().encode(out, dgeHeader);

		// write header
		List<String> header = new ArrayList<>(cellBarcodes.size()+1);
		header.add("GENE");
		for (String c: cellBarcodes)
               header.add(c);
		String h = StringUtils.join(header, "\t");
		out.write(h);
		out.newLine();

		// write body
		for (String gene: genes) {
               List<String> line = new ArrayList<>(cellBarcodes.size()+1);
               line.add(gene);
               double [] expressionByGene = this.getExpression(gene);
               for (double exp: expressionByGene)
				line.add(formatExpressionValue(exp, formatAsInteger));
               String b = StringUtils.join(line, "\t");
               out.write(b); out.newLine();
           }
		out.close();
	} catch (IOException e) {
		throw new RuntimeException("Trouble writing " + output.getAbsolutePath(), e);
	}
}
 

private void writeReport (final BottomUpCollapseResult result, final BottomUpCollapseResult resultClean, final UMIsPerCellResult umiResult) {
	PrintStream outReport = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.OUTPUT_REPORT));

	// write comments section, each line starts with a "#"
	outReport.println("# FILTER_AMBIGUOUS="+FILTER_AMBIGUOUS);
	outReport.println("# MIN_UMIS_PER_CELL="+MIN_UMIS_PER_CELL);
	outReport.println("# UMI_BIAS_THRESHOLD="+MIN_UMIS_PER_CELL);
	outReport.println("# EDIT_DISTANCE="+MIN_UMIS_PER_CELL);
	outReport.println("#");
	outReport.println("# TOTAL_BARCODES_TESTED="+umiResult.getNumCellBarocodesTested());
	outReport.println("# BARCODES_COLLAPSED="+result.getUnambiguousSmallBarcodes().size());
	outReport.println("# ESTIMATED_UMIS_COLLAPSED="+getTotalAmbiguousUMIs(result.getUnambiguousSmallBarcodes(), umiResult.getUmisPerCell()));
	outReport.println("# AMBIGUOUS_BARCODES="+result.getAmbiguousBarcodes().size());
	outReport.println("# ESTIMATED_AMBIGUOUS_UMIS="+getTotalAmbiguousUMIs(result.getAmbiguousBarcodes(), umiResult.getUmisPerCell()));
	outReport.println("# POLY_T_BIASED_BARCODES="+umiResult.getPolyTBiasedBarcodes());
	outReport.println("# POLY_T_BIASED_BARRCODE_UMIS="+umiResult.getPolyTBiasedUMIs());
	outReport.println("# POLY_T_POSITION="+umiResult.getPolyTPosition());

	/// write header
	String [] header= {"intended_barcode", "neighbor_barcode", "intended_size", "neighbor_size", "position", "intended_base", "neighbor_base", "repaired"};
	outReport.println(StringUtil.join("\t", header));

	ObjectCounter<String> umiCounts=umiResult.getUmisPerCell();
	Iterator<String> smalls = result.getUnambiguousSmallBarcodes().iterator();
	while (smalls.hasNext()) {
		String small=smalls.next();
		String large = result.getLargerRelatedBarcode(small);
		BarcodeSubstitutionPair p = new BarcodeSubstitutionPair(large, small);
		String cleanLarger = resultClean.getLargerRelatedBarcode(small);
		boolean repaired = cleanLarger!=null;

		String [] body = {large, small, Integer.toString(umiCounts.getCountForKey(large)), Integer.toString(umiCounts.getCountForKey(small)),
				Integer.toString(p.getPosition()+1), p.getIntendedBase(), p.getNeighborBase(), Boolean.toString(repaired).toUpperCase()};
		outReport.println(StringUtil.join("\t", body));
	}

	CloserUtil.close(outReport);
}
 

@Override
protected int doWork() {
	// primer detection if requested.
	if (this.PRIMER_SEQUENCE!=null)
		this.detectPrimerTool = new DetectPrimerInUMI(this.PRIMER_SEQUENCE);

	PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(OUTPUT_STATS));

	UMIIterator iterator = prepareUMIIterator();
	BiasedBarcodeCollection biasedBarcodeCollection = findBiasedBarcodes(iterator, out, this.SUMMARY, UMI_BIAS_BASE);
	Map<String, BeadSynthesisErrorData> errorBarcodesWithPositions = biasedBarcodeCollection.getBiasedBarcodes();

	ObjectCounter<String> umisPerCell = biasedBarcodeCollection.getUMICounts();
	Map<String, Double> umiBias = biasedBarcodeCollection.getUMIBias();

	// Group T-biased barcodes into groups of neighbors.  The key is the padded sequence the neighbors share
	Map<String, BarcodeNeighborGroup> barcodeNeighborGroups = buildBarcodeNeighborGroups(errorBarcodesWithPositions.values(), this.EXTREME_BASE_RATIO);

	// find intended sequences for these neighbors
	Map<String, String> intendedSequenceMap = findIntendedBarcodes(barcodeNeighborGroups, umisPerCell.getKeys(), umiBias);

	// write report about which sequences were collapsed and why
	writeReport(umisPerCell, barcodeNeighborGroups.values(), intendedSequenceMap, umiBias, this.REPORT);


	// clean up the BAM if desired.
	if (this.OUTPUT!=null)
		cleanBAM(errorBarcodesWithPositions, intendedSequenceMap);
	return 0;
}
 

@Test
public void testFileListParsing() throws IOException {
    // In the temp dir, create 2 files matching the pattern "testFileListParsing.pattern.*.txt"
    final String PATTERN = "testFileListParsing.pattern.*.txt";
    final File patternFile1 = TestUtils.getTempReportFile("testFileListParsing.pattern.", ".txt");
    final File patternFile2 = TestUtils.getTempReportFile("testFileListParsing.pattern.", ".txt");
    final File noPatternFile1 = TestUtils.getTempReportFile("testFileListParsing.no.pattern.", ".txt");
    final File noPatternFile2 = TestUtils.getTempReportFile("testFileListParsing.no.pattern.", ".txt");

    final File tempDir = patternFile1.getParentFile();

    final File fileList = TestUtils.getTempReportFile("testFileListParsing.", ".file_list");
    final PrintStream out = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(fileList));
    out.println(noPatternFile2.getAbsolutePath());
    out.close();

    List<File> expandedList = FileListParsingUtils.expandFileList(Arrays.asList(new File(tempDir, PATTERN), noPatternFile1, fileList));
    final List<File> expected = Arrays.asList(
            patternFile1,
            patternFile2,
            noPatternFile1,
            noPatternFile2
    );

    Collections.sort(expected);
    Collections.sort(expandedList);
    Assert.assertEquals(expected, expandedList);
}
 

@Override
protected int doWork() {
	ObjectCounter<ContigResult> contigResults = null;
	// key = original gene name, value = collection of mappings.
	// Map<String, GeneResult> geneResults = null;
	INPUT_1=FileListParsingUtils.expandFileList(INPUT_1);
	INPUT_2=FileListParsingUtils.expandFileList(INPUT_2);
	INPUT_1.stream().forEach(x-> IOUtil.assertFileIsReadable(x));
	INPUT_2.stream().forEach(x-> IOUtil.assertFileIsReadable(x));
	// if (this.CONTIG_REPORT==null && this.GENE_REPORT==null) log.error ("What's the point of running me if there's no output defined?");
	if (this.CONTIG_REPORT!=null) {
		IOUtil.assertFileIsWritable(this.CONTIG_REPORT);
		contigResults = new ObjectCounter<>();
	}
	
	PrintStream readWriter=null;
	if (this.READ_REPORT!=null) {
		IOUtil.assertFileIsWritable(this.READ_REPORT);
		readWriter = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.READ_REPORT));	
		writeReadHeader(this.TAG_NAMES, readWriter);
	}
	
	/*
	if (this.GENE_REPORT!=null) {
		IOUtil.assertFileIsWritable(this.GENE_REPORT);
		if (this.GENE_EXON_TAG==null) {
			log.error("If GENE_REPORT is specified, a GENE_EXON_TAG must also be specified");
			return 1;
		}

		geneResults = new HashMap<>();
	}
	*/
	
	log.info("Reading INPUT_1");
	PeekableIterator<List<SAMRecord>> oIter = getReadIterator(this.INPUT_1, this.READ_QUALITY, true);
	log.info("Reading INPUT_2");
	PeekableIterator<List<SAMRecord>> nIter = getReadIterator(this.INPUT_2, null, false);

	QueryNameJointIterator qnji = new QueryNameJointIterator(oIter, nIter);

	int counter =0;
	while (qnji.hasNext()) {
		if (counter%1000000==0) log.info("Number of reads processed [" + counter + "]");
		JointResult jr = qnji.next();
		counter++;

		List<SAMRecord> r1 = jr.getOne();
		List<SAMRecord> r2 = jr.getTwo();

		if (contigResults!=null) {
			ContigResult cr = evaluateByContig(r1, r2);
			contigResults.increment(cr);
		}
		if (readWriter!=null) {
			writeReadReportLine(r1, r2, this.TAG_NAMES, readWriter);
		}
		/*
		if (geneResults!=null)
			geneResults=evaluateByGene(r1, r2, geneResults, this.GENE_EXON_TAG);
		*/

	}
	log.info("Number of reads processed [" + counter + "]");

	if (this.CONTIG_REPORT!=null) writeContigReport(this.CONTIG_REPORT, contigResults);
	if (readWriter!=null) readWriter.close();
	
	// if (this.GENE_REPORT!=null) writeGeneReport (this.GENE_REPORT, geneResults);
	return 0;
}
 

@Override
protected int doWork() {
	int vc = validateCommands();
	if (vc>0) return vc;

	if (this.COUNT_TAGS_EDIT_DISTANCE>0) this.medUMI = new MapBarcodesByEditDistance(false);

	med = new MapBarcodesByEditDistance(false, this.NUM_THREADS, 0);
	
	PrintStream outMetrics = null;
	if (this.ADAPTIVE_ED_METRICS_FILE!=null) {
		outMetrics = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.ADAPTIVE_ED_METRICS_FILE));
		writeAdaptiveEditDistanceMetricsHeader(this.ADAPTIVE_ED_METRICS_ED_LIST, outMetrics);
	}
	
	if (this.MUTATIONAL_COLLAPSE_METRICS_FILE!=null) {
		med = new MapBarcodesByEditDistance(true, this.NUM_THREADS, 1000);
		outMetrics = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.MUTATIONAL_COLLAPSE_METRICS_FILE));
		writeMutationalCollapseMetricsHeader(this.ADAPTIVE_ED_METRICS_ED_LIST, outMetrics);
	}

	IOUtil.assertFileIsReadable(INPUT);
       IOUtil.assertFileIsWritable(OUTPUT);

       SamReader reader = SamReaderFactory.makeDefault().open(INPUT);
       SAMFileHeader header =  reader.getFileHeader();
       SortOrder sortOrder= header.getSortOrder();
       
       SAMFileWriter writer = getWriter (reader);
       final ObjectSink<SAMRecord> recSink = new SamWriterSink(writer);
               
       PeekableGroupingIterator<SAMRecord> groupingIter = orderReadsByTagsPeekable(reader, this.COLLAPSE_TAG, this.CONTEXT_TAGS, this.READ_MQ, this.OUT_TAG, recSink);

       log.info("Collapsing tag and writing results");

       if (!LOW_MEMORY_MODE) 
       	fasterIteration(groupingIter, writer, outMetrics);
       else
       	lowMemoryIteration(groupingIter, writer, outMetrics, header);
               
       log.info("Re-sorting output BAM in "+ sortOrder.toString()+ " if neccesary");
       CloserUtil.close(groupingIter);
       CloserUtil.close(reader);
       writer.close();
       if (outMetrics!=null) CloserUtil.close(outMetrics);
       log.info("DONE");
	return 0;
}