This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Hard Masked + Raw Fasta To Soft Masked Fasta

Hello,

I downloaded from ENSEMBL either unmasked genomic sequences or hard masked (Ns for repeats) genomic sequences. There is no soft masked sequence in sight on ENSEMBL ftp server.

Flipping upper cases to lower cases for each position is quite straightforward (see below), but it takes long time for mammalian genome. Simple question: how to do it faster, be it in Python or any other language? 

#!/usr/bin/env python

from pyfasta import Fasta

masked_fasta = Fasta('test10k.rm.fa')
unmask_fasta = Fasta('test10k.fa')


for seqid in unmask_fasta.keys():
  print ">" + seqid

  unmasked_seq = unmask_fasta[seqid]
  masked_seq   = masked_fasta[seqid]

  output_seq = ""
  for position in range(0, len(unmasked_seq)):
    if masked_seq[position] == "N":
      base =  unmasked_seq[position].lower()
    else:
      base  =  unmasked_seq[position]
    output_seq += base
  print output_seq
python code
• 5,197 views • link
updated by Gustavo Gustavo • written by Darked89 Darked89

2 answers

pyfasta lets you turn the sequence into a numpy array so you can avoid loops. I think this is close to what you want to do:

from pyfasta import Fasta
import numpy as np
import sys

fa = Fasta(sys.argv[1]) # sequence with N's
fb = Fasta(sys.argv[2]) # sequence with replacements

for seqid, aseq in fa.iteritems():

    aseq = np.array(aseq)
    bseq = np.array(str(fb[seqid]).lower(), dtype="c")

    cseq = np.where(aseq == "N", bseq, aseq)

    print ">%s\n%s" % (seqid, cseq.tostring())
• 0 views • link
written by brentp brentp

This really speed things up. I went from ca 4hrs to 40mins. Thank you!

• 0 views • link
written by Darked89 Darked89

Slightly sidestepping your question, and assuming your genome of interest is among those provided...

--> http://repeatmasker.org/PreMaskedGenomes.html

• 1 views • link
written by Gustavo Gustavo

Thanks, indeed it is there, but it is an older version.

• 0 views • link
written by Darked89 Darked89

Log in to answer this question.

More posts like this