Java SAMFileHeader类代码示例

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
private boolean checkHeader(SAMFileHeader header){
List<SAMSequenceRecord> sequences = header.getSequenceDictionary().getSequences();
HashSet<String> map = new HashSet<String>();

//load kourami panel sequence names
BufferedReader br;
try{
    br = new BufferedReader(new InputStreamReader(new GZIPInputStream(new FileInputStream(HLA.MSAFILELOC + File.separator + "All_FINAL_with_Decoy.fa.gz"))));
    String curline = "";
    while((curline = br.readLine())!=null){
	if(curline.charAt(0) == ('>'))
	    map.add(curline.substring(1));
    }
    br.close();
}catch(IOException ioe){
    ioe.printStackTrace();
}

//check if input bam has sequences to kourami panel
for(SAMSequenceRecord ssr : sequences){
    if(!map.contains(ssr.getSequenceName()))
	return false;
}
return true;
   }
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
public static SAMFileHeader transfer(SamHeaderInfo samHeaderInfo){
    SAMFileHeader header = new SAMFileHeader();
    // read groups
    List<ReadGroupInfo> readGroupInfoList = CollectionConverter.asJavaList(samHeaderInfo.getReadGroupInfos());
    List<SAMReadGroupRecord> samReadGroupRecords = readGroupInfoList.stream()
            .map(SAMReadGroupRecordTransfer::transfer)
            .collect(Collectors.toList());
    header.setReadGroups(samReadGroupRecords);
    // Sequence records
    // TODO 获取id的方法不通用 fix me
    int contigCount = samHeaderInfo.getRefContigInfo().getContigIds().size();
    for(int id = 0; id < contigCount; id ++) {
        header.addSequence(new SAMSequenceRecord(
                samHeaderInfo.getRefContigInfo().getName(id),
                samHeaderInfo.getRefContigInfo().getLength(id)));
    }
    if(samHeaderInfo.sorted()) {
        header.setAttribute("SO", "coordinate");
    }
    return header;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
/**
 * setup read groups for the specified read groups and sample names
 *
 * @param header       the header to set
 * @param readGroupIDs the read group ID tags
 * @param sampleNames  the sample names
 * @return the adjusted SAMFileHeader
 */
public static SAMFileHeader createEnumeratedReadGroups(SAMFileHeader header, List<String> readGroupIDs, List<String> sampleNames) {
    if (readGroupIDs.size() != sampleNames.size()) {
        throw new ReviewedGATKException("read group count and sample name count must be the same");
    }

    List<SAMReadGroupRecord> readGroups = new ArrayList<SAMReadGroupRecord>();

    int x = 0;
    for (; x < readGroupIDs.size(); x++) {
        SAMReadGroupRecord rec = new SAMReadGroupRecord(readGroupIDs.get(x));
        rec.setSample(sampleNames.get(x));
        readGroups.add(rec);
    }
    header.setReadGroups(readGroups);
    return header;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
/**
 * Calculates the distance between two genomeLocs across contigs (if necessary).
 * <p>
 * Returns minDistance(other) if in same contig.
 * Works with intervals!
 * Uses the SAMFileHeader to extract the size of the contigs and follows the order in the dictionary.
 *
 * @param other         the genome loc to compare to
 * @param samFileHeader the contig information
 * @return the sum of all the bases in between the genomeLocs, including entire contigs
 */
public long distanceAcrossContigs(GenomeLoc other, SAMFileHeader samFileHeader) {
    if (onSameContig(other))
        return minDistance(other);

    // add the distance from the first genomeLoc to the end of it's contig and the distance from the
    // second genomeLoc to the beginning of it's contig.
    long distance = 0;
    if (contigIndex < other.contigIndex) {
        distance += samFileHeader.getSequence(contigIndex).getSequenceLength() - stop;
        distance += other.start;
    } else {
        distance += samFileHeader.getSequence(other.contigIndex).getSequenceLength() - other.stop;
        distance += start;
    }

    // add any contig (in its entirety) in between the two genomeLocs
    for (int i = Math.min(this.contigIndex, other.contigIndex) + 1; i < Math.max(this.contigIndex, other.contigIndex); i++) {
        distance += samFileHeader.getSequence(i).getSequenceLength();
    }
    return distance;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
public static StaticData getStaticDataInstance(RefContigInfo refContigInfo,
                                               boolean useGVCF,
                                               SamHeaderInfo samHeaderInfo,
                                               VcfHeaderInfo vcfHeaderInfo) {
    StaticData data = new StaticData();
    SAMSequenceDictionary samSequenceDictionary = SAMSequenceDictTransfer.transfer(refContigInfo);
    data.genomeLocParser = new GenomeLocParser(samSequenceDictionary);
    samHeaderInfo.addReadGroupInfo(ReadGroupInfo.apply("rg1", "sample1"));
    SAMFileHeader header = SAMHeaderTransfer.transfer(samHeaderInfo);
    List<SAMReadGroupRecord> readGroupInfos = header.getReadGroups();
    List<String> samples = new ArrayList<>();
    for (SAMReadGroupRecord readGroup : readGroupInfos) {
        samples.add(readGroup.getSample());
    }
    data.haplotypeCaller = new HaplotypeCaller(data.genomeLocParser, samples, useGVCF);
    data.basic2SAMRecordTransfer = new Basic2SAMRecordTransfer(header);

    VCFCodec codec = new VCFCodec();
    VCFHeaderLineIterable headerLineIterable = new VCFHeaderLineIterable(vcfHeaderInfo);
    data.vcfFileHeader = (VCFHeader) codec.readActualHeader(headerLineIterable);
    data.codec = codec;
    return data;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
protected static java.util.List<SAMRecord> getReads(String filePath, SamHeaderInfo headerInfo) {

        String realPath = TestRealignerTargetCreator.class.getResource(filePath).getFile();

        SAMFileHeader header = SAMHeaderTransfer.transfer(headerInfo);
        SAMLineParser parser = new SAMLineParser(header);
        java.util.List<SAMRecord> result = new java.util.ArrayList<>();
        try (BufferedReader reader = new BufferedReader(new FileReader(new File(realPath)))) {
            String line = reader.readLine();
            while (line != null) {
                if (line.length() > 0 && !line.startsWith("@")) {
                    result.add(parser.parseLine(line));
                }
                line = reader.readLine();
            }
        } catch (Exception e) {
            e.printStackTrace();
        }

        return result;
    }
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的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.SAMFileHeader; //导入依赖的package包/类
public static SAMFileHeader getHeader(Configuration conf) {
	if (conf.get(BAM_HEADER_FILE_NAME) == null)
		return null;

	SAMFileHeader header = null;
	try {
		Path headerPath = new Path(conf.get(BAM_HEADER_FILE_NAME));
		HdfsHeaderLineReader reader = new HdfsHeaderLineReader(headerPath,
				conf);
		header = SamFileHeaderCodec.readHeader(reader);
	} catch (IOException e) {
		throw new RuntimeException(e.getMessage());
	}

	return header;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
/**
 *  constructor
 * @param options
 * @param samFileHeader
 */
public VariantCallingEngine(GenotyperOptions options, SAMFileHeader samFileHeader) {
    this.options = options;
    genomeLocationParser = new GenomeLocationParser(samFileHeader.getSequenceDictionary());
    samples = new HashSet<>();
    for(SAMReadGroupRecord rg : samFileHeader.getReadGroups()) {
        samples.add(rg.getSample());
    }

    GenotypeLikelihoodCalculator.getGenotypeLikelihoodsCalculatorObject(options);
    GenotypeLikelihoodCalculator.getCalculators(options);

    this.N = samples.size() * options.getSamplePloidy();
    log10AlleleFrequencyPriorsSNPs = new double[N+1];
    log10AlleleFrequencyPriorsIndels = new double[N+1];
    computeAlleleFrequencyPriors(N, log10AlleleFrequencyPriorsSNPs, options.getHeterozygosity());
    computeAlleleFrequencyPriors(N, log10AlleleFrequencyPriorsIndels, options.getIndelHeterozygosity());
    filter.add(LOW_QUAL_FILTER_NAME);
    tracker = new VariantDataTracker();

    annotationEngine = new VariantAnnotatorEngine(options.getAnnotationGroups(), options.getAnnotations(), null);
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
private void initialize(String input, SAMFileHeader header) {
	tables = new HashMap<String, RecalibratorReportTable>();
	addTables(input);
	option.parse(tables.get(RecalibratorUtil.ARGUMENT_TABLE_NAME));

	covariates = CovariateUtil.initializeCovariates(option, header);
	for (int i = 2; i < covariates.length; i++) {
		String covName = covariates[i].getClass().getSimpleName().split("Covariate")[0];
		optionalIndex.put(covName, i - 2);
	}

	recalTable = new RecalibratorTable(covariates,
			readGroupSize(tables.get(RecalibratorUtil.RECALIBRATOR_TABLE_NAME[0])));
	// read group table parse
	readGroupParser(tables.get(RecalibratorUtil.RECALIBRATOR_TABLE_NAME[0]),
			recalTable.getTable(RecalibratorTable.Type.READ_GROUP_TABLE));
	
	// quality group table parse
	qualityScoreParser(tables.get(RecalibratorUtil.RECALIBRATOR_TABLE_NAME[1]),
			recalTable.getTable(RecalibratorTable.Type.QUALITY_SCORE_TABLE));
	// covariate tables parse
	covariateParser(tables.get(RecalibratorUtil.RECALIBRATOR_TABLE_NAME[2]), recalTable);
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
private static void writeToSlice(@NotNull final String path, @NotNull final SAMFileHeader header,
        @NotNull final CloseableIterator<SAMRecord> iterator) {
    final File outputBAM = new File(path);
    final SAMFileWriter writer = new SAMFileWriterFactory().setCreateIndex(true).makeBAMWriter(header, true, outputBAM);
    String contig = "";
    while (iterator.hasNext()) {
        final SAMRecord record = iterator.next();
        if (record.getContig() != null && !contig.equals(record.getContig())) {
            contig = record.getContig();
            LOGGER.info("Reading contig: {}", contig);
        }
        writer.addAlignment(record);
    }
    iterator.close();
    writer.close();
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
@NotNull
public static List<ChromosomeLength> create(@NotNull final SAMFileHeader header, @NotNull final Set<String> excludedSequences) {
    final List<ChromosomeLength> results = Lists.newArrayList();

    for (final SAMSequenceRecord samSequenceRecord : header.getSequenceDictionary().getSequences()) {
        final String chromosome = samSequenceRecord.getSequenceName();
        if (!excludedSequences.contains(chromosome)) {
            results.add(ImmutableChromosomeLength.builder()
                    .chromosome(chromosome)
                    .length(samSequenceRecord.getSequenceLength())
                    .build());
        }
    }

    return results;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
private static File queryNameSortedBAM(final SamReader reader, final QueryInterval[] intervals, final String name) throws IOException {

        final SAMFileHeader header = reader.getFileHeader().clone();
        header.setSortOrder(SAMFileHeader.SortOrder.queryname);

        final File file = File.createTempFile(name, ".bam");
        final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(header, false, file);

        final SAMRecordIterator iterator = reader.queryOverlapping(intervals);
        while (iterator.hasNext()) {
            writer.addAlignment(iterator.next());
        }

        iterator.close();
        writer.close();

        return file;
    }
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
public static int findFirstOverlappingRegion(SAMFileHeader samHeader, SAMRecordWrapper read, int readStart, int readEnd, List<Feature> regions, int start) {
	if (start < 0) {
		start = 0;
	}

	for (int idx=start; idx<regions.size(); idx++) {
		Feature region = regions.get(idx);
		if ( (read.getSamRecord().getReferenceIndex() < samHeader.getSequenceDictionary().getSequenceIndex(region.getSeqname())) ||
			 (read.getSamRecord().getReferenceName().equals(region.getSeqname()) && readStart < region.getStart()) ) {
			
			// This read is in between regions
			// TODO: adjust start region here
			return -1;
		} else if (region.overlaps(read.getSamRecord().getReferenceName(), readStart, readEnd)) {
			return idx;
		}
	}
	
	// This read is beyond all regions
	return -1;
}
 

import htsjdk.samtools.SAMFileHeader; //导入依赖的package包/类
public static List<Integer> findAllOverlappingRegions(SAMFileHeader samHeader, SAMRecordWrapper read, List<Feature> regions, int start) {
	List<Integer> overlappingRegions = new ArrayList<Integer>();
	
	for (Span span : read.getSpanningRegions()) {
	
		int idx = findFirstOverlappingRegion(samHeader, read, span.start, span.end, regions, start);
		if (idx > -1) {
			overlappingRegions.add(idx);
			boolean isOverlap = true;
			idx += 1;
			
			while (isOverlap && idx < regions.size()) {
				Feature region = regions.get(idx);
				if (region.overlaps(read.getSamRecord().getReferenceName(), span.start, span.end)) {
					overlappingRegions.add(idx);
				} else {
					isOverlap = false;
				}
				idx += 1;
			}
		}
	}
	
	return overlappingRegions;
}