pyBigWig
0.3.23
Расширение Python, написанное в C, для быстрого доступа к файлам Bigbid и доступа к и созданию файлов BigWig. Это расширение использует Libbigwig для локального и удаленного доступа к файлам.
Вы можете установить это расширение прямо из GitHub с:
pip install pyBigWig
или с кондой
conda install pybigwig -c conda-forge -c bioconda
Следующие непитоны должны быть установлены:
curl-config )Заголовки и библиотеки для них требуются.
Основное использование следующим образом:
>>> import pyBigWig
Это будет работать, если ваш рабочий каталог является каталогом исходного кода Pybigwig.
>>> bw = pyBigWig.open("test/test.bw")
Обратите внимание, что если файл не существует, вы увидите сообщение об ошибке, и None не будет возвращен. Будьте по умолчанию, все файлы открываются для чтения, а не написания. Вы можете изменить это, передавая режим, содержащий w :
>>> bw = pyBigWig.open("test/output.bw", "w")
Обратите внимание, что файл, открытый для написания, не может быть запрошен для его интервалов или статистики, его можно записать только . Если вы откроете файл для написания, вам затем нужно добавить заголовок (см. Раздел по этому поводу ниже).
Также поддерживается местный и удаленный доступ
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
В то время как вы можете указать режим для файлов Bigbid, он игнорируется. Объект, возвращенный pyBigWig.open() одинаково независимо от того, открываете ли вы файл Bigwig или Bigbed.
Поскольку могут быть открыты файлы Bigwig и Bigbid, может потребоваться определить, указывает ли данный объект bigWigFile на файл BigWig или Bigbid. Для этого можно использовать функции isBigWig() и isBigBed() :
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
Объекты bigWigFile содержат словарь, содержащий длину хромосомы, к которому можно получить доступ с помощью аксессуара chroms() .
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
Вы также можете напрямую запросить определенную хромосому.
>>> bw.chroms("1")
195471971L
Длина хранится «длинным» целым типом, поэтому есть суффикс L Если вы указываете несуществующую хромосому, то ничего не выводит.
>>> bw.chroms("c")
>>>
Иногда полезно напечатать заголовок Bigwig. Это представлено здесь как словарь Python, содержащий: версию (обычно 4 ), количество уровней масштабирования ( nLevels ), количество описанных оснований ( nBasesCovered ), минимальное значение ( minVal ), максимальное значение ( maxVal ), о Сумма всех значений ( sumData ) и сумма всех квадратных значений ( sumSquared ). Последние два из них необходимы для определения среднего и стандартного отклонения.
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
Обратите внимание, что это также может присутствовать для файлов с большим подн., И будут присутствовать те же клавиши словаря. Такие записи, как maxVal , sumData , minVal и sumSquared , в значительной степени не значимы.
Файлы BigWig используются для хранения значений, связанных с позициями и диапазонами. Как правило, мы хотим быстро получить доступ к среднему значению в диапазоне, что очень просто:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
Предположим, вместо среднего значения мы вместо этого хотели максимальное значение:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
Другими опциями являются «мин» (минимальное значение), «покрытие» (доля охватываемых оснований) и «Std» (стандартное отклонение значений).
Часто это тот случай, когда мы хотели бы вычислить значения некоторого количества равномерно распределенных бункеров в данном интервале, что также является простым:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins по умолчанию до 1, так же, как type по mean .
Если начальные и конечные позиции опущены, то используется вся хромосома:
>>> bw.stats("1")
[1.3351851569281683]
Примечание для Lay Reader: этот раздел довольно технический и включен только ради полноты. Резюме заключается в том, что если ваши потребности требуют точного среднего/макс/и т. Д. Сводные значения для интервала или интервалов, и что небольшой компромисс в скорости приемлем, чтобы вы должны использовать опцию
exact=Trueв функцииstats().
По умолчанию существуют некоторые неинтуитивные аспекты для вычисления статистики диапазонов в файле BigWig. Формат BigWig был первоначально создан в контексте браузеров генома. Там вычисление точной сводной статистики для данного интервала менее важна, чем быстро, чтобы вычислить приблизительную статистику (в конце концов, браузеры должны быть в состоянии быстро отобразить ряд смежных интервалов и поддержки прокрутки/масштабирования). Из-за этого файлы BigWig содержат не только интервальные ассоциации, но и sum of values / sum of squared values / minimum value / maximum value / number of bases covered для бункеров одинакового размера различных размеров. Эти разные размеры называются «уровнями масштабирования». Наименьший уровень масштабирования имеет контейнеры, которые в 16 раз превышают средний размер интервала в файле, и каждый последующий уровень масштабирования имеет контейнеры в 4 раза больше, чем в предыдущем. Эта методология используется в инструментах Кента и, следовательно, вероятно, используется почти в каждом существующем в настоящее время существующем файле BigWig.
Когда файл BigWig запрашивается для сводной статистики, размер интервала используется для определения того, использовать ли уровень масштабирования и, если это так, какой. Оптимальный уровень масштабирования - это тот, который имеет самые большие бункеры не более половины ширины желаемого интервала. Если такой уровень масштабирования не существует, вместо этого исходные интервалы используются для расчета.
Для согласованности с другими инструментами Pybigwig принимает эту же методологию. Однако, поскольку это (а) неинтуированное и (б) нежелательно в некоторых приложениях, Pybigwig позволяет вычислять точную сводную статистику независимо от размера интервала (то есть, это позволяет игнорировать уровни увеличения). Это было первоначально предложено здесь, и пример ниже:
>>> import pyBigWig
>>> from numpy import mean
>>> bw = pyBigWig.open("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign75mer.bigWig")
>>> bw.stats('chr1', 89294, 91629)
[0.20120902053804418]
>>> mean(bw.values('chr1', 89294, 91629))
0.22213841940688142
>>> bw.stats('chr1', 89294, 91629, exact=True)
[0.22213841940688142]
В то время как метод stats() может использоваться для извлечения исходных значений для каждого базового (например, путем установки nBins на количество оснований), предпочтительнее вместо этого использовать аксессу values() .
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
Изготовленный список всегда будет содержать одно значение для каждой базы в указанном диапазоне. Если конкретная база не имеет связанного значения в файле BigWig, возвращаемое значение будет nan .
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
Иногда удобно получить все записи, перекрывающие некоторый диапазон. Это можно сделать с помощью функции intervals() :
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
То, что возвращается, - это список кортежей, содержащих: начальное положение, конечное положение и значение. Таким образом, приведенный выше пример имеет значения 0.1 , 0.2 и 0.3 в положениях 0 , 1 и 2 соответственно.
Если начальная и конечная позиция опущена, то возвращаются все интервалы на указанной хромосоме:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
В отличие от файлов bigwig, файлы Bigbid содержат записи, которые представляют собой интервалы со связанной строкой. Вы можете получить доступ к этим записям, используя функцию entries() :
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
>>> bb.entries('chr1', 10000000, 10020000)
[(10009333, 10009640, '61035t130t-t0.026t0.42t404'), (10014007, 10014289, '61047t136t-t0.029t0.42t404'), (10014373, 10024307, '61048t630t-t5.420t0.00t2672399')]
Выход - это список входных кортежей. Элементы кортежа являются start и end позицией каждой записи, за которой следует связанная string . Строка возвращается точно так же, как она удерживается в большом файле, поэтому анализ ее остается для вас. Чтобы определить, какие различные поля находятся в этой строке, обратитесь к строке SQL:
>>> bb.SQL()
table RnaElements
"BED6 + 3 scores for RNA Elements data"
(
string chrom; "Reference sequence chromosome or scaffold"
uint chromStart; "Start position in chromosome"
uint chromEnd; "End position in chromosome"
string name; "Name of item"
uint score; "Normalized score from 0-1000"
char[1] strand; "+ or - or . for unknown"
float level; "Expression level such as RPKM or FPKM. Set to -1 for no data."
float signif; "Statistical significance such as IDR. Set to -1 for no data."
uint score2; "Additional measurement/count e.g. number of reads. Set to 0 for no data."
)
Обратите внимание, что первые три записи в строке SQL не являются частью строки.
Если вам нужно только знать, где находятся записи, а не связанные с ними значения, вы можете сохранить память, дополнительно указав withString=False in entries() :
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
Если вы открыли файл для написания, вам нужно будет дать ему заголовок, прежде чем вы сможете добавить какие -либо записи. Заголовок содержит все хромосомы по порядку и их размеры. Если ваш геном имеет две хромосомы, Chr1 и CHR2, длины 1 и 1,5 миллиона оснований, то следующее добавило бы соответствующий заголовок:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Заголовки BigWig чувствительны к случаям, поэтому chr1 и Chr1 разные. Аналогично, 1 и chr1 не совпадают, поэтому вы не можете смешивать имена хромосомы Ensembl и UCSC. После добавления заголовка вы можете добавить записи.
По умолчанию до 10 «уровней масштабирования» построены для файлов BigWig. Вы можете изменить этот номер по умолчанию с помощью необязательного аргумента maxZooms . Распространенным использованием этого является создание файла BigWig, который просто удерживает интервалы, и без уровня масштабирования:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
Если вы установите maxTooms=0 , обратите внимание, что IGV и многие другие инструменты не будут работать, поскольку они предполагают, что будет присутствовать по крайней мере один уровень масштабирования. Вам рекомендуется использовать по умолчанию, если вы не ожидаете, что файлы BigWig будут использоваться другими пакетами.
Предполагая, что вы открыли файл для написания и добавили заголовок, вы можете добавить записи. Обратите внимание, что записи должны быть добавлены по порядку, так как файлы BigWig всегда содержат упорядоченные интервалы. Есть три формата, которые файлы Bigwig могут использовать внутренне для хранения записей. Наиболее часто наблюдаемый формат идентичен файлу постельного графа:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
Эти записи будут добавлены следующим образом:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
Каждая запись занимает 12 байтов перед сжатием.
Второй формат использует фиксированный пролет, но размер переменного шага между записями. Они могут быть представлены в файле покачивания как:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
Приведенные выше записи описывают (1) позиции 501-520, 601-620 и 636-655. Это будет добавлено следующим образом:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
Каждая запись этого типа занимает 8 байтов перед сжатием.
В конечном итоге используется фиксированный шаг и пролет для каждой записи, соответствующий формату фиксированного штифта:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
Вышеуказанные записи описывают (1) основы 901-920, 931-950 и 961-980 и будут добавлены следующим образом:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
Каждая запись этого типа занимает 4 байта.
Обратите внимание, что Pybigwig постарается предотвратить добавление записей в неправильном порядке. Это, однако, требует дополнительных над головой. Если это не будет приемлемо, вы можете просто указать validate=False при добавлении записей:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
Вы, очевидно, тогда отвечаете за обеспечение того, чтобы вы не добавляли записи вне порядка. Полученные файлы в противном случае не будут использоваться.
Файл может быть закрыт с помощью простого bw.close() , как это обычно делается с другими типами файлов. Для открытых файлов для написания, закрытие файла записывает любые буферизованные записи на диск, конструирует и записывает индекс файла, и строит уровни масштабирования. Следовательно, это может занять немного времени.
Начиная с версии 0.3.0, Pybigwig поддерживает ввод координат с использованием целых чисел и векторов Numpy в некоторых функциях , если Numpy был установлен до установки Pybigwig . Чтобы определить, был ли Pybigwig был установлен с поддержкой Numpy, проверив numpy Accessor:
>>> import pyBigWig
>>> pyBigWig.numpy
1
Если pyBigWig.numpy составляет 1 , то Pybigwig был составлен с поддержкой Numpy. Это означает, что addEntries() может принимать координаты Numpy:
>>> import pyBigWig
>>> import numpy
>>> bw = pyBigWig.open("/tmp/delete.bw", "w")
>>> bw.addHeader([("1", 1000)], maxZooms=0)
>>> chroms = np.array(["1"] * 10)
>>> starts = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90], dtype=np.int64)
>>> ends = np.array([5, 15, 25, 35, 45, 55, 65, 75, 85, 95], dtype=np.int64)
>>> values0 = np.array(np.random.random_sample(10), dtype=np.float64)
>>> bw.addEntries(chroms, starts, ends=ends, values=values0)
>>> bw.close()
Кроме того, values() могут непосредственно выводить вектор Numpy:
>>> bw = bw.open("/tmp/delete.bw")
>>> bw.values('1', 0, 10, numpy=True)
[ 0.74336642 0.74336642 0.74336642 0.74336642 0.74336642 nan
nan nan nan nan]
>>> type(bw.values('1', 0, 10, numpy=True))
<type 'numpy.ndarray'>
Если у вас нет установленного сгиба, Pybigwig будет установлен без возможности доступа к удаленным файлам. Вы можете определить, сможете ли вы получить доступ к удаленным файлам с помощью pyBigWig.remote . Если это возвращает 1, то вы можете получить доступ к удаленным файлам. Если он возвращает 0, то вы не можете.
Начиная с версии 0.3.5, Pybigwig может читать и писать файлы BigWig, в которых отсутствуют записи. Обратите внимание, что такие файлы, как правило, не совместимы с другими программами, поскольку нет определения того, как должен выглядеть файл BigWig без записей. Для такого файла аксессу intervals() не будет возвращать None , функция stats() вернет список None из желаемой длины, а values() вернет [] (пустой список). Как правило, это должно позволять программам, использующим Pybigwig продолжать без проблем.
Для тех, кто хочет имитировать функциональность Pybigwig/Libbigwig в этом отношении, обратите внимание, что он рассматривает количество охватываемых баз (как сообщается в заголовке файла), чтобы проверить «пустые» файлы.
Файлы Wiggle, Bigwig и Bigbid используют 0 на основе полуоткрытых координат, которые также используются этим расширением. Таким образом, чтобы получить доступ к значению для первой базы на chr1 , можно было бы указать начальную позицию как 0 и конечную позицию как 1 . Точно так же базы от 100 до 115 будут иметь начало 99 и конец 115 . Это просто ради согласованности с базовым файлом BigWig и может измениться в будущем.
Pybigwig также доступен в качестве пакета в Galaxy. Вы можете найти его в инструментах, и IUC в настоящее время принимает определение XML на GitHub.
