compare and trim fasta files by header ID

Question

[joie.cannon]

Hi there, first time poster here!

I'm trying to compare two series of multiple sequence alignments in fasta format by header ID. I've dug around the archives, and I haven't found anything quite like what I need to do. File1.fasta contains 65 sequences, and File2.fasta contains 188 sequences. I need to trim File2.fasta so that it only contains the 65 sequences that match the headers in File1.fasta. The headers will match, but the sequences themselves do not.

For example:

If File1.fasta looks like this:

>Species1|sequenceID1|1
GIYFEANGHGTILNETVGDAIADLLGT---WLNLYTDLPQRQLKVAVKDRTM---IQT-TDAERRC-
>Species2|m.5506|0
GVYAEANGHATFLNSKTGDAICDLFLI---WQATYTPLSYRQKKLTIPDRYA---VVT-TDADRRV-

And File2.fasta looks like this:

>Species1|sequenceID1|1
MSDIAKE-----KLA--ALMLEFPKPA-----------GVILGYGTAGFRSRAD----ILPWIMIRIGILASLRSKVKQ-A-----
>Species1|sequenceID2|0
MSDIAKE-----KLA--ALMLEFPKPA-----------GVILGYGTAGFRSRAD----ILPWIMIRIGILASLRSKVKQ-A-----
>Species1|sequenceID3|0
MTHIVKE-----KLM--TLLLEYPKPE-----------GIIMGYGTAGFRARAD----TLPWIMIRIGLLAALRSKVKR-A-----
>Species2|m.6218|1
--------------------------------------GVEFHYGTSGFRMHAD----RLDTVIFRMGLLSALRSQAIGGK-----
>Species2|m.5506|0
-----------------------------------------------------------------DIHFHYGTAGFRSPAK----

I want the output to be this:

>Species1|sequenceID1|1
MSDIAKE-----KLA--ALMLEFPKPA-----------GVILGYGTAGFRSRAD----ILPWIMIRIGILASLRSKVKQ-A-----
>Species2|m.5506|0
-----------------------------------------------------------------DIHFHYGTAGFRSPAK-----

It looks as though SeqIO in BioPerl or Biopython could do it tidily, but I'm not quite there with my programming skills.

Any help would be deeply appreciated!

[Ashutosh Pandey]

1. Grep grep > all the header lines from file 1 into a new file.
2. Use grep -A1 -F -f. Google these parameters to know what they do and how they are used.

[Matt Shirley]

Using grep is a fine solution, but I want to point out that grep -A1 will only return the first line after the matched header line. I'm assuming the OP has some multiline sequences in her FASTA file.

[Ashutosh Pandey]

You are right. Mine won't work if it is a multiline fasta file.

[joie.cannon]

Thanks! I was overcomplicating it and knew there should be a straightforward solution. I have already run the files through a program to remove newlines, so this should be ideal.

[Ashutosh Pandey]

What you mean by "remove newlines"? The fasta sequence should be contained within a single line for "-A 1" to work. That's it.

[Matt Shirley]

I think "remove newlines" was just in the context to the sequence portion of the file - to make multiline sequences into single lines.

Answer 1

[Matt Shirley]

$ pip install --user pyfaidx
$ faidx -i File1.fasta | cut -f1 | xargs faidx File2.fasta > File3.fasta

This command will grab a list of the sequence names in File1, one per line, then pass those back to faidx as regions to grab from File2, and send the results to File3.

pyfaidx also has an object-oriented interface you can use within a Python script for when you outgrow the "faidx" utility.

Answer 2

[venu]

Hello, go like this

1. grep '^>' file1.fasta > file1_headers.txt
2. Download 'faSomeRecords' from the link http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ and change its mode by chmod +x faSomeRecords
3. ./faSomeRecords file2.fasta file1_headers.txt newFile2.fasta

Now your newFile2.fasta contains 65 sequences with the id from file1.fasta and sequence information of file2.fasta

Answer 3

[5heikki]

This only works IF 1) there are no linebreaks in the seqs and 2) there are no tabs in the headers. The tab in -t " " is a literal one (ctrl+v+tab).

join -1 1 -2 1 -t "       " -o 2.1,2.2 <(cat File1.fasta | tr "\n" "\t" | tr ">" "\n" | sort -k1,1) <(cat File2.fasta | tr "\n" "\t" | tr ">" "\n" | sort -k1,1) | awk 'BEGIN {FS="\t"} {print ">"$1"\n"$2}' | tail -n +3

>Species1|sequenceID1|1
MSDIAKE-----KLA--ALMLEFPKPA-----------GVILGYGTAGFRSRAD----ILPWIMIRIGILASLRSKVKQ-A-----
>Species2|m.5506|0
-----------------------------------------------------------------DIHFHYGTAGFRSPAK----