How to remove duplicates from fasta file but keep at least one per group based on header

user_13995211 0
Viewed 492

I have a multifasta file that looks like this:

( all sequences are >100bp, more than one line, and same lenght )

>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage3_samplenameD
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA

I need to remove the duplicates BUT keep at least on sequence per lineage. So in this simple example (Notice samplenameA,C and D are identical) above I would want to remove only samplenameD or samplenameC but not both of them. In the end I want to get the same header information as in the original file.

Example output:

>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA

I found out a way that works to remove just the duplicates. Thanks to Pierre Lindenbaum.

sed -e '/^>/s/$/@/' -e 's/^>/#/'
file.fasta  |\
tr -d '\n' | tr "#" "\n" | tr "@"
"\t" |\
sort -u -t '  ' -f -k 2,2  |\
sed -e 's/^/>/' -e 's/\t/\n/'

Running this on my example above would result in:

>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG

—> so losing the lineage 3 sequence

Now I’m just looking for a quick solution to remove duplicates but keep at least one sequence per lineage based on the fasta header.

I’m new to scripting... any ideas in bash/python/R are welcome.

Thanks!!!

fasta python
1 Answers

In this case I can see two relatively good alternatives. A) look into existing tools (such as the Biopython library, or FASTX toolkit. I think both of them have good commands to do most of the work here, so it may be worthwhile to learn them. Or, B) write your own. In this case you may want to try (I'll stick to python):

loop over the file, line-by-line, and add the lineage/sequence data to a dictionary. I suggest using the sequence as a key. This way, you can easily know if you already encountered this key.

myfasta = {}
if myfasta[sequence]:
    myfasta[sequence].append(lineage_id)
else:
    myfasta[sequence] = [lineage_id]

This way your key (sequence) will hold the list of lineage_ids that have the same sequence. Note that the annoying bits of this solution will be to loop over the file, separate lineage-id from sequence, account for sequences that may extend to multiple lines, etc.

After that, you can loop over the dictionary, and write the sequences to file by using only the first lineage_id from the list within the dictionary.

user_43587850
Related