基因注释

// 加载基因注释文件,并构建区间树
final OverlapDetector<Gene> geneOverlapDetector = GeneAnnotationReader.loadAnnotationsFile(ANNOTATIONS_FILE, bamDict);
// 打开输出bam文件
SAMFileWriter writer= new SAMFileWriterFactory().makeSAMOrBAMWriter(header, true, OUTPUT);
// 遍历输入bam中的每条read
for (SAMRecord r: inputSam) {
    // 对完成比对的read,进行find overlapping操作并添加注释
    if (!r.getReadUnmappedFlag())
        r=setAnnotations(r, geneOverlapDetector);
    // 输出注释后的read
    writer.addAlignment(r);
}

// 用来解析gtf每行数据,提取需要的字段
final FilteringGTFParser parser = new FilteringGTFParser(gtfFlatFile);
// gene name相同的gtf行,代表它们是一个gene的数据,使用它们构建GeneFromGTF数据类型,
// 此类型继承自Interval类型
final GeneFromGTFBuilder geneBuilder = new GeneFromGTFBuilder(parser);
// 初始化interval tree
final OverlapDetector<GeneFromGTF> overlapDetector = new OverlapDetector<>(0, 0);
while (geneBuilder.hasNext()){
        // 将每条gene添加到interval tree,其内部按照chromosome进行分类
        GeneFromGTF gene = geneBuilder.next();
        overlapDetector.addLhs(gene, gene);
    }

public SAMRecord setAnnotations (final SAMRecord r, final OverlapDetector<Gene> geneOverlapDetector) {
        Map<Gene, LocusFunction> map = AnnotationUtils.getInstance().getLocusFunctionForReadByGene(r, geneOverlapDetector);
        Set<Gene> exonsForRead = AnnotationUtils.getInstance().getConsistentExons (r, map.keySet(), ALLOW_MULTI_GENE_READS);

        List<Gene> genes = new ArrayList<>();

        for (Gene g: exonsForRead) {
            LocusFunction f = map.get(g);
            if (f==LocusFunction.CODING || f==LocusFunction.UTR)
                genes.add(g);
        }

        List<LocusFunction> allPassingFunctions = new ArrayList<>();
        if (USE_STRAND_INFO) {
            // constrain gene exons to read strand.
            genes = getGenesConsistentWithReadStrand(genes, r);
            // only retain functional map entries that are on the correct strand.
            for (Gene g: map.keySet()) {
                boolean strandCheck=readAnnotationMatchStrand(g, r);
                if (strandCheck) allPassingFunctions.add(map.get(g));
            }
        } else
            allPassingFunctions=new ArrayList<>(map.values());

        // if strand tag is used, only add locus function values for passing genes.
        for (Gene g: genes)
            allPassingFunctions.add(map.get(g));

        LocusFunction f = AnnotationUtils.getInstance().getLocusFunction(allPassingFunctions, false);

        if (genes.size()>1 && this.ALLOW_MULTI_GENE_READS==false)
            log.error("There should only be 1 gene assigned to a read for DGE purposes.");

        String finalGeneName = getCompoundGeneName(genes);
        String finalGeneStrand = getCompoundStrand(genes);

        if (f!=null)
            r.setAttribute(this.FUNCTION_TAG, f.toString());
        if (finalGeneName!=null && finalGeneStrand!=null) {
            r.setAttribute(this.TAG, finalGeneName);
            r.setAttribute(this.STRAND_TAG, finalGeneStrand);
        } else {
            r.setAttribute(this.TAG, null);
            r.setAttribute(this.STRAND_TAG, null);
        }
        return (r);
    }

private LocusFunction getLocusFunctionForRead (final SAMRecord rec, final Gene g) {
    List<AlignmentBlock> alignmentBlocks = rec.getAlignmentBlocks();

    LocusFunction [] blockSummaryFunction = new LocusFunction[alignmentBlocks.size()];
    Set<Gene> temp = new HashSet<>();
    temp.add(g);

    for (int i=0; i<alignmentBlocks.size(); i++) {
        AlignmentBlock alignmentBlock =alignmentBlocks.get(i);

        LocusFunction [] blockFunctions=getLocusFunctionsByBlock(alignmentBlock, temp);
        LocusFunction blockFunction = getLocusFunction(blockFunctions, false);
        blockSummaryFunction[i]=blockFunction;
    }

    LocusFunction readFunction = getLocusFunction(blockSummaryFunction, false);
    return readFunction;
}