используя htsjdk определенных классов от Jython или заводной

Question

Я пытаюсь методы доступа, предоставляемых htsjdk.j здесь: 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 2.7.0 после ввода в Jython реестра:

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.

Но я потерял, когда дело доходит до подкласса классов Java в Jython. Задача состоит в том, чтобы получить список блоков (ов) файла BAM, соответствующих заданной геномной позиции. Для тестового BAM/Bai файлов, то есть для CHRM 10000-15000 (хромосомные, start_pos, end_pos) Я получаю 11 сопоставляются чтения, используя с полки samtools автономную программу вместо htsjdk:

samtools view index_test.bam chrM:10000-15000 

Большое спасибо за вашу помощь

Дарек

Edit: повторно Отлич- 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() 

Выше код работает в хорошем состоянии. Моя проблема заключается не в том, что этот код не работает, но я не знаю, как правильно обращаться к методам из классов Java, касающихся формата файла BAI. Я искал AbstractBAMFileIndex в htsjdk/src/java/htsjdk/samtools/* java-файлах (git clone из repo @ github), но до сих пор неясно, что мне нужно делать.

Answer 1

Существует один пример на их 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 

} 

Answer 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