This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Obtain sequences that did not appear in your fasta file

Good day, I am trying to identify which sequences did not appear from my list. I have file1.txt which is a list of genes

#gene1
#gene2
#gene3
#gene4
#gene5

I also have file2.fa

#>gene1
ACTAGA
#>gene3
ACATGA
#>gene6
AGATA

I want to be able to identify the genes that are not found in file2.fa based on file1.txt list sample output would be

#gene2
#gene4
#gene5

I tried for i in $(cat file1.txt); do perl -ne '/$i/ && print' file2.fa > output.txt; done it gives everything that appeared in the list. I tried diffirent iterations to get whats not on the list but I wasn't able to. Hope someone could help me with this. Thanks!

genome ngs perl unix bash
• 2,632 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by sicat.paolo20 • sicat.paolo20

I assume those # are not in your real data since that will break the fasta format.

• 1 views • link
written by GenoMax GenoMax

yup sorry. It was showing as something else in my laptop before I added the #

• 1 views • link
written by sicat.paolo20 • sicat.paolo20

Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time. You would not need to add # in that case.
code_formatting

Thank you!

• 1 views • link
written by GenoMax GenoMax

Thank you very much genomax. Will do that next time.

• 1 views • link
written by sicat.paolo20 • sicat.paolo20

2 answers 

comm -23 <(sort file1.txt ) <(grep '>' file2.fa | cut -c 2- | sort)
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

thank you very much. follow up question. What if the header in the fasta file has >gene3 <additional information="">. they were skipped. What would be the easiest way to include them?

• 1 views • link
written by sicat.paolo20 • sicat.paolo20

What would be the easiest way to include them?

insert a cut command before sort

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks for the help! PS. Big fan of your one liner scripts btw.

• 1 views • link
written by sicat.paolo20 • sicat.paolo20

If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted.

Upvote|Bookmark|Accept

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

output:

$ grep \> file2.fa | sed 's/>//g' | grep -vf - file1.txt 
gene2
gene4
gene5

input:

 $ cat file2.fa 
>gene1
ACTAGA
>gene3
ACATGA
>gene6
AGATA

$ cat file1.txt 
gene1
gene2
gene3
gene4
gene5
• 1 views • link
written by cpad0112 cpad0112

using comm,sed and grep:

$ grep \> file2.fa| sed 's/>//g'  | sort | comm -13 - <(sort file1.txt)

gene2
gene4
gene5
• 1 views • link
written by cpad0112 cpad0112

Thanks for the help. This one worked. The first one had an error message about invalid range end. Thanks!!!

• 1 views • link
written by sicat.paolo20 • sicat.paolo20

Log in to answer this question.

More posts like this