# How to remove duplicates from a fasta file using python

I am using the following command:

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(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
print (">"+list.group(1))


And this is providing me with duplicates of the same gene:

>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2
>H149xcV_Fge342_r3_h2
>H149xcV_Fge342_r3_h2


How can I reformat my command so that I only receive unique gene id's:

>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

• Why do you need to use python for this? Nov 24 '20 at 0:21
• I am using nano within commandline to write my python command. @RamRS Nov 24 '20 at 0:23
• I'm asking why, not how. There are command line tools that de-dup fasta, so why reinvent the wheel? Nov 24 '20 at 2:28
• @RamRS The only command-line tools that I am aware of are sort|uniq but I am not sure as to how I would use it? I must apologize, I am kind of new to bioinformatics, so I am still getting the hang of it. I know I can output my python into another file, but is there a way to obtain the unique id's afterwards? Thanks! Nov 24 '20 at 2:46
• You can use seqkit rmdup instead reinvent the wheel. bioinf.shenwei.me/seqkit/usage/#rmdup Nov 24 '20 at 9:47


from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
if len(myfile) > 250:
gene_id = myfile.id
if gene.id not in unique:
unique.append(gene.id)
mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
print (">"+mylist.group(1))

$$$$

• Yeah, I would do something along these lines too. Though, I would use a set or dict insead of list if the fasta file is big... Jan 27 at 15:55

Your question isn't well formatted, for example list is used as a variable name where it should be reserved as a keyword...

print (">"+list.group(1))


However, the answer you're after is going to be something like:

print (">"+set(name_of_list.group(1)))


Anytime you want to get the unique elements of a list, just wrap it in set()

Perhaps itertools might be an option? From the itertools page / moreitertools project https://docs.python.org/3/library/itertools.html#itertools-recipes

def unique_everseen(iterable, key=None):
"List unique elements, preserving order. Remember all elements ever seen."
# unique_everseen('AAAABBBCCDAABBB') --> A B C D
# unique_everseen('ABBCcAD', str.lower) --> A B C D
seen = set()
if key is None:
for element in filterfalse(seen.__contains__, iterable):
yield element
else:
for element in iterable:
k = key(element)
if k not in seen:
yield element


I would use a dict so that the lookup time is constant. If you are creating an array with all the unique gene IDs, then checking if a gene is in that array, your runtime is O(N2). When you create an entry in the dict, the time to access is constant, and if the entry is already present, it will just overwrite it. Something like this should do it for you:

from Bio import SeqIO

genes = {}
fasta_file = sys.argv[1]
for entry in SeqIO.parse( fasta_file, "fasta" ):
genes[ entry.id ] = entry.seq



Then you can write a new, unique fasta file:

with open('unique_genes.fasta', 'w') as out:
for key,value in genes.items():
out.write( '>%s\n%s\n'%(key,value) )


I use the following code. It doesn't only remove repeated sequences but similar sequences (with a quite simplistic approach). I use it just sporadically and I guess it can be largely depurated. Remove the debug prints if these annoy you.

import sys
import requests
import shlex, subprocess
import random
from Bio import SeqIO
import editdistance

if len(sys.argv)<3:
print("Uso: "+sys.argv[0]+" input_fasta output_fasta")
sys.exit(0)

savefile = sys.argv[2]
content_file = sys.argv[1]

class Sequence:
id=""
sequence=""

def are_similar(seq1, seq2):
l1 = len(seq1.sequence)
l2 = len(seq2.sequence)
if abs(l1-l2) < l1/40:
d = editdistance.eval(seq1.sequence[0:250], seq2.sequence[0:250])
#d = editdistance.eval(seq1.sequence, seq2.sequence)
print(seq1.id+" "+seq2.id+" -> "+str(d))
#if d < l1/20:
if d < 50:
return True
return False

def find_repeats(sequences, seq, repeated):
b=False
for seq2 in sequences:
if(b and not seq2.id in repeated):
if are_similar(seq, seq2):
if seq2.id==seq.id:
b=True

sequences=[]

res=""

for record in SeqIO.parse(content_file, "fasta"):
seq=Sequence()
seq.id=record.id
seq.sequence=record.seq
sequences.append(seq)

repeated_sequences=set()
for seq in sequences:
find_repeats(sequences, seq, repeated_sequences)

print(repeated_sequences)

def mySortFunc(e):
return len(e.sequence)

sequences.sort(key=mySortFunc)

for seq in sequences:
b= seq.id in repeated_sequences
if not b:
res+=">"+seq.id+"\n"
res+=str(seq.sequence)+"\n"

output_file = open(savefile, "w")

output_file.write(res)

output_file.close()
`

EDIT: I see that's not exactly what OP is asking ...