0
$\begingroup$

I want to remove sequence VRE32514 – it doesn’t belong and thus is the reason it lacks additional metadata. However I tried implementing this code from a similar question:

awk 'BEGIN{while((getline<">VRE32514")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' *.fasta | grep -v ">VRE32514" *.fasta

With no success and the proposed solution in R using seqinr doesn't seem to work either for myself.

#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("test.fasta")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), "|"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c(">VRE32514")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fasta")

Below is an example of what my Multi-FASTA appears as and the desired output:

>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT
>VRE32514
AATCAGATTGTTCGCCATGAGCTGCTGAGAATGATATCATGGAAGGGCGGCATCGAAACA
GGCTTTAAATTAAGTGTAGGCAAGAACTATAAGTTTATTGAAAGGTATATATCCGAGGAT

Desired Output:

>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT
Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

You can use samtools faidx to first index your file and then retrieve only those sequences you specify. For this grep all the sequence headers you want and then remove the unwanted ones, followed by retrieval of the fasta entries with faidx:

#/ index:
samtools faidx your.fa

#/ run:
grep '^>' test.fa \             # grep names
| grep -v '>VRE32514' \         # remove unwanted
| tr -d '>' \                   # remove ">" from name
| xargs samtools faidx your.fa  # get the entries you want

Output:

>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

In R one could use Biostrings:

library(Biostrings)

fasta <- readDNAStringSet("path/to/your.fa")
fasta_new <- fasta[!names(fasta) %in% c("VRE32514")]

#/ save to disk:
writeXStringSet(fasta_new, "path/to/outputfile.fa")
Improve this answer
$\endgroup$
2
  • $\begingroup$ Your solution in R using Biostrings works perfectly, thank you! However I couldn't understand how to get your UNIX` solution to work. When I enter grep '^>' example.fasta | grep -v '>VRE32514' | tr -d '>' | xargs samtools faidx example.fasta I get the following error message: \n not found in FASTA file, returning empty sequenceood|2016-12-07 >VRE32491|736|PUH-10C|Blood|2016-12-07 [faidx] Failed to fetch sequence in VRE32491|736|PUH-10C|Blood|2016-12-07 $\endgroup$
    – MicroBiostat
    Jul 14, 2021 at 1:13
  • $\begingroup$ @MicroBiostat It works fine on my machine. Maybe something odd wirh your fasta, e.g. empty lines somewhere? $\endgroup$
    – user3051
    Jul 14, 2021 at 12:58
0
$\begingroup$

Here's a solution using Python:

from Bio import SeqIO

sequences = SeqIO.parse('test.fasta', 'fasta')
filtered = [seq for seq in sequences if 'VRE32514' not in seq.id]

with open('output.fasta', 'wt') as output:
    SeqIO.write(filtered, output, 'fasta')
```
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.