I've found several threads on this (rather simple) topic but none quite simple enough, which is to remove entries in a fasta file based on their one liner >name, which in my case is numeric (gi).
Based on Pierre Lindenbaum's posting on other comments, you would linearise the sequences and then sort by column 1 (as opposed to column 2 if you wanted to sort by sequence). And then you'd employ sort unique and sed?
>123456
AAAGTGTGTAGGAAGATGTGATGCCTCGAGATGC
>123456
AAAGTGTGTAGGAAGATGTGATGCCTCGAGATGC
There are no spaces between characters or lines in my file.
•
written
by
oliver.bayfield •
8 answers
Explanation:
/^>/: when matches line start with>do this block{f=!d[$1];d[$1]=1}.f=!d[$1]: f is false only if sequence name does exists ind[$1], and keeps being false until new sequence.d[$1]=1: register sequence name;fis true, print line.
•
written
by
biocyberman
Try this command:
awk 'BEGIN{RS=">"}NR>1{sub("\n","\t"); gsub("\n",""); print RS$0}' file.fa | awk '!seen[$1]++'
In the first part the fasta is converted to tabular format and in the second part the duplicated ID's are removed.
You'll need to transform to fasta format again, but it's not complicated ;)
airan thanks that seems to work on a test file - any chance you could guide me doing the conversion back to fasta? I'm rather inexperienced here. Thanks
•
written
by
oliver.bayfield •
Extending airan solution for converting it to a fasta format:
awk 'BEGIN{RS=">"}NR>1{sub("\n","\t"); gsub("\n",""); print RS$0}' biostarhelp.txt | awk '!seen[$1]++' | awk -v OFS="\n" '{print $1,$2}'
And in case your headers have spaces and you also want to consider the sequence
awk 'BEGIN{RS=">"}NR>1{sub("\n","\t"); gsub("\n",""); print RS$0}' biostarhelp.txt | awk '!seen[$0]++' | awk -v OFS="\n" '{for(i=2;i<NF;i++) head = head " " $i; print $1 " " head,$NF; head = ""}'
Thanks for your comments. I went with the orignal suggestion of Pier's:
sed -e '/^>/s/$/@/' -e 's/^>/#/' filein.fa | tr -d '\n' | tr "#" "\n" | tr "@" "\t" | sort -u -k1,1 | sed -e 's/^/>/' -e 's/\t/\n/' > fileout.faThe first line of the output contains a '>' only, which I deleted manually.
Please see: Best bioinfo one-liners? for a better way to linearize fasta. At the end you can convert the lines to fasta using tr "\t" "\n"
•
written
by
Pierre Lindenbaum
This would work assuming you have sequence in one line and it is not split over multiple lines.
#!/usr/bin/perl-w
use strict;
use warnings;
my %id2seq=();
my $key = '';
while(<>){
chomp;
if($_ =~ /^>(.+)/){
$key = $1;
}else{
$id2seq{$key} = $_;
}
}
foreach(keys %id2seq){
print join("\n",">".$_,$id2seq{$_}),"\n";
}
@ biolab
It doesn't work if 2 same id's are adjacent some thing like this
>123456
AAAGTGTGTAGGAAGATGTGATGCCTCGAGATGC
CCGG
>123456
AAAGTGTGTAGGAAGATGTGATGCCTCGAGATGC
CCGG
IT will append it for the same header
•
written
by
Varun Gupta
Hi, Varun Gupta, your point is very good. I further modified the script to allow multi-line sequences. If a header is duplicated, only the 1st is outputted.
#!/usr/bin/perl
use strict;
use warnings;
my (%id2seq, %seen);
my ($key, $duplicate);
while(<>) {
chomp;
if($_ =~ /^>(.+)/){
$key = $1;
if (exists $seen{$key}) {
print STDERR "Attention: header $key duplicated.\n";
$duplicate = 1;
} else {
$seen{$key} = 1;
$duplicate = 0;
}
} else {
($duplicate == 1) ? (next) : ($id2seq{$key} .= $_);
}
}
foreach(keys %id2seq) {
print join("\n",">".$_,$id2seq{$_}),"\n";
}
Assuming dup.fasta has one >defline followed by one line of sequence
awk '/^>/{id=$0;getline;arr[id]=$0}END{for(id in arr)printf("%s\n%s\n",id,arr[id])}' dup.fasta > uniq.fasta
There is no requirement duplicates be adjacent, there is no guarantee the output order is related to the input order.
Here is my free program on Github Sequence database curator
It is a very fast program and it can deal with:
- Nucleotide sequences
- Protein sequences
It can work under Operating systems:
- Windows
- Mac
- Linux
It also works for:
- Fasta format
- Fastq format
Best Regards
•
written
by
Eslam Samir
If you use perl, you could put the sequences in a hash with the name/number as the key and the sequence as the value. Then you could print it to a new file, or empty your file and print it back in. Every duplicate you come across will just replace the previous one.
Could you give an example @anp375?
use Bio::SeqIO;
use strict;
my $protein_fasta = "/Somefilepath/protein.txt";
my $protein_out = ">/Somefilepath/reduced.txt";
my $seq_in = Bio::SeqIO->new(-file => "$protein_fasta", -format =>'Fasta');# This is a constructor for SeqIO
# The filename specifies whether it is read or write
# The big arrows are called fat commas. They are just commas used for documenting. -format => 'Fasta' is
# just -format, 'Fasta' but the arrow shows the relationship for readability
my $seq_out = Bio::SeqIO->new(-file => "$protein_out", -format => 'Fasta');# SeqIO out because of >
# Those constructors set up the filehandles. The method to get the sequences inside is called 'next_seq'
# It returns a generically formatted sequence rather than Fasta
my %seqs; # This is a hash. I'm assuming the name you are using is an accession number, so I'll make the keys accession numbers. Only one of each will be left.
while(my$seq = $seq_in -> next_seq){
$seqs{$seq->accession_number} = $seq;
}
foreach(values %seqs){
$seq_out->write_seq($_);# write_seq converts the generic sequence to a Fasta sequence and writes
# it to the file
}
I hope that works. If the number is not the accesion number, there is a list of other ids here: http://search.cpan.org/dist/BioPerl/Bio/Seq.pm#accession_number
Copied from my other post: How to remove duplicate sequences in fasta file using python?
Learn to use Biopython library. It's handy as hell. You can use any format as in/out
from Bio import SeqIO
with open('output.fasta', 'a') as outFile:
record_ids = list()\
for record in SeqIO.parse('input.fasta', 'fasta'):
if record.id not in record_ids:
record_ids.append( record.id )
SeqIO.write(record, outFile, 'fasta')
Log in to answer this question.
-
error in Genome Mepping by BWA tools in Linux
written by drkuldeepbiotech •
-
Prioritizing/rank ordering genes from two lists: averaging values?
written by kelen
-
single cell ATAC-seq cicero gene activities in Seurat
written by el24 •
-
Remove duplicates from fasta keeping at least one sequence per group based on header
written by vixelaa •
-
unix, awk or sed question on splitting lines
written by geneart$$ •
-
Easy way to filter gene list based on chromosome of origin?
written by steve
-
Simple clustering of start/stop nt positions in a text file column?
written by oliver.bayfield •
-
Mapping microRNA reads to Genome and non-coding RNAs
written by Zammy •
-
Should We Be Using Much More Json In Our Delimited Data Formats?
written by Casbon •
-
Taxonomy Of Blast Hits
written by Darked89
linerarize,
sortusing options-k1,1 -u, move back to fasta usingtrIs this correct?
Update, the above (from Pierre Lindenbaum) does the job. Very good.
The only thing, I dropped the
-t ' 'and-fflags in sort (didn't seem necessary?). And the first line in the output file gives a single line of>, which I manually deleted.