0
$\begingroup$

I am working on a project and am having issues with the following code that I have written in nano:

from Bio import SeqIO
import sys
import re 


    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
       print (">"+list.group(1)) 

This is the error I receive when I execute my command on command-line:

File "mpo.py", line 7, in <module>
    gene_id = myfile.id
NameError: name 'myfile' is not defined

I have a fasta file with the format

>H149xcV_Fge342_r3_h2_d1 len=210 path=[0:0-206]
ACTATACATGAGGAGAACATAGAACAAAAATGGGACCATAGATATATAACAATAGAAGATATAGAGAACACAATAGACAACTTATTAGGAAAGAGGTGTGTCGTCATGGAGCTGATGTTCGAGGATACTTTGCATGGTCATTCTTGGATAATTTTGAGTGGGCTATGGGATACACCAAGAGGTTTGGCATTGTTTATGTTGATTATAAGAACGGGC

>H149xcV_ytR1oP_r3_h2_d1 len=306 path=[0:0-207]
ATTAGAGTCTGAGAGAGTCTTGATTTGTCGTCGTCGAGAAATATAGGAGATCTGATTAGAGGAGAGAGCGGCCTAGGCGATGCGCGATATAGCGCTATATAGGCCTAGAGGAGAGTCTCTCTCTTTTAGAAGAGATAATATATATATATATATGGCTCTCCGGCGGGGCCGCGCGAGAGCTCGATCGATCGATATTAGCTGTACGATGCTAGCTAGCTTATATTCGATCGATTATAGCTTAGATCTCTCTCTAAAGGTCGATATCGCTTATGCGCGCGTATATCG

I wish to reformat my file so that it only provides me with the unique gene id's for genes with a length greater than 250 bp.

My desired output is:

>H149xcV_Fge342_r3_h2
>H149xcV_ytR1oP_r3_h2
>H149xcV_DPN78333_r3_h2
>H149xcV_AgV472_r3_h2
>H149xcV_DNP733_r3_h2
$\endgroup$
3
  • $\begingroup$ note that after your edit, your code now does not match your error. You should update the error to the new error, otherwise we can't help you very well. $\endgroup$ Nov 23 '20 at 20:13
  • $\begingroup$ Hi @MaximilianPress I have updated the error! $\endgroup$ Nov 23 '20 at 20:18
  • $\begingroup$ I updated the post to attempt to trap the bugs. $\endgroup$
    – M__
    Nov 23 '20 at 20:49
1
$\begingroup$

You might consider changing the code to make the whitespace a little more consistent (as it was, it failed when i tried to run it).

Also, why are you using regular expressions? You have already extracted the gene_id using Biopython. The following works for me:

(updated to handle the underscore delimiting)

from Bio import SeqIO
import sys
import re 

fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile.seq) > 250:  # note .seq here!
        gene_id = myfile.id
        # unclear why using regex
        #match_list = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
        #print (">"+match_list.group(0)) 
        print("_".join(gene_id.split("_")[0:-1]))
$\endgroup$
9
  • $\begingroup$ Hi @MaximilianPress , how would I change the regex so that it only provides me with >H149xcV_Fge342_r3_h2 and not>H149xcV_Fge342_r3_h2_d1. I would not like the very end part that consists of d1. So basically, how can I format it so it it only provides this >H149xcV_Fge342_r3_h2. If I am using regex would I not need to use gene.id? $\endgroup$ Nov 23 '20 at 20:43
  • $\begingroup$ @andnowmywatchbegins wants to remove the stuff after the space in the fasta header. In phylogeny this is either removed or spaces substituted because the space in phylip format denotes the start of the sequence $\endgroup$
    – M__
    Nov 23 '20 at 21:13
  • $\begingroup$ @Michael I've tried reformatting the entire regular expression and I am still getting the full >H149xcV_Fge342_r3_h2_d1, is there a way to only receive >H149xcV_Fge342_r3_h2 as in output? Thanks! $\endgroup$ Nov 23 '20 at 21:20
  • 1
    $\begingroup$ Can you put this as a separate question please? .. Yes you will but thats a different part of regex $\endgroup$
    – M__
    Nov 23 '20 at 21:23
  • $\begingroup$ @Michael stuff after the whitespace is already stripped out by biopython when you use the .id. this is one of the more annoying behaviors of biopython IMO, but in this instance it is useful. I've confirmed that my example code does this. It admittedly does not handle the underscore delimitation, which I missed. But that can be done with split much more simply. $\endgroup$ Nov 23 '20 at 22:45
3
$\begingroup$

You wrote,

mylist = re.match(H149xcV\_\w+\_\w+\_\w+)

try...

myregex = re.compile(r'H149xcV\_\w+\_\w+\_\w+')
mylist = myregex.match(myfile.id)

or

mylist = re.match(r'H149xcV\_\w+\_\w+\_\w+',myfile.id)

Watch your indentation it should be 4 spaces. Note you are losing the last '_\w+', but don't think it matters.

I think it is better to use,

mylist, _, _ = re.split(r' ', myfile.id)

The problem with your code is list will just return true ... and that is not what you want. If you are using split carefully watch the number of ' ' in your code you can also use r'\s'

You can also set maxsplit as a parameter, which might be useful.


myfile not defined.

Please try,

fasta_file = (sys.argv[0])

Generally your argv command isn't cool and not just [0] versus [1]. I forget the argv module everyone uses .. there's two of them, they are a bit fiddle but are very cool of options. You way in is okay, okayish, sort of okay ... Anyway you appear to be requesting the second input into the commandline, wherea you are probably placing the fasta file in the first position.

Also

print (fasta_file)

Check something happens that you picked up the file.

The other way in is,

mypath = '/User/username/location/fasta.fa'
for myfile in SeqIO.parse(mypath, "fasta"):

where mypath is a OS X style, you're probably using Linux so its... /home/username/pathtofile/fastafile.fa (something like that) or

for myfile in SeqIO.parse('/User/username/location/fasta.fa', "fasta"):

I forgot .... don't ever use list as a variable, by mylist is fine. Again I've not debugged my code.

$\endgroup$
4
  • $\begingroup$ Hi @Michael I tried doing that, but received a NameError indicating that the name 'myfile' is not defined. I appreciate your help! $\endgroup$ Nov 23 '20 at 19:24
  • $\begingroup$ @andnowmywatchbegins can you update your answer with this information? $\endgroup$ Nov 23 '20 at 20:11
  • 2
    $\begingroup$ list is used in Python for creating lists. I would recommend naming this variable something else. $\endgroup$ Nov 23 '20 at 20:16
  • $\begingroup$ list .... Yes!! Use mylist as a variable. instead $\endgroup$
    – M__
    Nov 23 '20 at 20:34

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.