wig file from fasta
2
0
Entering edit mode
7.3 years ago
boczniak767 ▴ 870

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?

wig fasta • 2.8k views
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode
7.3 years ago

Normally one makes such tracks with GEM.

ADD COMMENT
0
Entering edit mode
7.3 years ago
boczniak767 ▴ 870

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)
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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"):'

ADD REPLY

Login before adding your answer.

Traffic: 1423 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6