2

我正在尝试从这里访问 htsjdk.jar 提供的方法： https ://samtools.github.io/htsjdk/

并记录在这里： https ://samtools.github.io/htsjdk/javadoc/htsjdk/index.html

使用 jython。我需要访问/查询 BAM 文件索引（BAI 文件）以获取二进制 BAM 文件中的起始位置的方法。测试 BAM 和 BAI 文件可以从以下网址获得： https ://github.com/samtools/htsjdk/tree/master/testdata/htsjdk/samtools/BAMFileIndexTest

在放入 Jython 注册表后的 jython 2.7.0 中：

python.security.respectJavaAccessibility = false
#I did in the jython comandline:
import sys
sys.path.append("/usr/local/soft/picard_1.138/htsjdk-1.138.jar")
from htsjdk.samtools import *
from java.io import File
#the BAM index file + BAM files 
bai_fh = File("./index_test.bam.bai")
mydict = SAMSequenceDictionary()
bai_foo = DiskBasedBAMFileIndex(bai_fh, mydict)

我可以访问一些方法，例如 bai_foo.getNumberOfReferences() 等，但所需的方法
getBinsOverlapping(int referenceIndex, int startPos, int endPos) 在 BrowseableBAMIndex 接口中。

但是当谈到在 Jython 中对 Java 类进行子类化时，我迷失了方向。任务是获取与给定基因组位置相对应的 BAM 文件块列表。对于测试 BAM/BAI 文件，即对于 chrM 10000-15000（染色体、start_pos、end_pos），我使用现成的 samtools 独立程序而不是 htsjdk 获得 11 个映射读取：

samtools view index_test.bam  chrM:10000-15000

非常感谢您的帮助

达雷克

编辑：重新成为 groovy 部分 Groovy 版本：2.4.4

groovy -cp libs/htsjdk-1.138.jar test_htsjdk.groovy

#!/usr/bin/env groovy

import htsjdk.samtools.*

File bam_fh =  new File("./A.bam")
File bai_fh =  new File("./A.bam.bai")

def mydict = new SAMSequenceDictionary()
def bai_foo = new DiskBasedBAMFileIndex(bai_fh, mydict)
println bai_foo.getNumberOfReferences()

上面的代码在 groovy 中工作。我的问题不是这段代码不起作用，而是我不知道从处理 BAI 文件格式的 Java 类中访问方法的正确方法。我确实在 htsjdk/src/java/htsjdk/samtools/*java 文件（来自 repo@github 的 git clone）中搜索了 AbstractBAMFileIndex，但仍然不清楚我需要做什么。

4

2 回答 2

3

他们的github上有一个示例，我附上了一小部分，但还有更多示例说明如何使用他们的 SamReader

    /**
     * Broken down
     */
    final SamReaderFactory factory =
            SamReaderFactory.makeDefault().enable(SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS).validationStringency(ValidationStringency.LENIENT);

    final SamInputResource resource = SamInputResource.of(new File("/my.bam")).index(new URL("http://broadinstitute.org/my.bam.bai"));

    final SamReader myReader = factory.open(resource);

    for (final SAMRecord samRecord : myReader) {
        System.err.print(samRecord);
    }

无论格式如何，Groovy 都非常适合与一堆文件进行交互。

new File('my/directory/with/many/files').eachFile{File readFile ->
    final SamInputResource resource = SamInputResource.of(readFile).index(new URL('IGotLazy.bai'))
    final SamReader myReader = ... etc

}
于 2015-09-09T22:25:03.950 回答
2

我不知道要为参考索引放什么，但作为一个从不处理此类数据的人，打印输出对我来说已经足够好了：

import sys
sys.path.append("/home/shackle/sam-tools/samtools-htsjdk-f650176/dist/htsjdk-1.138.jar")
from htsjdk.samtools import *
from java.io import File
#the BAM index file + BAM files 
indexFile = File("/home/shackle/index_test.bam.bai")
bamFile = File("/home/shackle/index_test.bam")
sfr = SAMFileReader(bamFile, indexFile)
sfr.enableIndexCaching(True)
bbi = sfr.getBrowseableIndex()
for i in range(0,256):
    print "i = " , i
    bl = bbi.getBinsOverlapping(i, 10000, 15000)
    count = 0
    for bin in bl:
        print "bin.getBinNumber() = " , bin.getBinNumber()
        count = count + 1
        print "count = ", count
于 2015-09-09T23:33:49.403 回答