This is a test version of Biostars. For the public version, visit https://www.biostars.org.
remove sequences with non-canonical nucleotides from fasta file

I want to print sequences form fasta file which do not have non-canonical nucleotides. Example fasta is:

>1
ATAcctcatctaGTGTG
ATGCTGCTAGTZ
>2
agagagagagagagag

My code is 

from Bio import SeqIO
for record in SeqIO.parse("test.fasta", "fasta") :
    if set(record.seq) <= "ATCGatcg":
                print record

Instead of print the sequence of >2, it prints both.

What am I doing wrong? Thanks

seqio fasta
• 3,395 views • link
updated by Eric Lim Eric Lim • written by grant.hovhannisyan grant.hovhannisyan

2 answers

linearize and filter with awk:

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' input.fa |\
awk -F '\t' '($2 ~ /^[ATGCatgc]+$/)' |\
tr "\t" "\n"

using bioalcidaejdk:

$ java -jar dist/bioalcidaejdk.jar -e 'stream().filter(F->java.util.regex.Pattern.matches("^[ATGCatgc]+$",F)).forEach(S->println(">"+S.getName()+"\n"+S));' input.fa
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you Pierre, I accept your answer. But anyways, what is wrong with SeqIO?

• 0 views • link
written by grant.hovhannisyan grant.hovhannisyan 
from Bio import SeqIO
from Bio.Alphabet.IUPAC import IUPACUnambiguousDNA
for record in SeqIO.parse("test.fasta", "fasta"):
    if set(record.seq.upper()) <= set(IUPACUnambiguousDNA.letters):
       print(record)
• 0 views • link
written by Eric Lim Eric Lim

Log in to answer this question.

More posts like this