This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Extract Sequences/Accession Numbers That Are Shared By A Number Of Fasta Files

Hi everybody,

I have 8 fasta files each containing 100 sequences and I want to extract those sequences that are present in all 8 files, or thus eliminating those sequences that are only present in a subset of the 8 files. Sequences are identified by their genbank accession number, so I'm guessing it should be possible by extracting the accession numbers that are shared.

I was wondering whether there is an existing Perl script to do this?

Kind regards,

Sam

perl fasta script genbank
• 5,218 views • link
updated by cts cts • written by samlambrechts299 samlambrechts299

It may help if you could post a small sample of your input and what you expect as output.

• 0 views • link
written by SES SES

Are you confident that the sequences sharing the same identifiers are the same sequence? If you aren't, you may need to calculate a checksum for each sequence (using something like md5) to be sure...

• 0 views • link
written by sarahhunter sarahhunter

1 answer

you can determine the set of sequences shared using unix tools:

grep -oP '(?<=\>)\S+' *fa | sort | uniq -c | awk '{if($1 == 8)print $2}'

you could then extract these sequences from anyone of the files using a variety of scripts. I use one called contig_extractor.pl that in its simplest form takes a list of identifiers and a fasta/fastq database and returns a subset of those sequences

• 46 views • link
written by cts cts

Log in to answer this question.

More posts like this