Converting .bw files to .fq (fastq)

Question

I am able to convert .bw files to .fq(fastq) manually in Shell, but I would like to automate the process coz I have hundreds of .bw file that I need to convert. Till now, I could think of following script, but unable to solve the problem. The .bw files are in H3K27ac folder.

import subprocess

path_bw = "~/Documents/H3K27ac"
import os

files = os.path.expanduser(path_bw)
for file in files:
    if file.endswith('.gz'):
        files.remove(file) 
files = [x[:-3] if x.endswith('.bw') else x for x in files]


# for .bw files

for file in files:
    bigwigtowig = "bigWigToWig " + path_bw + file + ".bw " + path_bw + file +".wig"
    subprocess.call(bigwigtowig, shell=True)

    wigtobed = "wig2bed < " + path_bw + file + ".wig > " + path_bw + file + ".bed"
    subprocess.call(wigtobed, shell=True)

    getfasta = "bedtools getfasta -fi " + path_bw + "hg19.fa -bed " + path_bw + file + ".bed -fo " + path_bw + file + ".fa"
    subprocess.call(getfasta, shell=True)

    perl = "perl " + path_bw + "fasta_to_fastq.pl" + path_bw + file + ".fa > " + path_bw + file + ".fq"
    #os.system(perl)
    subprocess.run(perl, shell=True)

Manually, I used the following to convert the .bw files:

1. # bigWigToWig file1.bw fileout.wig
2. # wig2bed <fileout.wig> fileout.bed
3. # bedtools getfasta -fi hg19.fa -bed fileout.bed -fo fileout.fa
4. # perl fasta_to_fastq.pl fileout.fa > file1.fq

All the above packages were preinstalled via shell and works perfectly. hg19 is the genom that i am using. fasta_to_fastq.pl is openly available in github and helps to convert .fa to .fq.

MAin thing i require is to convert all the files from .bw to .fq and this is the only way in my research that i could find. I need all the names to be the same as output. For example file1.bw shoud be file1.fq , file2.bw --> file2.fq and so on.

Do you have any ideas if this is possible or are there easier ways to convert it? It would be awesome to know your view point.

Answer 1

What you are doing makes no sense at all. A bigwig file contains intervals and then indicates the coverage of that interval. That coverage again comes from the alignments covering that region. There is no direct connection between the interval and the fastq file other than that the fastq file is the basis for the coverage calculation. Converting the interval to a fasta and then forcing it into a fastq file is technically possible, but has absolutely no practical meaning. You cannot recover sequencing reads from a bigwig, that information is lost once you convert your alignments (usually a BAM file) into a bigwig. You should open a new question describing the problem that this is actually based on.

Answer 2

The code is verbose and I ain't at all sure some of it will work, at least in two places.

The first point is os.path.expanduser ... admittedly I prefer pathlib, however that is definitely not how I understand that command. My understanding of expanduser is if ~ is being used. You definitely want to look at globb-ing for the opening part of the code.

from pathlib import Path
import os

mypath = Path(os.path.expanduser('~/Documents/H3K27ac'))
files = mypath.glob('*.bw') # or use ('*' + '.bw') 

I'm pretty certain that is what you want and its a lot easier than lines 1 to 10 in your code. The files list will then go straight into the loop.

Also, please check

files = [x[:-3] if x.endswith('.bw') else x for x in files]

Comprehensions are never logical in Python but convention is always to place the if statement after the list of interest.

Thus the following is much closer:

files = [x for x in files if not x.endswith('.bw') else x[:-3]]

Note not tested, so this may throw a bug, but its better. However, using glob bing you don't need this statement.

The other issue is the sub.process call ... Normally they are written:

os.chdir(path_bw)
file1 = 'file1.bw' 
file2 = 'fileout.wig'
subprocess.call(["bigWigToWig",  file1, file2], shell=True)

The way you've written is:

• Really hard to read
• I don't think it would work because of the absence of the [], certainly if you are using popen it will not work

Also check the following variant:

subprocess.call(["bigWigToWig",  file1, file2], text=True)

I know its seems a bit weird, but it very likely work (don't ask why)

I agree subprocess.call is a good idea in a Python pipeline. Issuing a hard path within subprocess call I'm not a fan hence using os to get Python to move into the cwd is a better approach.

A further point is the way you've used 'file' throughout your pipeline line, hmmmm... okay, but its difficult to read. Its not generally how a pipeline is written, but just saying.

nextflow is an increasingly popular pipeline language and you may consider learning that rather than going down a Python route (Python is very good for data science), of course this particularly question could be very easily solved in bash (but for some reason it ain't trendy).