wig file from fasta

Question

[boczniak767]

Hi, I need to produce fixed-step (1bp) wig file for use with F-seq.

bffBuilder takes as input a fixed step, step=1 wig file that contains a header and one value on each succeeding line for each base. For assemblies with multiple chromosomes, each chromosome should have a bff file generated separately. Each line of the wig file should contain a value between 0 and 1 that reflects the uniqueness of sequences aligning to that base on either strand. For example, if the sequence tag length is 35 bases, and the 35bp sequence starting at that position on the positive strand occurs 3 times in the genome, and the 35bp sequence starting on the negative strand occurs 5 times in the genome, this could be assigned the value 1/8 or 0.125.

How could it be done?

[boczniak767]

Thanks, I've tried it before, there was a number of issues with running its scripts.

I've tried it (both version GEM-binaries-Linux-x86_64-core_i3-20130406-045632 and GEM-binaries-Linux-x86_64-core_i3-20121106-022124) again. On Intel(R) Xeon(R) with 24 cores it gives me this error

Creating sequence and location files...Fatal error: exception Failure("Command 'echo "-------gem-indexer_fasta2meta+cont" > genome_idxgem.log && gem-indexer_fasta2meta+cont -i /home/sekwoja/agpv3/Zea_mays_Ensembl_AGPv3/Zea_mays/Ensembl/AGPv3/Sequence/WholeGenomeFasta/genome.fa -c dna --filter-function iupac-dna --st

So I've tried it on my second machine with AMD Phenom(tm) II X4 965. This time I've got

memory security violation

(or something like that, its my translation from Polish) error.

I've used the simplest syntax gem-indexer -i /home/mj/bin/data/genomes/zea_b73/genome.fa -o genome_gemidx but it gives the same error as commands with T and c options set.

So if there are alternative to GEM? It is abadoned and even old versions are beta.

[GenoMax]

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized.

This comment should have gone against @Devon's answer below.

Answer 1

[Devon Ryan]

Normally one makes such tracks with GEM.

Answer 2

[boczniak767]

I'ts turned-out that gem works, when all executables are in the the $PATH. However, gem-mappability and gem-2-wig can't create one-step wig file. One have to provide k-mer length and the file generated is variable step. So running:

gem-indexer -T 22 -i ~/agpv3/Zea_mays_Ensembl_AGPv3/Zea_mays/Ensembl/AGPv3/Sequence/Chromosomes/1.fa -o chr1_idxgem
gem-mappability -T 22 -I chr1_idxgem.gem -l 20 -o mapp_gem
gem-2-wig -I chr1_idxgem.gem -i chr1mapp.mappability -o chr1mapp

produce file like this:

variableStep    chrom=1 span=32
1       1
variableStep    chrom=1
33      0.0116169
variableStep    chrom=1 span=10
34      1
variableStep    chrom=1
44      0.00613059
variableStep    chrom=1 span=8
45      1

Which is not acceptable by bffBuilder ~/bin/bffBuilder/bin/bffBuilder chr1mapp.wig

Exception in thread "main" java.io.IOException: Bad '.wig' format for file /home/sekwoja/agpv3/gem_lipiec17/chr1mapp.wig
        at edu.duke.igsp.bffWriter.wigReader.badFile(wigReader.java:102)
        at edu.duke.igsp.bffWriter.wigReader.read(wigReader.java:34)
        at edu.duke.igsp.bffWriter.bffWriterApp.main(bffWriterApp.java:33)

[Devon Ryan]

You can reformat that with a bit of python/perl/whatever to make it fixed step.

[boczniak767]

Unfortunately, my programming skills will not suffice, is there any popular program which could do that?

[Devon Ryan]

This is such an obscure thing to want to do that there won't be any program for it. The following python code should work:

chrom = None
for line in open("foo.wig"):
    cols = line.strip().split()
    if cols[0] == "variableStep":
        span = 1
        for col in cols[1:]:
            if col.startswith("chrom="):
                if chrom is None or chrom != col[6:]:
                    print("fixedStep chrom={} start=1".format(col[6:]))
                chrom = col[6:]
            elif col.startswith("span="):
                span = int(col[5:])
    else:
        for i in xrange(span):
            print("{}".format(cols[1]))

Make sure to rename foo.wig to your file name and redirect the output to a new file (it's printed to the screen otherwise).

[boczniak767]

I'm somehow returning to this data, and tried your script. It gives error (my translation from Polish)

./wig1step_devonryan.py: line 1: chrom: command not found ./wig1step_devonryan.py: line 2: syntax error in unexpected mark (' ./wig1step_devonryan.py: line 2:for line in open("Zea_mays_mappa125.wig"):'