htsjdk.samtools.SAMFileWriter Java Examples

The following examples show how to use htsjdk.samtools.SAMFileWriter. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example #1
Source File: SamUtilsTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 6 votes vote down vote up
public static void convertSamToBam(File outFile, File indexFile, File samFile) throws IOException {
  try (RecordIterator<SAMRecord> sam = new SkipInvalidRecordsIterator(samFile.getPath(), new SamClosedFileReader(samFile, null, null, SamUtils.getSingleHeader(samFile)))) {
    final SAMFileHeader header = sam.header().clone();
    header.setSortOrder(SAMFileHeader.SortOrder.coordinate);
    SamUtils.addProgramRecord(header);
    SamUtils.updateRunId(header);

    final SAMFileWriterFactory fact = new SAMFileWriterFactory();
    try (SAMFileWriter writer = fact.makeBAMWriter(header, true, outFile)) {
      while (sam.hasNext()) {
        final SAMRecord record = sam.next();
        try {
          writer.addAlignment(record);
        } catch (final IllegalArgumentException iae) {
          throw new NoTalkbackSlimException(iae.getMessage().replaceAll(Pattern.quote("SAMFileWriterImpl.addAlignment for "), ""));
        }
      }
    }
  }
  try {
    BamIndexer.saveBamIndex(outFile, indexFile);
  } catch (final UnindexableDataException e) {
    Diagnostic.warning("Cannot create BAM index: " + e.getMessage());
  }
}
 
Example #2
Source File: SamBamUtils.java    From chipster with MIT License 6 votes vote down vote up
public static void sortSamBam(File samBamFile, File sortedBamFile) {
	
	SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
	SAMFileReader reader = new SAMFileReader(IOUtil.openFileForReading(samBamFile));
	SAMFileWriter writer = null;
	try {
		
		reader.getFileHeader().setSortOrder(SAMFileHeader.SortOrder.coordinate);
		writer = new SAMFileWriterFactory().makeBAMWriter(reader.getFileHeader(), false, sortedBamFile);
		Iterator<SAMRecord> iterator = reader.iterator();
		while (iterator.hasNext()) {
			writer.addAlignment(iterator.next());
		}
		
	} finally {
		closeIfPossible(reader);
		closeIfPossible(writer);
	}
}
 
Example #3
Source File: FilterBamByTag.java    From Drop-seq with MIT License 6 votes vote down vote up
/**
 * Just work through the reads one at a time.
 * @param in
 * @param out
 */
void processUnpairedMode (final SamReader in, final SAMFileWriter out, final Set<String> values, final File summaryFile, Integer mapQuality ,boolean allowPartial) {
	FilteredReadsMetric m = new FilteredReadsMetric();
	ProgressLogger progLog = new ProgressLogger(log);
	for (final SAMRecord r : in) {
		progLog.record(r);
		boolean filterFlag = filterRead(r, this.TAG, values, this.ACCEPT_TAG, mapQuality, allowPartial);
		if (!filterFlag) { 
			out.addAlignment(r);
			m.READS_ACCEPTED++;
		}
		else {
			m.READS_REJECTED++;
		}
	}
	writeSummary(summaryFile, m);
	CloserUtil.close(in);
	out.close();
	reportAndCheckFilterResults(m);
}
 
Example #4
Source File: PreprocessingTools.java    From halvade with GNU General Public License v3.0 6 votes vote down vote up
public int callElPrep(String input, String output, String rg, int threads, 
        SAMRecordIterator SAMit,
        SAMFileHeader header, String dictFile, boolean updateRG, boolean keepDups, String RGID) throws InterruptedException, QualityException {
    
    SAMRecord sam;
    SAMFileWriterFactory factory = new SAMFileWriterFactory();
    SAMFileWriter Swriter = factory.makeSAMWriter(header, true, new File(input));
    
    int reads = 0;
    while(SAMit.hasNext()) {
        sam = SAMit.next();
        if(updateRG)
            sam.setAttribute(SAMTag.RG.name(), RGID);
        Swriter.addAlignment(sam);
        reads++;
    }
    Swriter.close();
    
    String customArgs = HalvadeConf.getCustomArgs(context.getConfiguration(), "elprep", "");  
    String[] command = CommandGenerator.elPrep(bin, input, output, threads, true, rg, null, !keepDups, customArgs);
    long estimatedTime = runProcessAndWait("elPrep", command);
    if(context != null)
        context.getCounter(HalvadeCounters.TIME_ELPREP).increment(estimatedTime);
    
    return reads;
}
 
Example #5
Source File: CollectGcBiasMetricsTest.java    From picard with MIT License 6 votes vote down vote up
public File build (final List<SAMRecordSetBuilder> setBuilder, final File unsortedSam, final SAMFileHeader header) throws IOException {
    final File sortedSam = VcfTestUtils.createTemporaryIndexedFile("CollectGcBias", ".bam");

    final SAMFileWriter writer = new SAMFileWriterFactory()
            .setCreateIndex(true).makeBAMWriter(header, false, unsortedSam);

    for (final SAMRecordSetBuilder subSetBuilder : setBuilder) {
        for (final SAMRecord record : subSetBuilder) {
            writer.addAlignment(record);
        }
    }
    writer.close();

    final SortSam sorter = new SortSam();
    final String[] args = new String[] {
            "INPUT=" + unsortedSam.getAbsolutePath(),
            "OUTPUT=" + sortedSam.getAbsolutePath(),
            "SORT_ORDER=coordinate"
    };

    sorter.instanceMain(args);

    return sortedSam;
}
 
Example #6
Source File: CollapseBarcodesInPlace.java    From Drop-seq with MIT License 6 votes vote down vote up
public void processOnlyPrimary () {
       final SamHeaderAndIterator inputs = openInputs();
	CloseableIterator<SAMRecord> inputSam = inputs.iterator;
	SAMFileHeader header = inputs.header;
	header.addComment("Edit distance collapsed tag " +  this.PRIMARY_BARCODE + " to new tag " + this.OUT_BARCODE+ " with edit distance "+ this.EDIT_DISTANCE);
       SAMFileWriter writer= new SAMFileWriterFactory().makeSAMOrBAMWriter(header, true, this.OUTPUT);

	// gather up the barcodes that exist in the BAM
       final SamHeaderAndIterator inputs2 = openInputs();
	ObjectCounter<String> barcodes = new BamTagHistogram().getBamTagCounts(inputs2.iterator, this.PRIMARY_BARCODE,this.MINIMUM_MAPPING_QUALITY, this.FILTER_PCR_DUPLICATES);
       CloserUtil.close(inputs2.iterator);

	// filter barcodes by #reds in each barcode.
	barcodes=filterBarcodesByNumReads(barcodes, this.MIN_NUM_READS_NONCORE);

	// collapse them
	Map<String, String> childParentBarcodes=collapseBarcodes(this.MIN_NUM_READS_CORE, this.NUM_CORE_BARCODES, barcodes, this.FIND_INDELS, this.EDIT_DISTANCE);
	// iterate through the reads and retag with the proper reads.
	// log.info("STUFF");
	retagReads(inputSam, writer, childParentBarcodes, this.PRIMARY_BARCODE, this.OUT_BARCODE);
	// collapsed.size();

	CloserUtil.close(inputSam);
	writer.close();
}
 
Example #7
Source File: SamBamUtils.java    From chipster with MIT License 6 votes vote down vote up
public String printSamBam(InputStream samBamStream, int maxRecords) throws IOException {
	SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
	SAMFileReader in = new SAMFileReader(samBamStream);
	SAMFileHeader header = in.getFileHeader();
	ByteArrayOutputStream buffer = new ByteArrayOutputStream();
	SAMFileWriter out = new SAMFileWriterFactory().makeSAMWriter(header, true, buffer);
	int i = 0;
	try {
		for (final SAMRecord rec : in) {
			if (i > maxRecords) {
				break;
			}
			out.addAlignment(rec);
			i++;
		}
	} finally {
		closeIfPossible(out);
	}

	if (i > maxRecords) {
		buffer.write("SAM/BAM too long for viewing, truncated here!\n".getBytes());
	}
	
	return buffer.toString();
}
 
Example #8
Source File: CollapseTagWithContext.java    From Drop-seq with MIT License 6 votes vote down vote up
/**
 * If the number of records exceeds the number of records allowed in memory, spill to disk.
 * @param groupingIter
 * @param writer
 * @param outMetrics
 * @param header
 */
private void lowMemoryIteration (PeekableGroupingIterator<SAMRecord> groupingIter,									 
								 SAMFileWriter writer, PrintStream outMetrics, SAMFileHeader header) {
	log.info("Running (slower) memory efficient mode");				
       while (groupingIter.hasNext()) {
       	// for this group, get a SortingCollection.  Note that this is not used for sorting.  It is merely
		// an unsorted collection if there might be more objects than can fit in RAM.
       	SortingCollection<SAMRecord> sortingCollection = SortingCollection.newInstance(SAMRecord.class, new BAMRecordCodec(header), NO_OP_COMPARATOR, this.MAX_RECORDS_IN_RAM);

       	// you have to grab the next element, in case it's the first of the group but not the first group!
       	sortingCollection.add(groupingIter.next()); 
       	
       	// spool the reads for a whole group into the sorting collection to operate on - the code uses a multi-pass approach so we can't just iterate over the grouping iterator.
       	while (groupingIter.hasNextInGroup())         		
       		sortingCollection.add(groupingIter.next());
       	
       	// wrap up the sorting collection for adding records.
       	sortingCollection.doneAdding();
       	sortingCollection.setDestructiveIteration(false);
       	
       	processContext(sortingCollection, writer, false, outMetrics);        	
       }	
}
 
Example #9
Source File: CollectDuplicateMetricsTester.java    From picard with MIT License 6 votes vote down vote up
@Override
public File getOutput() {
    final File output;
    try {
        output = File.createTempFile("CollectDuplicateMetrics", ".sam");
    } catch (IOException e) {
        throw new PicardException("problems creating output file", e);
    }
    output.deleteOnExit();

    try(SAMFileWriter writer = new SAMFileWriterFactory().makeWriter(samRecordSetBuilder.getHeader(), true, output.toPath(), this.fastaFiles.get(samRecordSetBuilder.getHeader()))) {
        samRecordSetBuilder.getRecords().forEach(a -> {
            final SAMRecord b = a.deepCopy();
            b.setDuplicateReadFlag(duplicateFlags.get(samRecordToDuplicatesFlagsKey(a)));
            writer.addAlignment(b);
        });
    }

    return output;
}
 
Example #10
Source File: CollapseTagWithContext.java    From Drop-seq with MIT License 6 votes vote down vote up
/**
 * With this method, we keep all the records in memory 
 * @param groupingIter
 * @param mapQualityPredicate
 * @param requiredTagPredicate
 * @param writer
 * @param outMetrics
 */
private void fasterIteration (PeekableGroupingIterator<SAMRecord> groupingIter,	SAMFileWriter writer, PrintStream outMetrics) {
	log.info("Running fast single iteration mode");
	int maxNumInformativeReadsInMemory=1000; // the starting value is just for reporting purposes.
       while (groupingIter.hasNext()) {
       	
       	List<SAMRecord> informativeRecs = new ArrayList<>();
       	// you have to grab the next element, in case it's the first of the group but not the first group!
       	informativeRecs.add(groupingIter.next()); 
       	
       	while (groupingIter.hasNextInGroup())         		
       		informativeRecs.add(groupingIter.next());        	        	        	
       	
       	// you have all the informative reads.
       	// do some additional logging if the number of reads is bigger than what you've seen before.
       	boolean verbose = false;
       	if (informativeRecs.size()>maxNumInformativeReadsInMemory) {
       		maxNumInformativeReadsInMemory=informativeRecs.size();
       		log.info("Max informative reads in memory [" + maxNumInformativeReadsInMemory +"]");
       		verbose=true;
       	}
       	
       	// get context.
       	processContext(informativeRecs, writer, verbose, outMetrics);    	
       }		
}
 
Example #11
Source File: CollapseTagWithContext.java    From Drop-seq with MIT License 6 votes vote down vote up
private void processContext (Iterable<SAMRecord> i, SAMFileWriter writer, boolean verbose, PrintStream outMetrics) {
	PeekableIterator<SAMRecord> iter = new PeekableIterator<>(i.iterator());
   	if (!iter.hasNext()) return;

	// get context.
	String context = getContextString(iter.peek(), this.CONTEXT_TAGS);

	// get barcode counts.
	ObjectCounter<String> barcodeCounts = getBarcodeCounts (iter, this.COLLAPSE_TAG, this.COUNT_TAGS, this.COUNT_TAGS_EDIT_DISTANCE);
	if (this.MIN_COUNT > 1 & !this.MUTATIONAL_COLLAPSE) barcodeCounts.filterByMinCount(this.MIN_COUNT);
	
	Map<String, String> collapseMap = collapseBarcodes(barcodeCounts, this.FIND_INDELS, this.EDIT_DISTANCE, this.ADAPTIVE_ED_MIN, this.ADAPTIVE_ED_MAX, this.MIN_COUNT, this.MUTATIONAL_COLLAPSE_PATH_ED, verbose, outMetrics, context, this.ADAPTIVE_ED_METRICS_ED_LIST);
	iter = new PeekableIterator<>(i.iterator());
	retagBarcodedReads(iter, barcodeCounts, collapseMap, this.DROP_SMALL_COUNTS, writer, this.COLLAPSE_TAG, this.OUT_TAG);        	

}
 
Example #12
Source File: CollapseTagWithContext.java    From Drop-seq with MIT License 6 votes vote down vote up
private void retagBarcodedReads (Iterator<SAMRecord> informativeRecs, ObjectCounter<String> barcodeCounts, Map<String, String> collapseMap, boolean dropSmallCounts, SAMFileWriter writer,
		String collapseTag, String outTag) {
	
	Set<String> expectedBarcodes = null;
	// already validated that if dropSmallCounts is true, then the minNumObservations > 1.
	if (dropSmallCounts)
		// use all the remaining barcodes that have counts.
		expectedBarcodes = new HashSet<>(barcodeCounts.getKeys());
	while (informativeRecs.hasNext()) {
		SAMRecord r = informativeRecs.next();
		String tagValue = r.getStringAttribute(collapseTag);
		// if the tag was not set, then don't set it.
		if (tagValue!=null) {
			// tag was set,
			// if the tagValue is not in the expected barocde list and the barcode list is populated, then don't add this read and short circuit to next read.
			if (expectedBarcodes!=null && !expectedBarcodes.contains(tagValue))
				continue;
			// tag was set, set it.
			// is it in the map?  If so, update the tag value.
			if (collapseMap.containsKey(tagValue))
				tagValue = collapseMap.get(tagValue);
			r.setAttribute(outTag, tagValue);
		}
		writer.addAlignment(r);
	}		
}
 
Example #13
Source File: DetectBeadSynthesisErrors.java    From Drop-seq with MIT License 6 votes vote down vote up
/**
 * For each problematic cell, replace cell barcodes positions with N.
 * Take the replaced bases and prepend them to the UMI, and trim the last <X> bases off the end of the UMI.
 */
private void cleanBAM (final Map<String, BeadSynthesisErrorData> errorBarcodesWithPositions, final Map<String, String> intendedSequenceMap) {
	log.info("Cleaning BAM");
       final SamHeaderAndIterator headerAndIterator = SamFileMergeUtil.mergeInputs(INPUT, true);
	SamHeaderUtil.addPgRecord(headerAndIterator.header, this);

	SAMFileWriter writer= new SAMFileWriterFactory().setCreateIndex(CREATE_INDEX).makeSAMOrBAMWriter(headerAndIterator.header, true, OUTPUT);
	ProgressLogger pl = new ProgressLogger(log);
	for (SAMRecord r: new IterableAdapter<>(headerAndIterator.iterator)) {
		pl.record(r);
		r=padCellBarcodeFix(r, errorBarcodesWithPositions, intendedSequenceMap, this.CELL_BARCODE_TAG, this.MOLECULAR_BARCODE_TAG, this.EXTREME_BASE_RATIO);
		if (r!=null)
			writer.addAlignment(r);
	}
	CloserUtil.close(headerAndIterator.iterator);
	writer.close();
}
 
Example #14
Source File: DetectBeadSubstitutionErrors.java    From Drop-seq with MIT License 6 votes vote down vote up
private void repairBAM (final BottomUpCollapseResult result) {

		final SamHeaderAndIterator headerAndIterator = SamFileMergeUtil.mergeInputs(this.INPUT, true);
		SamHeaderUtil.addPgRecord(headerAndIterator.header, this);
		headerAndIterator.header.addComment("Bottom-up edit distance collapse tag " + this.CELL_BARCODE_TAG +" with edit distance " + this.EDIT_DISTANCE+ " filtering ambiguous neighbors=" + this.FILTER_AMBIGUOUS);
		SAMFileWriter writer= new SAMFileWriterFactory().setCreateIndex(CREATE_INDEX).makeSAMOrBAMWriter(headerAndIterator.header, true, OUTPUT);

		ProgressLogger pl = new ProgressLogger(log);
		log.info("Repairing BAM");
		for (SAMRecord r: new IterableAdapter<>(headerAndIterator.iterator)) {
			pl.record(r);
			r=repairBarcode(r, result);
			if (r!=null)
				writer.addAlignment(r);
		}
		CloserUtil.close(headerAndIterator.iterator);
		writer.close();
		log.info("Repair Complete");
	}
 
Example #15
Source File: BamSlicerApplication.java    From hmftools with GNU General Public License v3.0 6 votes vote down vote up
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();
}
 
Example #16
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
private void addAlignmentsForBestFragmentMapqStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final String sequence, final int[] mapqs) {
    boolean reverse = false;
    int alignmentStart = 1;
    for (final int mapq : mapqs) {
        final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
        alignedRecord.setReadName(unmappedRecord.getReadName());
        alignedRecord.setReadBases(unmappedRecord.getReadBases());
        alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
        alignedRecord.setReferenceName(sequence);
        alignedRecord.setAlignmentStart(alignmentStart);
        alignmentStart += 10; // Any old position will do
        alignedRecord.setReadNegativeStrandFlag(reverse);
        reverse = !reverse;
        alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
        alignedRecord.setMappingQuality(mapq);
        alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
        alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
        alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
        alignedRecord.setMateUnmappedFlag(true);
        writer.addAlignment(alignedRecord);
    }
}
 
Example #17
Source File: ReplaceSamHeader.java    From picard with MIT License 6 votes vote down vote up
private void standardReheader(final SAMFileHeader replacementHeader) {
    final SamReader recordReader = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).validationStringency(ValidationStringency.SILENT).open(INPUT);
    if (replacementHeader.getSortOrder() != recordReader.getFileHeader().getSortOrder()) {
        throw new PicardException("Sort orders of INPUT (" + recordReader.getFileHeader().getSortOrder().name() +
                ") and HEADER (" + replacementHeader.getSortOrder().name() + ") do not agree.");
    }
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(replacementHeader, true, OUTPUT);

    final ProgressLogger progress = new ProgressLogger(Log.getInstance(ReplaceSamHeader.class));
    for (final SAMRecord rec : recordReader) {
        rec.setHeader(replacementHeader);
        writer.addAlignment(rec);
        progress.record(rec);
    }
    writer.close();
    CloserUtil.close(recordReader);
}
 
Example #18
Source File: SortSam.java    From picard with MIT License 6 votes vote down vote up
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    final SamReader reader = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);
    ;
    reader.getFileHeader().setSortOrder(SORT_ORDER.getSortOrder());
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(reader.getFileHeader(), false, OUTPUT);
    writer.setProgressLogger(
            new ProgressLogger(log, (int) 1e7, "Wrote", "records from a sorting collection"));

    final ProgressLogger progress = new ProgressLogger(log, (int) 1e7, "Read");
    for (final SAMRecord rec : reader) {
        writer.addAlignment(rec);
        progress.record(rec);
    }

    log.info("Finished reading inputs, merging and writing to output now.");

    CloserUtil.close(reader);
    writer.close();
    return 0;
}
 
Example #19
Source File: RevertSam.java    From picard with MIT License 6 votes vote down vote up
RevertSamWriter(
        final boolean outputByReadGroup,
        final Map<String, SAMFileHeader> headerMap,
        final Map<String, File> outputMap,
        final SAMFileHeader singleOutHeader,
        final File singleOutput,
        final boolean presorted,
        final SAMFileWriterFactory factory,
        final File referenceFasta) {

    this.outputByReadGroup = outputByReadGroup;
    if (outputByReadGroup) {
        singleWriter = null;
        for (final Map.Entry<String, File> outputMapEntry : outputMap.entrySet()) {
            final String readGroupId = outputMapEntry.getKey();
            final File output = outputMapEntry.getValue();
            final SAMFileHeader header = headerMap.get(readGroupId);
            final SAMFileWriter writer = factory.makeWriter(header, presorted, output, referenceFasta);
            writerMap.put(readGroupId, writer);
        }
    } else {
        singleWriter = factory.makeWriter(singleOutHeader, presorted, singleOutput, referenceFasta);
    }
}
 
Example #20
Source File: SamFormatConverter.java    From picard with MIT License 6 votes vote down vote up
/**
 * Convert a file from one of sam/bam/cram format to another based on the extension of output.
 *
 * @param input             input file in one of sam/bam/cram format
 * @param output            output to write converted file to, the conversion is based on the extension of this filename
 * @param referenceSequence the reference sequence to use, necessary when reading/writing cram
 * @param createIndex       whether or not an index should be written alongside the output file
 */
public static void convert(final File input, final File output, final File referenceSequence, final Boolean createIndex) {
    IOUtil.assertFileIsReadable(input);
    IOUtil.assertFileIsWritable(output);
    final SamReader reader = SamReaderFactory.makeDefault().referenceSequence(referenceSequence).open(input);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeWriter(reader.getFileHeader(), true, output, referenceSequence);
    if (createIndex && writer.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
        throw new PicardException("Can't CREATE_INDEX unless sort order is coordinate");
    }

    final ProgressLogger progress = new ProgressLogger(Log.getInstance(SamFormatConverter.class));
    for (final SAMRecord rec : reader) {
        writer.addAlignment(rec);
        progress.record(rec);
    }
    CloserUtil.close(reader);
    writer.close();
}
 
Example #21
Source File: GatherBamFiles.java    From picard with MIT License 6 votes vote down vote up
/**
 * Simple implementation of a gather operations that uses SAMFileReaders and Writers in order to concatenate
 * multiple BAM files.
 */
private static void gatherNormally(final List<File> inputs, final File output, final boolean createIndex, final boolean createMd5,
                                   final File referenceFasta) {
    final SAMFileHeader header;
    {
        header = SamReaderFactory.makeDefault().referenceSequence(referenceFasta).getFileHeader(inputs.get(0));
    }

    final SAMFileWriter out = new SAMFileWriterFactory().setCreateIndex(createIndex).setCreateMd5File(createMd5).makeSAMOrBAMWriter(header, true, output);

    for (final File f : inputs) {
        log.info("Gathering " + f.getAbsolutePath());
        final SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceFasta).open(f);
        for (final SAMRecord rec : in) out.addAlignment(rec);
        CloserUtil.close(in);
    }

    out.close();
}
 
Example #22
Source File: SmartSamWriter.java    From rtg-tools with BSD 2-Clause "Simplified" License 6 votes vote down vote up
/**
 * Constructor
 * @param out where records should eventually be sent
 */
public SmartSamWriter(SAMFileWriter out) {
  super(BUFFER_SIZE, new SAMRecordCoordinateComparator() {
    @Override
    public int compare(final SAMRecord samRecord1, final SAMRecord samRecord2) {
      final int cmp = super.compare(samRecord1, samRecord2);
      if (cmp != 0) {
        return cmp;
      }
      final SAMReadGroupRecord rg1 = samRecord1.getReadGroup();
      final SAMReadGroupRecord rg2 = samRecord2.getReadGroup();
      if (rg1 == rg2) {
        return 0;
      }
      return rg1 == null ? -1 : rg2 == null ? 1 : rg1.getReadGroupId().compareTo(rg2.getReadGroupId());
    }
  });
  mOut = out;
}
 
Example #23
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
private void addAlignmentForMostStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final MostDistantStrategyAlignmentSpec spec,
        final boolean reverse) {
    final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
    alignedRecord.setReadName(unmappedRecord.getReadName());
    alignedRecord.setReadBases(unmappedRecord.getReadBases());
    alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
    alignedRecord.setReferenceName(spec.sequence);
    alignedRecord.setAlignmentStart(spec.alignmentStart);
    alignedRecord.setReadNegativeStrandFlag(reverse);
    alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
    alignedRecord.setMappingQuality(spec.mapQ);
    alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
    alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
    alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
    alignedRecord.setMateUnmappedFlag(true);
    writer.addAlignment(alignedRecord);
}
 
Example #24
Source File: CleanSam.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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;
}
 
Example #25
Source File: BAMTestUtil.java    From Hadoop-BAM with MIT License 5 votes vote down vote up
public static File writeBamFileWithLargeHeader() throws IOException {
  SAMRecordSetBuilder samRecordSetBuilder =
      new SAMRecordSetBuilder(true, SAMFileHeader.SortOrder.queryname);
  for (int i = 0; i < 1000; i++) {
    int chr = 20;
    int start1 = (i + 1) * 1000;
    int start2 = start1 + 100;
    samRecordSetBuilder.addPair(String.format("test-read-%03d", i), chr, start1,
        start2);
  }

  final File bamFile = File.createTempFile("test", ".bam");
  bamFile.deleteOnExit();
  SAMFileHeader samHeader = samRecordSetBuilder.getHeader();
  StringBuffer sb = new StringBuffer();
  for (int i = 0; i < 1000000; i++) {
    sb.append("0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789");
  }
  samHeader.addComment(sb.toString());
  final SAMFileWriter bamWriter = new SAMFileWriterFactory()
      .makeSAMOrBAMWriter(samHeader, true, bamFile);
  for (final SAMRecord rec : samRecordSetBuilder.getRecords()) {
    bamWriter.addAlignment(rec);
  }
  bamWriter.close();

  return bamFile;
}
 
Example #26
Source File: CollectHsMetricsTest.java    From picard with MIT License 5 votes vote down vote up
/** Writes the contents of a SAMRecordSetBuilder out to a file. */
private File writeBam(final SAMRecordSetBuilder builder, final File f) {
    try (final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(builder.getHeader(), false, f)) {
        builder.forEach(out::addAlignment);
    }
    return f;
}
 
Example #27
Source File: RevertSam.java    From picard with MIT License 5 votes vote down vote up
void addAlignment(final SAMRecord rec) {
    final SAMFileWriter writer;
    if (outputByReadGroup) {
        writer = writerMap.get(rec.getReadGroup().getId());
    } else {
        writer = singleWriter;
    }
    writer.addAlignment(rec);
}
 
Example #28
Source File: SamBamUtils.java    From chipster with MIT License 5 votes vote down vote up
public void normaliseBam(File bamFile, File normalisedBamFile) {

		// Read in a BAM file and its header
		SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
		SAMFileReader reader = new SAMFileReader(IOUtil.openFileForReading(bamFile));
		SAMFileWriter writer = null;
		try {
			SAMFileHeader normalisedHeader = reader.getFileHeader();

			// Alter the chromosome names in header's SAMSequenceDictionary
			SAMSequenceDictionary normalisedDictionary = new SAMSequenceDictionary();
			for (SAMSequenceRecord sequenceRecord : normalisedHeader.getSequenceDictionary().getSequences()) {

				// Normalise chromosome
				String sequenceName = chromosomeNormaliser.normaliseChromosome(sequenceRecord.getSequenceName());
				normalisedDictionary.addSequence(new SAMSequenceRecord(sequenceName, sequenceRecord.getSequenceLength()));
			}
			normalisedHeader.setSequenceDictionary(normalisedDictionary);

			// Write new BAM file with normalised chromosome names
			writer = new SAMFileWriterFactory().makeBAMWriter(normalisedHeader, true, normalisedBamFile);
			for (final SAMRecord rec : reader) {
				rec.setHeader(normalisedHeader);
				writer.addAlignment(rec);
			}
			
		} finally {
			closeIfPossible(reader);
			closeIfPossible(writer);
		}
	}
 
Example #29
Source File: SamBamUtils.java    From chipster with MIT License 5 votes vote down vote up
private static void closeIfPossible(SAMFileWriter writer) {
	if (writer != null) {
		try {
			writer.close();
		} catch (Exception e) {
			// Ignore
		}
	}
}
 
Example #30
Source File: Cram2Bam.java    From cramtools with Apache License 2.0 5 votes vote down vote up
private static SAMFileWriter createSAMFileWriter(Params params, CramHeader cramHeader,
		SAMFileWriterFactory samFileWriterFactory) throws IOException {
	/*
	 * building sam writer, sometimes we have to go deeper to get to the
	 * required functionality:
	 */
	SAMFileWriter writer = null;
	if (params.outputFastq) {
		if (params.cramURL == null) {
			writer = new FastqSAMFileWriter(System.out, null, cramHeader.getSamFileHeader());
		} else {
			writer = new FastqSAMFileWriter(Utils.getFileName(params.cramURL), false, cramHeader.getSamFileHeader());

		}
	} else if (params.outputFastqGz) {
		if (params.cramURL == null) {
			GZIPOutputStream gos = new GZIPOutputStream(System.out);
			PrintStream ps = new PrintStream(gos);
			writer = new FastqSAMFileWriter(ps, null, cramHeader.getSamFileHeader());
		} else {
			writer = new FastqSAMFileWriter(Utils.getFileName(params.cramURL), true, cramHeader.getSamFileHeader());

		}
	} else if (params.outputFile == null) {
		OutputStream os = new BufferedOutputStream(System.out);
		if (params.outputBAM) {
			writer = new SAMFileWriterFactory().makeBAMWriter(cramHeader.getSamFileHeader(), true, os);
		} else {
			writer = Utils.createSAMTextWriter(samFileWriterFactory, os, cramHeader.getSamFileHeader(),
					params.printSAMHeader);
		}
	} else {
		writer = samFileWriterFactory.makeSAMOrBAMWriter(cramHeader.getSamFileHeader(), true, params.outputFile);
	}
	return writer;
}