Java SamReader.getFileHeader方法代码示例

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
private void parseBam(final SamReader reader) {
    Assert.notNull(reader, getMessage(RESOURCE_NOT_FOUND));
    final SAMFileHeader samFileHeader = reader.getFileHeader();

    //check we can read this Bam-file
    //get list of chromosome
    final List<SAMSequenceRecord> list = samFileHeader.getSequenceDictionary().getSequences();
    Assert.notEmpty(list, getMessage(MessagesConstants.WRONG_HEADER_BAM_FILE_EMPTY_FILE));

    //get first chromosome and make a request to the file with this chromosome
    final SAMSequenceRecord samSequenceRecord = list.get(0);
    SAMRecordIterator iterator = reader.query(samSequenceRecord.getSequenceName(),
            Constants.BAM_START_INDEX_TEST, Math.min(Constants.MAX_BAM_END_INDEX_TEST,
                    samSequenceRecord.getSequenceLength()), false);
    Assert.notNull(iterator);
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public static void main(String[] args) throws IOException {
	SamReader reader = SamReaderFactory
			.makeDefault()
			.validationStringency(ValidationStringency.SILENT)
			.open(new File(
					"./data/miniCaviar_IDT_NEB.runA.NA12878.bwa.chrom1.bam"));

	SAMFileHeader samFileHeader = reader.getFileHeader();

	System.out.println(samFileHeader.getTextHeader());

	int i = 0;
	for (SAMRecord record : reader) {
		if (i < 100)
			System.out.println(record.getAlignmentStart());
		i++;
	}
	reader.close();
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * Prepare to iterate through the given SAM records, skipping non-primary alignments
 *
 * @param samReader    must be coordinate sorted
 * @param intervalList Either the list of desired intervals, or null.  Note that if an intervalList is
 *                     passed in that is not coordinate sorted, it will eventually be coordinated sorted by this class.
 * @param useIndex     If true, do indexed lookup to improve performance.  Not relevant if intervalList == null.
 *                     It is no longer the case the useIndex==true can make performance worse.  It should always perform at least
 *                     as well as useIndex==false, and generally will be much faster.
 */
public SamLocusIterator(final SamReader samReader, final IntervalList intervalList, final boolean useIndex) {
    //if (samReader.getFileHeader().getSortOrder() == null || samReader.getFileHeader().getSortOrder() == SAMFileHeader.SortOrder.unsorted) {
    //   LOG.warn("SamLocusIterator constructed with samReader that has SortOrder == unsorted.  ", "" +
    //            "Assuming SAM is coordinate sorted, but exceptions may occur if it is not.");
    //} else if (samReader.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
    //    throw new SAMException("SamLocusIterator cannot operate on a SAM file that is not coordinate sorted.");
    //}
    this.samReader = samReader;
    this.useIndex = useIndex;
    if (intervalList != null) {
        intervals = intervalList.uniqued().getIntervals();
        this.referenceSequenceMask = new IntervalListReferenceSequenceMask(intervalList);
    } else {
        intervals = null;
        this.referenceSequenceMask = new WholeGenomeReferenceSequenceMask(samReader.getFileHeader());
    }
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@Test
public void SamToReadToSamTest() {
  String filePath = "src/test/resources/com/google/cloud/genomics/utils/conversion_test.sam";
  File samInput = new File(filePath);
  SamReader reads = SamReaderFactory.makeDefault().open(samInput);
  SAMFileHeader header = reads.getFileHeader();

  int numReads = 0;
  for (SAMRecord sam : reads){
    Read read = ReadUtils.makeReadGrpc(sam);
    SAMRecord newSam = ReadUtils.makeSAMRecord(read, header );
    final String originalSamString = sam.getSAMString();
    final String postConversionString = newSam.getSAMString();
    assertEquals(originalSamString, postConversionString);
    numReads++;
  }
  assertEquals(19, numReads);//sanity check to make sure we actually read the file
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@Test
public void SamToReadToSamTest() throws IOException {
  String filePath = "src/test/resources/com/google/cloud/genomics/utils/conversion_test.sam";
  File samInput = new File(filePath);
  SamReader reads = SamReaderFactory.makeDefault().open(samInput);
  SAMFileHeader header = reads.getFileHeader();

  int numReads = 0;
  for (SAMRecord sam : reads){
    Read read = ReadUtils.makeRead(sam);
    SAMRecord newSam = ReadUtils.makeSAMRecord(read, header );
    assertEquals(newSam.getSAMString(), sam.getSAMString());
    numReads++;
  }
  assertEquals(19, numReads);//sanity check to make sure we actually read the file
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * Create the BAI index.
 * @param conf the Hadoop configuration
 * @param bamFile the BAM file
 * @param indexFile the BAI file
 * @throws IOException if an error occurs while creating the index
 */
static void createIndex(final Configuration conf, final Path bamFile,
    final Path indexFile) throws IOException {

  final InputStream in = FileSystem.get(conf).open(bamFile);

  final SamReader reader = SamReaderFactory.makeDefault()
      .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS)
      .validationStringency(ValidationStringency.DEFAULT_STRINGENCY)
      .open(SamInputResource.of(in));

  final BAMIndexer indexer =
      new BAMIndexer(indexFile.getFileSystem(conf).create(indexFile),
          reader.getFileHeader());

  for (SAMRecord rec : reader) {
    indexer.processAlignment(rec);

  }

  indexer.finish();
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * Read the SAM header of a SAM file.
 * @param is input stream
 * @return a String with the SAM header
 */
public static String readSAMHeader(final InputStream is) {

  if (is == null) {
    throw new NullPointerException("The input stream is null.");
  }

  // Read SAM file header
  final SamReader reader =
      SamReaderFactory.makeDefault().open(SamInputResource.of(is));
  final SAMFileHeader header = reader.getFileHeader();

  // Close reader
  // reader.close();

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

  return headerTextBuffer.toString();
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public void canFilterFromIntFlag(){
	
	int f_incl= 131;
	int F_excl= 72;
	String chrom= "chrY";
	int from= 1;
	int to= 100;
	
	SamReaderFactory srf=SamReaderFactory.make();
	SamReader samReader= srf.open(new File("test_data/mjb050_oxBS.bam"));
	SAMFileHeader fh= samReader.getFileHeader();
	IntervalList il= new IntervalList(fh);
	Interval interval= new Interval(chrom, from, to);
	il.add(interval);
	
	List<SamRecordFilter> filters= FlagToFilter.flagToFilterList(f_incl, F_excl);

	SamLocusIterator samLocIter= new SamLocusIterator(samReader, il, true);
	samLocIter.setSamFilters(filters);
	while(samLocIter.hasNext()){
		LocusInfo locus= samLocIter.next();
		if(locus.getRecordAndPositions().size() > 0){
			//System.out.println(locus.getPosition() + " " + locus.getRecordAndPositions().size() + " " +
			//		locus.getRecordAndPositions().get(0).getRecord().getFlags());
		}
	}
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@Test
  public void testIteratorWithIntron() {
  	// See https://github.com/samtools/htsjdk/issues/838
  	// Prepare a sam header
String sqHeader = "@HD\tSO:coordinate\tVN:1.0\n"
                      + "@SQ\tSN:chr1\tAS:HG18\tLN:10000000\n";

// Prepare one read with a 500,000 bases skipped
String cigar= "18M500000N18M";
String s1 = "read1\t0\tchr1\t1\t255\t" + cigar + "\t*\t0\t0\tACCTACGTTCAATATTACAGGCGAACATACTTACTA\t*\n";

// Prepare sam input and samReader
String exampleSam = sqHeader + s1 + s1;
ByteArrayInputStream inputStream = new ByteArrayInputStream(exampleSam.getBytes());
SamReader samReader= SamReaderFactory.makeDefault().open(SamInputResource.of(inputStream));

// A small interval to iterate over:
IntervalList il= new IntervalList(samReader.getFileHeader());
il.add(new Interval("chr1", 1, 100));

SamLocusIterator sli= new SamLocusIterator(samReader, il, true);

// Iterate 
long t0= System.currentTimeMillis();
int n= 0;
for (SamLocusIterator.LocusInfo li : sli) {
	n++;
}
long t1= System.currentTimeMillis();
System.err.println("Time to iterate " + n + " loci: " + (t1-t0) / 1000.0 + " sec");
sli.close();
  }
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/** reads local disks or GCS -> header, and PCollection */
private PCollection<GATKRead> ingestReadsAndGrabHeader(final Pipeline pipeline, String filename) throws IOException {

    // input reads
    if (BucketUtils.isCloudStorageUrl(filename)) {
        // set up ingestion on the cloud
        // but read the header locally
        GcsPath path = GcsPath.fromUri(filename);
        InputStream inputstream = Channels.newInputStream(new GcsUtil.GcsUtilFactory().create(pipeline.getOptions())
                .open(path));
        SamReader reader = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT).open(SamInputResource.of(inputstream));
        header = reader.getFileHeader();

        final SAMSequenceDictionary sequenceDictionary = header.getSequenceDictionary();
        final List<SimpleInterval> intervals = IntervalUtils.getAllIntervalsForReference(sequenceDictionary);
        return new ReadsDataflowSource(filename, pipeline).getReadPCollection(intervals, ValidationStringency.SILENT, false);
    } else {
        // ingestion from local file
        try( ReadsDataSource readsSource = new ReadsDataSource(new File(filename)) ) {
            header = readsSource.getHeader();
            List<GATKRead> reads = new ArrayList<>();
            for ( GATKRead read : readsSource ) {
                reads.add(read);
            }
            return pipeline.apply("input ingest",Create.of(reads));
        }
    }
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public static HeaderInfo getHeaderFromBAMFile(Storage.Objects storage, String BAMPath, Iterable<Contig> explicitlyRequestedContigs) throws IOException {
  HeaderInfo result = null;

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

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

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

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

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
private static List<QueryInterval> getAllChunksBalanced(SamReader bam, int countPerContig) {
    List<QueryInterval> ret = new ArrayList<>();
    SAMFileHeader header = bam.getFileHeader();
    for (SAMSequenceRecord s : header.getSequenceDictionary().getSequences()) {
        ret.addAll(getChunksBalanced(bam, s.getSequenceIndex(), countPerContig));
    }
    return ret;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
private static List<QueryInterval> getChunksBalanced(SamReader bam, int sequenceIndex, int retCount) {
    List<QueryInterval> ret = new ArrayList<>();
    BAMIndex index = bam.indexing().getIndex();
    SAMFileHeader header = bam.getFileHeader();
    SAMSequenceRecord s = header.getSequence(sequenceIndex);
    long totalLength = chunksLength(getChunks(index, sequenceIndex, 1, s.getSequenceLength() + 1));
    if (totalLength == 0) {
        return ret;
    }
    int sofar = 0;
    long targetLength = totalLength / retCount;
    int end = s.getSequenceLength();
    int step = s.getSequenceLength() / (100 * retCount);
    if (step < 1) step = 1;
    int start = 1;
    for (int j = step; j < end; j += step) {
        if (j > end) j = end;
        List<Chunk> candidate = getChunks(index, sequenceIndex, start, j);
        long size = chunksLength(candidate);
        if (size < targetLength) {
            // not big enough yet
            continue;
        }
        if (size > targetLength * 2) {
            // too large, search for a good separation point
            // TODO
        }
        // good, emit.
        ret.add(new QueryInterval(sequenceIndex, start, j + 1));
        start = j;
        sofar += size;
        if (ret.size() < retCount) {
            targetLength = (totalLength - sofar) / (retCount - ret.size());
        } else {
            targetLength = totalLength / retCount;
        }

    }
    return ret;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public boolean importIntoHdfs(String weburl, FileSystem fileSystem,
		String path) throws IOException {

	// check if bai is available else download whole file

	final SamReader reader = SamReaderFactory.makeDefault().validationStringency(htsjdk.samtools.ValidationStringency.SILENT).open(
			SamInputResource.of(new URL(weburl)).index(
					new URL(weburl + ".bai")));

	// path in hdfs
	String[] tiles = weburl.split("/");
	String name = tiles[tiles.length - 1];

	String target = HdfsUtil.path(path, name);

	SAMFileHeader header = reader.getFileHeader();
	SAMSequenceDictionary seqDictionary = header.getSequenceDictionary();

	String referenceName = null;

	for (SAMSequenceRecord record : seqDictionary.getSequences()) {

		if (record.getSequenceLength() == REFERENCE_LENGTH) {
			referenceName = record.getSequenceName();
		}
	}

	if (referenceName == null) {
		reader.close();
		error = "No mitochondrial contig found in " + weburl + ".";
		return false;
	}

	FSDataOutputStream out = fileSystem.create(new Path(target));
	SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(
			reader.getFileHeader(), true, out);

	SAMRecordIterator reads = reader.query(referenceName, 0, 0, false);
	int good = 0;
	int bad = 0;
	int written = 0;
	SAMRecord read = null;
	while (reads.hasNext()) {
		try { // hansi style solution TODO!
			read = reads.next();
			good++;
		} catch (Exception e) {
			// e.printStackTrace(s);
			bad++;
		}
		writer.addAlignment(read);
		written++;
	}

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

	System.out.println("Bad reads: " + bad);
	System.out.println("Good reads: " + good);
	System.out.println("Written reads: " + written);
	
	return true;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@Override
public void run() {
	// TODO Auto-generated method stub
	final SamReaderFactory factory =
			SamReaderFactory.makeDefault()
			.enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS, 
					SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
			.validationStringency(ValidationStringency.SILENT);
	final SamReader inputSam = factory.open(new File(mySamFile));
	
	samHeader = inputSam.getFileHeader();
	samHeader.setSortOrder(SortOrder.unsorted);
	SAMRecordIterator iter=inputSam.iterator();
	Set<Entry<String, String>> attr = samHeader.getAttributes();
	
	List<SAMReadGroupRecord> rgs = samHeader.getReadGroups();
	
	SAMReadGroupRecord rg = new SAMReadGroupRecord("cz1");
	rg.setSample("cz1");
	samHeader.addReadGroup(rg);
	
	
	//samHeader.setAttribute("RG", "cz1");
	final SAMFileWriter outSam = new SAMFileWriterFactory().
			makeSAMOrBAMWriter(samHeader,
					true, new File(myOutput));
	
	for(int i=0; i<100; i++) {
		SAMRecord record = iter.next();
		List<SAMTagAndValue> tags = record.getAttributes();
		record.setAttribute("RG", "cz1");
		List<SAMTagAndValue> tags2 = record.getAttributes();
		outSam.addAlignment(record);
	}
	myLogger.info("exit...");
	
	try {
		inputSam.close();
	} catch (IOException e) {
		// TODO Auto-generated catch block
		e.printStackTrace();
	}
	outSam.close();
}