htsjdk.samtools.SAMTextHeaderCodec Java Examples

@Test
public void testAddCommentsToBam() throws Exception {
    final File outputFile = File.createTempFile("addCommentsToBamTest.", BamFileIoUtils.BAM_FILE_EXTENSION);
    outputFile.deleteOnExit();
    runIt(INPUT_FILE, outputFile, commentList);

    final SAMFileHeader newHeader = SamReaderFactory.makeDefault().getFileHeader(outputFile);

    // The original comments are massaged when they're added to the header. Perform the same massaging here,
    // and then compare the lists
    final List<String> massagedComments = new LinkedList<String>();
    for (final String comment : commentList) {
        massagedComments.add(SAMTextHeaderCodec.COMMENT_PREFIX + comment);
    }

    Assert.assertEquals(newHeader.getComments(), massagedComments);
    outputFile.delete();
}
 

private void writeHeader(final SAMFileHeader header) {
	binaryCodec.writeBytes("BAM\001".getBytes(Charset.forName("UTF8")));

	final Writer sw = new StringWriter();
	new SAMTextHeaderCodec().encode(sw, header);

	binaryCodec.writeString(sw.toString(), true, false);

	final SAMSequenceDictionary dict = header.getSequenceDictionary();

	binaryCodec.writeInt(dict.size());
	for (final SAMSequenceRecord rec : dict.getSequences()) {
		binaryCodec.writeString(rec.getSequenceName(), true, true);
		binaryCodec.writeInt   (rec.getSequenceLength());
	}
}
 

/**
 * Update the comments to represent the new field names for the locatable information.
 *
 * @param finalContigColumnName Never {@code null} and never blank.
 * @param finalStartColumnName Never {@code null} and never blank.
 * @param finalEndColumnName Never {@code null} and never blank.
 * @param currentComments Never {@code null}.  Can be empty.
 * @return list of the comments.  Never {@code null} and never empty.
 */
private static List<String> updateStructuredComments(final String finalContigColumnName, final String finalStartColumnName, final String finalEndColumnName, final List<String> currentComments) {

    Utils.nonNull(currentComments);
    Utils.validateArg(!StringUtils.isEmpty(finalContigColumnName), "Contig column name was null or blank, which is invalid.");
    Utils.validateArg(!StringUtils.isEmpty(finalStartColumnName), "Start column name was null or blank, which is invalid.");
    Utils.validateArg(!StringUtils.isEmpty(finalEndColumnName), "End column name was null or blank, which is invalid.");

    // Remove old structured comments, if present.
    final List<String> commentsToWrite = currentComments.stream()
            .filter(c -> !c.startsWith(SAMTextHeaderCodec.COMMENT_PREFIX + CONTIG_COL_COMMENT))
            .filter(c -> !c.startsWith(SAMTextHeaderCodec.COMMENT_PREFIX + START_COL_COMMENT))
            .filter(c -> !c.startsWith(SAMTextHeaderCodec.COMMENT_PREFIX + END_COL_COMMENT)).collect(Collectors.toList());

    // Write out the column headers as a comment
    commentsToWrite.add(CONTIG_COL_COMMENT + finalContigColumnName);
    commentsToWrite.add(START_COL_COMMENT + finalStartColumnName);
    commentsToWrite.add(END_COL_COMMENT + finalEndColumnName);
    return commentsToWrite;
}
 

/**
 * Given an InputStream connected to a fasta dictionary, returns its sequence dictionary
 *
 * Note: does not close the InputStream it's passed
 *
 * @param fastaDictionaryStream InputStream connected to a fasta dictionary
 * @return the SAMSequenceDictionary from the fastaDictionaryStream
 */
public static SAMSequenceDictionary loadFastaDictionary( final InputStream fastaDictionaryStream ) {
    // Don't close the reader when we're done, since we don't want to close the client's InputStream for them
    final BufferedLineReader reader = new BufferedLineReader(fastaDictionaryStream);

    final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();
    final SAMFileHeader header = codec.decode(reader, fastaDictionaryStream.toString());

    // Make sure we have a valid sequence dictionary before continuing:
    if (header.getSequenceDictionary() == null || header.getSequenceDictionary().isEmpty()) {
        throw new UserException.MalformedFile (
                "Could not read sequence dictionary from given fasta stream " +
                        fastaDictionaryStream
        );
    }

    return header.getSequenceDictionary();
}
 

@Override
public void traverse() {
    final SAMFileHeader bamHeader = getHeaderForReads();

    try ( final PrintWriter outputWriter = new PrintWriter(outputFile) ) {
        final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();
        codec.encode(outputWriter, bamHeader);
    } catch (FileNotFoundException e ) {
        throw new UserException.CouldNotCreateOutputFile("Error writing reads header to " + outputFile, e);
    }

    logger.info("Successfully wrote reads header to destination: " + outputFile);
}
 

private static byte[] toByteArray(SAMFileHeader samFileHeader) {
	ExposedByteArrayOutputStream headerBodyOS = new ExposedByteArrayOutputStream();
	OutputStreamWriter outStreamWriter = new OutputStreamWriter(headerBodyOS);
	(new SAMTextHeaderCodec()).encode(outStreamWriter, samFileHeader);

	try {
		outStreamWriter.close();
	} catch (IOException var8) {
		throw new RuntimeException(var8);
	}

	ByteBuffer buf = ByteBuffer.allocate(4);
	buf.order(ByteOrder.LITTLE_ENDIAN);
	buf.putInt(headerBodyOS.size());
	buf.flip();
	byte[] bytes = new byte[buf.limit()];
	buf.get(bytes);
	ByteArrayOutputStream headerOS = new ByteArrayOutputStream();

	try {
		headerOS.write(bytes);
		headerOS.write(headerBodyOS.getBuffer(), 0, headerBodyOS.size());
	} catch (IOException var7) {
		throw new RuntimeException(var7);
	}

	return headerOS.toByteArray();
}
 

/**
 * Writes a header to a BAM file. Might need to regenerate the String
 * version of the header, if one already has both the samFileHeader and the
 * String, use the version of this method which takes both.
 */
static void writeHeader(final BinaryCodec outputBinaryCodec, final SAMFileHeader samFileHeader) {
	// Do not use SAMFileHeader.getTextHeader() as it is not updated when
	// changes to the underlying object are made
	final String headerString;
	final Writer stringWriter = new StringWriter();
	new SAMTextHeaderCodec().encode(stringWriter, samFileHeader, true);
	headerString = stringWriter.toString();

	writeHeader(outputBinaryCodec, samFileHeader, headerString);
}
 

public int streamElPrep(Reducer.Context context, String output, String rg, 
            int threads, SAMRecordIterator SAMit, 
            SAMFileHeader header, String dictFile, boolean updateRG, boolean keepDups, String RGID) throws InterruptedException, IOException, QualityException {
        long startTime = System.currentTimeMillis();
        String customArgs = HalvadeConf.getCustomArgs(context.getConfiguration(), "elprep", "");  
        String[] command = CommandGenerator.elPrep(bin, "/dev/stdin", output, threads, true, rg, null, !keepDups, customArgs);
//        runProcessAndWait(command);
        ProcessBuilderWrapper builder = new ProcessBuilderWrapper(command, null);
        builder.startProcess(true);        
        BufferedWriter localWriter = builder.getSTDINWriter();
        
        // write header
        final StringWriter headerTextBuffer = new StringWriter();
        new SAMTextHeaderCodec().encode(headerTextBuffer, header);
        final String headerText = headerTextBuffer.toString();
        localWriter.write(headerText, 0, headerText.length());
        
        
        SAMRecord sam;
        int reads = 0;
        while(SAMit.hasNext()) {
            sam = SAMit.next();
            if(updateRG)
                sam.setAttribute(SAMTag.RG.name(), RGID);
            String samString = sam.getSAMString();
            localWriter.write(samString, 0, samString.length());
            reads++;
        }
        localWriter.flush();
        localWriter.close();
                
        int error = builder.waitForCompletion();
        if(error != 0)
            throw new ProcessException("elPrep", error);
        long estimatedTime = System.currentTimeMillis() - startTime;
        Logger.DEBUG("estimated time: " + estimatedTime / 1000);
        if(context != null)
            context.getCounter(HalvadeCounters.TIME_ELPREP).increment(estimatedTime);
        return reads;
    }
 

@Test(expectedExceptions = GATKException.ShouldNeverReachHereException.class)
public void testErrorCreateWithInconsistentRegionList() {
    final AnnotatedInterval a1 = new AnnotatedInterval(new SimpleInterval("1", 100, 200),
            ImmutableSortedMap.of("Foo", "5"));
    final AnnotatedInterval a2 = new AnnotatedInterval(new SimpleInterval("1", 100, 200),
            ImmutableSortedMap.of("Foo", "15"));
    final AnnotatedInterval a3 = new AnnotatedInterval(new SimpleInterval("1", 100, 200),
            ImmutableSortedMap.of("Bar", "15"));
    final LineReader reader = BufferedLineReader.fromString("@HD\tVN:1.5\n");
    final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();

    AnnotatedIntervalCollection.create(Arrays.asList(a1,a2,a3), codec.decode(reader, null), Arrays.asList("Foo"));
}
 

/**
 * DEV NOTE: If you are adding attributes that are neither RealMatrix nor a primitive,
 * you must follow the pattern in the constructor (i.e. the Lazy loading pattern).
 * Otherwise, some operations will hang.
 */
private HDF5SVDReadCountPanelOfNormals(final HDF5File file) {
    IOUtils.canReadFile(file.getFile());
    this.file = file;
    sequenceDictionary = new Lazy<>(() -> {
        final String sequenceDictionaryString = file.readStringArray(SEQUENCE_DICTIONARY_PATH)[0];
        return new SAMTextHeaderCodec()
                .decode(BufferedLineReader.fromString(sequenceDictionaryString), file.getFile().getAbsolutePath())
                .getSequenceDictionary();
    });
    originalIntervals = new Lazy<>(() -> HDF5Utils.readIntervals(file, ORIGINAL_INTERVALS_PATH));
    panelIntervals = new Lazy<>(() -> HDF5Utils.readIntervals(file, PANEL_INTERVALS_PATH));
}
 

/**
 * Creates a copy from the SAM File header or an empty list if no sam file header and no comments.
 *
 * Dev note:  This method will be used more when locatable fields are written from the structured comments.
  */
public List<String> getComments() {
    if (getSamFileHeader() == null) {
        return Collections.emptyList();
    } else {
        return getSamFileHeader().getComments().stream().map(c -> c.replaceFirst(SAMTextHeaderCodec.COMMENT_PREFIX, "")).collect(Collectors.toList());
    }
}
 

/**
 * DEV NOTE: If you are adding attributes that are neither RealMatrix nor a primitive,
 * you must follow the pattern in the constructor (i.e. the Lazy loading pattern).
 * Otherwise, some operations will hang.
 */
HDF5SimpleCountCollection(final HDF5File file) {
    Utils.nonNull(file);
    this.file = file;
    sampleName = new Lazy<>(() -> file.readStringArray(SAMPLE_NAME_PATH)[0]);
    sequenceDictionary = new Lazy<>(() -> {
        final String sequenceDictionaryString = file.readStringArray(SEQUENCE_DICTIONARY_PATH)[0];
        return new SAMTextHeaderCodec()
                .decode(BufferedLineReader.fromString(sequenceDictionaryString), file.getFile().getAbsolutePath())
                .getSequenceDictionary();
    });
    intervals = new Lazy<>(() -> HDF5Utils.readIntervals(file, INTERVALS_GROUP_NAME));
    counts = new Lazy<>(() -> new Array2DRowRealMatrix(file.readDoubleMatrix(COUNTS_PATH)));
}
 

/**
 * Private helper method for {@link #convertHeaderlessHadoopBamShardToBam} that takes a SAMFileHeader and writes it
 * to the provided `OutputStream`, correctly encoded for the BAM format and preceded by the BAM magic bytes.
 *
 * @param samFileHeader SAM header to write
 * @param outputStream stream to write the SAM header to
 */
private static void writeBAMHeaderToStream( final SAMFileHeader samFileHeader, final OutputStream outputStream ) {
    final BlockCompressedOutputStream blockCompressedOutputStream = new BlockCompressedOutputStream(outputStream, (File)null);
    final BinaryCodec outputBinaryCodec = new BinaryCodec(new DataOutputStream(blockCompressedOutputStream));

    final String headerString;
    final Writer stringWriter = new StringWriter();
    new SAMTextHeaderCodec().encode(stringWriter, samFileHeader, true);
    headerString = stringWriter.toString();

    outputBinaryCodec.writeBytes(ReadUtils.BAM_MAGIC);

    // calculate and write the length of the SAM file header text and the header text
    outputBinaryCodec.writeString(headerString, true, false);

    // write the sequences binarily.  This is redundant with the text header
    outputBinaryCodec.writeInt(samFileHeader.getSequenceDictionary().size());
    for (final SAMSequenceRecord sequenceRecord: samFileHeader.getSequenceDictionary().getSequences()) {
        outputBinaryCodec.writeString(sequenceRecord.getSequenceName(), true, true);
        outputBinaryCodec.writeInt(sequenceRecord.getSequenceLength());
    }

    try {
        blockCompressedOutputStream.flush();
    } catch (final IOException ioe) {
        throw new RuntimeIOException(ioe);
    }
}
 

@VisibleForTesting
static SAMFileHeader createSamFileHeaderWithCommentsOnly(final List<String> comments) {
    final LineReader reader = BufferedLineReader.fromString("");
    final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();
    final SAMFileHeader result = codec.decode(reader, null);

    // We grab all of the comments from the first decode, just in case the codec has some default comments.
    final List<String> finalComments = new ArrayList<>();
    finalComments.addAll(result.getComments());
    finalComments.addAll(comments);
    result.setComments(finalComments);

    return result;
}
 

/**
 * @return copy of the sam file header created from the given list of strings.  {@code null} is not possible
 */
@VisibleForTesting
static SAMFileHeader createSamFileHeader(final List<String> samFileHeaderAsStrings) {

    final LineReader reader = BufferedLineReader.fromString(StringUtils.join(samFileHeaderAsStrings, "\n"));
    final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();
    return codec.decode(reader, null);
}
 

@Override
public SampleLocatableMetadata readActualHeader(final LineIterator reader) {
    final List<String> samHeaderLines = new ArrayList<>(SAM_HEADER_LINES_INITIAL_CAPACITY);
    //we check that the SAM header lines and the column header line are present in the correct order, then return the mandatory column header
    boolean isSAMHeaderPresent = false;
    while (reader.hasNext()) {
        final String line = reader.peek();
        if (line.startsWith(CopyNumberFormatsUtils.COMMENT_PREFIX)) {
            isSAMHeaderPresent = true;
            samHeaderLines.add(line);
            reader.next();
        } else {
            if (!isSAMHeaderPresent) {
                throw new UserException.MalformedFile("SAM header lines must be at the beginning of the file.");
            } else if (!line.startsWith(COLUMN_HEADER_STRING)) {
                throw new UserException.MalformedFile("File does not have a column header.");
            } else {
                //we just peeked at the column header line, so we need to advance past it
                reader.next();
                break;
            }
        }
    }
    final SAMFileHeader samFileHeader = new SAMTextHeaderCodec()
            .decode(BufferedLineReader.fromString(StringUtils.join(samHeaderLines, System.lineSeparator())), null);
    return MetadataUtils.fromHeader(samFileHeader, Metadata.Type.SAMPLE_LOCATABLE);
}
 

/**
 * Converts a SAM header to text form
 * @param header the header
 * @return the text
 */
public static String getHeaderAsString(SAMFileHeader header) {
  final StringWriter writer = new StringWriter();
  new SAMTextHeaderCodec().encode(writer, header);
  return writer.toString();
}
 

private void writeSequenceDictionary(final String path, final SAMSequenceDictionary sequenceDictionary) {
    final StringWriter stringWriter = new StringWriter();
    new SAMTextHeaderCodec().encode(stringWriter, new SAMFileHeader(sequenceDictionary));
    file.makeStringArray(path, stringWriter.toString());
}
 

private SAMFileHeader getHeader() {
    try (final LineReader lineReader = new BufferedLineReader(BucketUtils.openFile(path.toString()))) { //use BucketUtils.openFile to handle *.gz files
        return new SAMTextHeaderCodec().decode(lineReader, getSource());
    }
}
 

private void writeSequenceDictionary(final SAMSequenceDictionary sequenceDictionary) {
    final StringWriter stringWriter = new StringWriter();
    new SAMTextHeaderCodec().encode(stringWriter, new SAMFileHeader(sequenceDictionary));
    file.makeStringArray(SEQUENCE_DICTIONARY_PATH, stringWriter.toString());
}
 

public void writeHeader(final SAMFileHeader header) {
  final StringWriter headerTextBuffer = new StringWriter();
  new SAMTextHeaderCodec().encode(headerTextBuffer, header);
  final String headerText = headerTextBuffer.toString();
  super.writeHeader(headerText);
}
 

static SAMFileHeader readHeader(final BinaryCodec stream, final ValidationStringency validationStringency,
		final String source) throws IOException {

	final byte[] buffer = new byte[4];
	stream.readBytes(buffer);
	if (!Arrays.equals(buffer, "BAM\1".getBytes())) {
		throw new IOException("Invalid BAM file header");
	}

	final int headerTextLength = stream.readInt();
	final String textHeader = stream.readString(headerTextLength);
	final SAMTextHeaderCodec headerCodec = new SAMTextHeaderCodec();
	headerCodec.setValidationStringency(validationStringency);
	final SAMFileHeader samFileHeader = headerCodec.decode(new StringLineReader(textHeader), source);

	final int sequenceCount = stream.readInt();
	if (samFileHeader.getSequenceDictionary().size() > 0) {
		// It is allowed to have binary sequences but no text sequences, so
		// only validate if both are present
		if (sequenceCount != samFileHeader.getSequenceDictionary().size()) {
			throw new SAMFormatException("Number of sequences in text header ("
					+ samFileHeader.getSequenceDictionary().size() + ") != number of sequences in binary header ("
					+ sequenceCount + ") for file " + source);
		}
		for (int i = 0; i < sequenceCount; i++) {
			final SAMSequenceRecord binarySequenceRecord = readSequenceRecord(stream, source);
			final SAMSequenceRecord sequenceRecord = samFileHeader.getSequence(i);
			if (!sequenceRecord.getSequenceName().equals(binarySequenceRecord.getSequenceName())) {
				throw new SAMFormatException("For sequence " + i
						+ ", text and binary have different names in file " + source);
			}
			if (sequenceRecord.getSequenceLength() != binarySequenceRecord.getSequenceLength()) {
				throw new SAMFormatException("For sequence " + i
						+ ", text and binary have different lengths in file " + source);
			}
		}
	} else {
		// If only binary sequences are present, copy them into
		// samFileHeader
		final List<SAMSequenceRecord> sequences = new ArrayList<SAMSequenceRecord>(sequenceCount);
		for (int i = 0; i < sequenceCount; i++) {
			sequences.add(readSequenceRecord(stream, source));
		}
		samFileHeader.setSequenceDictionary(new SAMSequenceDictionary(sequences));
	}

	return samFileHeader;
}
 

/** Writes out a HaplotypeMap file with the contents of this map. */
public void writeToFile(final File file) {
    try {
        final BufferedWriter out = new BufferedWriter(new OutputStreamWriter(IOUtil.openFileForWriting(file)));
        final FormatUtil format = new FormatUtil();

        // Write out the header
        if (this.header != null) {
            final SAMTextHeaderCodec codec = new SAMTextHeaderCodec();
            codec.encode(out, this.header);
        }

        // Write the header for the entries.
        out.write("#CHROMOSOME\tPOSITION\tNAME\tMAJOR_ALLELE\tMINOR_ALLELE\tMAF\tANCHOR_SNP\tPANELS");
        out.newLine();

        final List<HaplotypeMapFileEntry> entries = new ArrayList<>();
        for (final HaplotypeBlock block : this.getHaplotypes()) {
            String anchor = null;
            final SortedSet<Snp> snps = new TreeSet<>(block.getSnps());

            for (final Snp snp : snps) {
                entries.add(new HaplotypeMapFileEntry(snp.getChrom(), snp.getPos(), snp.getName(),
                        snp.getAllele1(), snp.getAllele2(), snp.getMaf(), anchor, snp.getFingerprintPanels()));

                if (anchor == null) {
                    anchor = snp.getName();
                }
            }
        }

        Collections.sort(entries);
        for (final HaplotypeMapFileEntry entry : entries) {
            out.write(entry.chromosome + "\t");
            out.write(format.format(entry.position) + "\t");
            out.write(entry.snpName + "\t");
            out.write((char)entry.majorAllele + "\t");
            out.write((char)entry.minorAllele + "\t");
            out.write(format.format(entry.minorAlleleFrequency) + "\t");
            if (entry.anchorSnp != null) {
                out.write(entry.anchorSnp);
            }
            out.write("\t");
            if (entry.getPanels() != null) {
                out.write(entry.getPanels());
            }
            out.newLine();
        }
        out.flush();
        out.close();
    }
    catch (IOException ioe) {
        throw new PicardException("Error writing out haplotype map to file: " + file.getAbsolutePath(), ioe);
    }
}
 

/**
 * This is factored out of doWork only for unit testing.
 */
int writeSamText(PrintStream printStream) {
    try {
        final CloseableIterator<SAMRecord> samRecordsIterator;
        final SamReader samReader = SamReaderFactory.makeDefault()
                .referenceSequence(REFERENCE_SEQUENCE)
                .open(SamInputResource.of(INPUT));

        // if we are only using the header or we aren't using intervals, then use the reader as the iterator.
        // otherwise use the SamRecordIntervalIteratorFactory to make an interval-ed iterator
        if (HEADER_ONLY || INTERVAL_LIST == null) {
            samRecordsIterator = samReader.iterator();
        } else {
            IOUtil.assertFileIsReadable(INTERVAL_LIST);

            final List<Interval> intervals = IntervalList.fromFile(INTERVAL_LIST).uniqued().getIntervals();
            samRecordsIterator = new SamRecordIntervalIteratorFactory().makeSamRecordIntervalIterator(samReader, intervals, samReader.hasIndex());
        }
        final AsciiWriter writer = new AsciiWriter(printStream);
        final SAMFileHeader header = samReader.getFileHeader();
        if (!RECORDS_ONLY) {
            if (header.getTextHeader() != null) {
                writer.write(header.getTextHeader());
            } else {
                // Headers that are too large are not retained as text, so need to regenerate text
                new SAMTextHeaderCodec().encode(writer, header, true);
            }
        }
        if (!HEADER_ONLY) {
            while (samRecordsIterator.hasNext()) {
                final SAMRecord rec = samRecordsIterator.next();

                if (printStream.checkError()) {
                    return 1;
                }

                if (this.ALIGNMENT_STATUS == AlignmentStatus.Aligned && rec.getReadUnmappedFlag()) continue;
                if (this.ALIGNMENT_STATUS == AlignmentStatus.Unaligned && !rec.getReadUnmappedFlag()) continue;

                if (this.PF_STATUS == PfStatus.PF && rec.getReadFailsVendorQualityCheckFlag()) continue;
                if (this.PF_STATUS == PfStatus.NonPF && !rec.getReadFailsVendorQualityCheckFlag()) continue;
                writer.write(rec.getSAMString());
            }
        }
        writer.flush();
        if (printStream.checkError()) {
            return 1;
        }
        CloserUtil.close(writer);
        CloserUtil.close(samRecordsIterator);
        return 0;
    } catch (IOException e) {
        throw new PicardException("Exception writing SAM text", e);
    }
}