1
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ Can you please explain why you want to do this? It's a very unusual situation, and it's likely that there are better ways to approach what you want to do. It is often better to ask questions about your actual goal, rather than intermediate steps (see here for more information). $\endgroup$
    – gringer
    Jun 27 at 6:18

2 Answers 2

4
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ I'll add slightly to this excellent answer and mention that if this is public data (highly likely), there is almost certainly a link somewhere to download the underlying fastq files. $\endgroup$
    – Devon Ryan
    Jun 26 at 19:37
  • 1
    $\begingroup$ Excellent answer that misuses the word "convert" like most people on these forums do. One does not convert BAM to (big)wig, one extracts/computes wig information from a BAM file. "Convert" kind of implies information content preservation (convert sam to bam, etc.), which as you say in your answer is not true when getting WIG from BAM. $\endgroup$
    – Ram RS
    Jul 5 at 21:27
1
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.