How to detect bedtools version used by pybedtools (in order to correctly preserve strand information when merging gtf records)

Question

I have issues working with pybedtools due to strand information being lost after extracting and merging transcript coordinates from a gtf file.

It seems that the solution to preserve strand information depends on the version of bedtools, so I would like to have my code detect bedtools version and adapt its behaviour accordingly.

On a machine with bedtools v2.27.0:

ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "cTel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "cTel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

merged = my_genes.sort().merge(s=True).saveas()
print(merged)

I       4118    10230
X       17714967        17718717

(Missing strand information)

merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

I       4118    10230   -
X       17714967        17718717        +

On a machine with bedtools v2.26.0:

ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_
source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "c
Tel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "c
Tel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

merged = my_genes.sort().merge(s=True).saveas()
print(merged)

I       4118    10230   -
X       17714967        17718717        +

merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

I       4118    10230   -       -
X       17714967        17718717        +       +

(Duplicate strand information)

Given that I work on linux systems, how can I detect bedtools version from within python?

Is there a better way to do what I want?

(Ideally, I would like to have my resulting bed consist in chrom, start, stop, name, score, strand fields, where the score would just be "." and the name would be the content of itvl["gene_id"])

---

Edit (05/05/2021): Regarding my last comment above, I found out that it is easy to convert gff into bed using pybedtools.featurefuncs.gff2bed (tested with bedtools 2.30):

from pybedtools import BedTool
from pybedtools.featurefuncs import gff2bed

ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).each(gff2bed, "gene_id").saveas()
print(my_genes)

I       4118    10230   WBGene00022277  .       -
X       17714967        17718714        WBGene00007068  .       +
X       17714967        17718717        WBGene00007068  .       +

merged = my_genes.sort().merge(s=True, c=[4, 5, 6], o="distinct").saveas()
print(merged)

I       4118    10230   WBGene00022277  .       -
X       17714967        17718717        WBGene00007068  .       +

This ensures proper bed format.

However, with bedtools 2.26, s=True makes strand information appear as a 4th column, so it is more complicated to obtain correct column order after merging.

Answer 1

pybedtools assumes that bedtools is in your path and bedtools itself will return the version with bedtools --version. So:

import subprocess
subprocess.check_output(['bedtools', '--version'], text=True).strip().split()[1]

For me that returns 'v2.30.0'.