Skip to content

Instantly share code, notes, and snippets.

@daler

daler/ex.py

Created August 9, 2017 14:17
Show Gist options
  • Save daler/9623fad6b3cfa1c599591b0538e9f4dd to your computer and use it in GitHub Desktop.
Save daler/9623fad6b3cfa1c599591b0538e9f4dd to your computer and use it in GitHub Desktop.
Download ZIP
bedtools issue 553
Raw
ex.py
from __future__ import print_function
import pybedtools
import subprocess as sp
print('pbt:', pybedtools.__version__)
pybedtools.debug_mode(True)
print(
sp.check_output(['which', 'bedtools'], universal_newlines=True).strip(),
'(',
sp.check_output(['bedtools', '--version'], universal_newlines=True).strip(),
') says:'
)
a = pybedtools.BedTool('variants_chr1.vcf.gz')
b = pybedtools.BedTool('reads_chr1.bam')
c = a.intersect(b, u=True, split=True, sorted=True)
for i in c:
print(i, end='')
print()
Raw
output.txt
===========================
Python = py2
BEDTools = dev
pbt: 0.7.10
pybedtools.logger [INFO]: Debug mode enabled. You may also want to set pybedtools.KEEP_TEMPFILES=True to prevent automatic deletion of files upon exit.
/tmp/tst/bedtools-dev/bin/bedtools ( bedtools v2.26.0-125-g52db6549 ) says:
pybedtools.logger [DEBUG]: BedTool.handle_kwargs() got these kwargs:
{'a': 'variants_chr1.vcf.gz',
'b': <BedTool(reads_chr1.bam)>,
'sorted': True,
'split': True,
'u': True}
pybedtools.logger [DEBUG]: helpers.call_bedtools(): input is filename, output is filename (/tmp/pybedtools.1ejW4r.tmp)
pybedtools.logger [DEBUG]: helpers.call_bedtools(): cmds=bedtools intersect -u -split -sorted -b reads_chr1.bam -a variants_chr1.vcf.gz
chr1 216 rs01 T G . . RS=1;RPOS=216
chr1 576 rs02 C A . . RS=2;RPOS=576
chr1 836 rs04 A C . . RS=4;RPOS=836
===========================
Python = py2
BEDTools = v2.26
pbt: 0.7.10
pybedtools.logger [INFO]: Debug mode enabled. You may also want to set pybedtools.KEEP_TEMPFILES=True to prevent automatic deletion of files upon exit.
/tmp/tst/bedtools-v2.26/bin/bedtools ( bedtools v2.26.0 ) says:
pybedtools.logger [DEBUG]: BedTool.handle_kwargs() got these kwargs:
{'a': 'variants_chr1.vcf.gz',
'b': <BedTool(reads_chr1.bam)>,
'sorted': True,
'split': True,
'u': True}
pybedtools.logger [DEBUG]: helpers.call_bedtools(): input is filename, output is filename (/tmp/pybedtools.h6plf4.tmp)
pybedtools.logger [DEBUG]: helpers.call_bedtools(): cmds=bedtools intersect -u -split -sorted -b reads_chr1.bam -a variants_chr1.vcf.gz
chr1 216 rs01 T G . . RS=1;RPOS=216
chr1 576 rs02 C A . . RS=2;RPOS=576
chr1 836 rs04 A C . . RS=4;RPOS=836
Traceback (most recent call last):
File "ex.py", line 16, in <module>
for i in c:
File "pybedtools/cbedtools.pyx", line 787, in pybedtools.cbedtools.IntervalIterator.__next__ (pybedtools/cbedtools.cxx:11123)
File "pybedtools/cbedtools.pyx", line 652, in pybedtools.cbedtools.create_interval_from_list (pybedtools/cbedtools.cxx:9208)
IndexError: list index out of range
===========================
Python = py3
BEDTools = dev
pbt: 0.7.10
pybedtools.logger [INFO]: Debug mode enabled. You may also want to set pybedtools.KEEP_TEMPFILES=True to prevent automatic deletion of files upon exit.
/tmp/tst/bedtools-dev/bin/bedtools ( bedtools v2.26.0-125-g52db6549 ) says:
pybedtools.logger [DEBUG]: BedTool.handle_kwargs() got these kwargs:
{'a': 'variants_chr1.vcf.gz',
'b': <BedTool(reads_chr1.bam)>,
'sorted': True,
'split': True,
'u': True}
pybedtools.logger [DEBUG]: helpers.call_bedtools(): input is filename, output is filename (/tmp/pybedtools.k8k_qt15.tmp)
pybedtools.logger [DEBUG]: helpers.call_bedtools(): cmds=bedtools intersect -u -split -sorted -b reads_chr1.bam -a variants_chr1.vcf.gz
chr1 216 rs01 T G . . RS=1;RPOS=216
chr1 576 rs02 C A . . RS=2;RPOS=576
chr1 836 rs04 A C . . RS=4;RPOS=836
===========================
Python = py3
BEDTools = v2.26
pbt: 0.7.10
pybedtools.logger [INFO]: Debug mode enabled. You may also want to set pybedtools.KEEP_TEMPFILES=True to prevent automatic deletion of files upon exit.
/tmp/tst/bedtools-v2.26/bin/bedtools ( bedtools v2.26.0 ) says:
pybedtools.logger [DEBUG]: BedTool.handle_kwargs() got these kwargs:
{'a': 'variants_chr1.vcf.gz',
'b': <BedTool(reads_chr1.bam)>,
'sorted': True,
'split': True,
'u': True}
pybedtools.logger [DEBUG]: helpers.call_bedtools(): input is filename, output is filename (/tmp/pybedtools.ysgkk0ri.tmp)
pybedtools.logger [DEBUG]: helpers.call_bedtools(): cmds=bedtools intersect -u -split -sorted -b reads_chr1.bam -a variants_chr1.vcf.gz
Raw
test.sh
#!/bin/bash
HERE=$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )
# build python 2 and python 3 environments
if [[ ! -e py2env ]]; then
conda create -y --prefix py2env pybedtools python=2 -c conda-forge -c bioconda
fi
if [[ ! -e py3env ]]; then
conda create -y --prefix py3env pybedtools python=3 -c conda-forge -c bioconda
fi
if [[ ! -e bedtools-v2.26/bin/bedtools ]]; then
git clone https://github.com/arq5x/bedtools2.git bedtools-v2.26
(cd bedtools-v2.26 && git checkout v2.26.0 && make)
fi
if [[ ! -e bedtools-dev/bin/bedtools ]]; then
git clone https://github.com/arq5x/bedtools2.git bedtools-dev
(cd bedtools-dev && make)
fi
wget -q -O- https://github.com/arq5x/bedtools2/files/1191844/variants_chr1.vcf.gz > variants_chr1.vcf.gz
wget -q -O- https://github.com/arq5x/bedtools2/files/1191857/reads_chr1.sam.txt > reads_chr1.sam
samtools view -Sb reads_chr1.sam > reads_chr1.bam
for py in 2 3; do
source activate py${py}env
for bt in dev v2.26; do
echo
echo "==========================="
echo "Python = py$py"
echo "BEDTools = $bt"
PATH="$HERE/bedtools-$bt/bin:$PATH" python ex.py
PATH="$HERE/bedtools-$bt/bin:$PATH" bedtools intersect \
-u -split -sorted -b reads_chr1.bam -a variants_chr1.vcf.gz \
> bedtools_results_$bt.vcf
done
source deactivate
done
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment