Pulling out gene names from a FASTA file with a list of IDs

Question

[jvire1]

Does anyone know of a simple scripting solution to take a list of accessions and pull out the gene name from the header of FASTA sequences?

For instance given the accession: XP_016469325.1

and given the FASTA entry:

XP_016469325.1 Nicotiana tabacum|C3H|C3H family protein MEEELLKRNTDCVYFLASPLTCKKGIECEYRHSEIARLNPRDCWYWLAESCLNPTCAFRH PPLESHAETSSESAPPQHKSAVPVNKTNVPCYFYFNGYCIKGERCSFLHGPDDGTTTWKS SKIASGVPDGPTAEKKTSVGSETGPASVEKPSNSSETGSKAAAHEYIKSQVDLISMTNDV GEQSASHETSGSPSEEATAVRLDSLVPAEGFTQGGSDLSPDWSSDEEVEDNVEREEWLES SPGFDVLVDDRIEGWSHKDDHSYLLQHDRECDERFAGYDFENNLEYDPAYPDMRIVSDEE LDDSYYSKVESHEVNEYAREIVIPAHGRQSIPHKRKFPREPGFCARGNVDLRDLLKKRRV IESDPPNYLSRRLDLSRFNAREQCRDRHRPQGSRWMPQSLASKLESNSSFSSGFVDATRL EGANQLKKLRQSHRSSYRQQHFKDRRRGRSQPFANETPRRMASRQRSTEVPKIFGGPKTL AQIREEKIKGREDGNSFERTVPSGGSEREDFSGPKPLSEILKDKRRLSSVVNFSN

I would like to have output the gene name "C3H" that is spanned by the "|"

This script I modified from a previous post can grab the gene names, however I'm not sure how to only get the gene names corresponding to a separate list of accessions (accessions.list).

with open('PlantTFDB_ALL_TF_pep.fas','r') as f:
    for line in f:
        if '>' in line:
            line = line.strip().split('|')
            print(line[1])

Answer 1

[shenwei356]

File acc.txt contains the accession IDs:

$ cat acc.txt 
XP_016469325.1

Firstly, getting the mapping relationship between acc and gene name and saving them to acc2gene.tsv

$ grep '>' seqs.fa  | \
    perl -ne 'next unless />(\S+).+\|(.+?)\|/;   print "$1\t$2\n";'
XP_016469325.1  C3H

Secondly, joining acc.txt and acc2gene.tsv using csvtk (http://bioinf.shenwei.me/csvtk/download/) or other tools：

$ csvtk join -H -t -k acc.txt acc2gene.tsv | csvtk cut -H -t -f 2

---

Original answer

Searching header line using grep -f, and capturing accession ID and gene name using regular expression in Perl.

$ grep '>' seqs.fa  | \
    grep -f acc.txt | \
    perl -ne 'next unless />(\S+).+\|(.+?)\|/;   print "$1\t$2\n";'
XP_016469325.1  C3H

If you just want the gene names, that would be easier:

$ grep '>' seqs.fa | grep -f acc.txt | cut -d '|' -f 2
C3H

[jvire1]

Thank you for this awk approach, unfortunately I am needing each line in the accession to be searched one at a time against the FASTA database, retaining duplicate instances. The idea here was to take the accessions and get back the gene name for each. For example an acc.txt file that looked like this:

Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
KZV22762.1
KZV22762.1
KZV22762.1

I would expect to see this:

NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
CAMTA
CAMTA
CAMTA

The output of this awk command however is this:

G2-like
GeBP
bZIP
ERF
B3
B3

which looks like it is just extracting the gene name for each FASTA entry regardless of the acc.txt:

loki@rnoyes-Precision-T3610:~/Desktop/PlantTFDB/apomict$ head -10 PlantTFDB_ALL_TF_pep.fas 
>KFK36254.1 Arabis alpina|G2-like|G2-like family protein
MYSAIRSLPVDGSMGEYSDGTNLPIEACLVLTTDPKPRLRWTTELHDRFVDAVTQLGGPD
KATPKTIMRTMGVKGLTLYHLKSHLQKFRLGRLPCKESTDNSKDVSSVAESQDTGSSTTS
SLRLAAQEQKESYQVTEALRAQMEVQRRLHEQLEVQQRLQLRIEAQGKYLQTILEKACKA
IEEQAVAFAGLEAAREELSELAIKVSNGCQSTFDTTKMTIPSLSELAVAIEHKNNCSAES
SLTSSTVGSPVSAAMMKKRHRGVFGNGDSVVVGHEAGWVMPSSSIG
>KFK25038.1 Arabis alpina|GeBP|GeBP family protein
MAPKHTDTIQNPPMPSSSSEEESASSGEDSDSSKKHNSASDSDPIKTKIVVTKPESSGSM
KTTTKSSVVVAEPESTTAKRPLKETEEDVEVKKKKMKTELVMKNQEPEEKVSGDETKKQL
MFQRLFSENDEIALLQGILDFTSTRGDPYEKMDVFCEYVKKLINFDANKSQLVTKIRRLK
loki@rnoyes-Precision-T3610:~/Desktop/PlantTFDB/apomict$

Answer 2

[st.ph.n]

Make another file with the acc #, and check if the header for the accesion number.

with open('acc_list.txt', 'r') as f:
    acc = [line.strip() for line in f]

with open('PlantTFDB_ALL_TF_pep.fas', 'r') as f:
    for line in f:
        if line.startswith(">"):
            if line.strip().split('>')[1].split(' ')[0] in acc:
                print line.strip().split('|')[1]

[jvire1]

After I looked at the output more closely I realized that it isn't doing quite what I want it to.

The beginning of the list of identifiers looks like this:

Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
Zpz_sc04830.1.g00020.1.sm.mk
KZV22762.1
KZV22762.1
KZV22762.1

However the beginning of the results look like this:

WRKY
B3
STAT
MYB
Dof
bHLH
NAC
HD-ZIP

I'm not very familiar with scripting yet, but if I understand the logic of the code:

with open('acc_list.txt', 'r') as f:
    acc = [line.strip() for line in f]

with open('PlantTFDB_ALL_TF_pep.fas', 'r') as f:
    for line in f:
        if line.startswith(">"):
            if line.strip().split('>')[1].split(' ')[0] in acc:
                print line.strip().split('|')[1]

It is reading through each FASTA entry of the 'PlantTFDB_ALL_TF_pep.fas' and checking if it exists in the 'acc_list.txt' and then printing the parsed out gene name?

Is there a way to instead search each identifier from the list against the FASTA file and then print each corresponding gene name (including duplicate instances) like so:

NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
CAMTA
CAMTA
CAMTA

[Joe]

Can you provide more of your input file? You've asked that the field delimited by | symbols be returned, based on the accession it's paired with in the fasta header as far as we can work out - correct?

If you aren't getting the results you expect, it's likely you fasta headers aren't all in exactly the same format, and the parsing of the headers is breaking down.

[jvire1]

Thank you for all your help. Here is three FASTA entries in the PlantTFDB FASTA file, with the first one corresponding to the first accession ID in the acc_list.txt:

>Zpz_sc04830.1.g00020.1.sm.mk Zoysia pacifica|NF-YA|NF-YA family protein
MEAKFFRFLKLVGVGFKARSESQGRELFLKLGYSHEVLFTAPPAVRVFCFKPNIICCTGI
DKDRVHNFAGAVRSCKPPEVYKGKGILYFDEVIKLKPGKKQKKQMNNCPRLPHLEGVQYP
GSTAEVPEGSATATPHMSVDNSGPSELIAGQRRALRLRQLPSLVPGRATIFVSAKQYEAI
LRRRTQRAKAGISRPSIKGGANTDISQATSSVPAEDFYGPFVATNEDEEHVGEDAFDSIL
NLNSPDFTTLLRLMMGDKYNTDISEDVTWILM
>Zpz_sc03704.1.g00020.1.am.mk Zoysia pacifica|MYB|MYB family protein
MGRSPCCEKAHTNKGTWTKEEDQRLTMGRSPCCEKAHTNKGTWTKEEDQRLVAYIKAHGE
GCWRSLPKAAGLQRCGKSCRLRWINYLRPDLKRGNFTDEEDDLIIKLHQVLGNKWSMIAG
RLPGRTDNEIKNYWNTHIKRKLLARGIDPKTHRPLSEAATAAASSSTHLQDHPARPSSSP
ETSCCAGQSGDEGSASGSLPETSTQRQLASCIIDLNLSIGPPQHQPSSSPPPPPTRPEAQ
AAGSSTTATAGTSSSSNSEKERICLPRNIQPSRAIFSPSCLMPASFSHTREAQSLPSSSC
RDLPAATPPSPAALRPATVPLAAAPPSVLLPCRSSPRALTRTAPCALTHDCTSLRS
>Zpz_sc00930.1.g00160.1.am.mk Zoysia pacifica|WRKY|WRKY family protein
LADEVAKQGAAAVVKPVVLPLASRMPPSPFLPSPRLLSDSADGPSMAAATRTGGWRIGPR
TTDPAAGTADRRQQRWIRRLGRVTGSRDAGPEVGTAARRTARGGEGRGKAEVRWASLPRL
AVVTTVLRLVRDQLLACSSCGGTAASARPSTATGPRGRRSPAGEQRRGGDGDAGDRGWGC
GLSRATGNNTSLTRIMAQISVHHGTKAPRPIPNLNAWINGGSFILPDGRWVKIRLSSGCC
CCPCDAAVSACSGKGTSAATNQRGGAGGDASSRTPAKKLRRHQIGSEARCETRGARALIK
TATRRTVGCRMHVALVQARATAAACLSSHPHAMEYRHPYSSGSASASANYGDASPAADFR
PAAFGAGRPAVRRAPLDVFDYLSDEGAQPAPAAVPAVPTFRSPPPRSMPAVQVVPDAAGY
GHAARPRGFMYLNIFFPRGAASVAVAAEEPTTMTDRIAFRMRSDKEVLDDGYRWRKYGKK
AVKNSPNPRGAASVAVAAEEPTTMTDRIAFRMRSDKEVLDDGYRWRKYGKKAVKNSPNPR
NYYRCSTEGCNVKKRVERDRDDPSYVVTMYEGVHNHASPGTVYYATQDAASGRFFVSGTH
HSSGPLKL

I have manually looked through the FASTA, which contains 320,370 protein entries and all I have observed have headers in the same format.

The script works great at pulling out the gene name from the headers, I think the issue might be with the:

if line.strip().split('>')[1].split(' ')[0] in acc:

I'm not sure if I understand this right, but it seems like the script is:

1. going through each FASTA entry one by one, 
2. checking if the header matches an entry in the acc_list.txt, 
3. printing the gene name if the FASTA header contains any accession ID in acc_list.txt.

What I'm aiming for is for each line in acc_list.txt to be searched for the matching accession in the FASTA and print out the gene name. The parsing out of gene names works fine, I'm just unsure of the logical argument that could take each line (including duplicates) one at a time in the acc_list.txt and search it against the FASTA file.

Answer 3

[Alex Reynolds]

Via awk:

$ awk '($0~/^>/){ n=split($0,a,"|"); print a[2]; }' in.fa > out.txt

Answer 4

[Joe]

I realise you've already got an answer, but if you want a more 'fully fledged' alternative, this is my day-to-day general purpose fasta retrieval script:

You can look for a fasta by a keyword within it's header, and it'll pull out the sequences. You can pass in a list or a string as the search query. To get your sequences back with just the |CH3| bit, you can just tweak the lines at the bottom of the script to output them as you wish.

	# Extract fasta files by their descriptors stored in a separate file.
	# Requires biopython
	from Bio import SeqIO
	import sys
	import argparse
	def getKeys(args):
	"""Turns the input key file into a list. May be memory intensive."""
	with open(args.keys, "r") as kfh:
	keys = []
	for line in kfh:
	line = line.rstrip('\n')
	line = line.lstrip('>')
	keys.append(line)
	return keys
	def main():
	"""Takes a string or list of strings in a text file (one per line) and retreives them and their sequences from a provided multifasta."""
	# Parse arguments from the commandline:
	try:
	parser = argparse.ArgumentParser(description='Retrieve one or more fastas from a given multifasta.')
	parser.add_argument(
	'-f',
	'--fasta',
	action='store',
	required=True,
	help='The multifasta to search.')
	parser.add_argument(
	'-k',
	'--keys',
	action='store',
	required=True,
	help='A string provided directly, or a file of header strings to search the multifasta for. Must be exact. Must be one per line.')
	parser.add_argument(
	'-o',
	'--outfile',
	action='store',
	default=None,
	help='Output file to store the new fasta sequences in. Just prints to screen by default.')
	parser.add_argument(
	'-v',
	'--verbose',
	action='store_true',
	help='Set whether to print the key list out before the fasta sequences. Useful for debugging.')
	parser.add_argument(
	'-i',
	'--invert',
	action='store_true',
	help='Invert the search, and retrieve all sequences NOT specified in the keyfile.')
	args = parser.parse_args()
	except:
	print('An exception occured with argument parsing. Check your provided options.')
	sys.exit(1)
	# Main code:
	# Call getKeys() to create the list of keys from the provided file:
	try:
	keys = getKeys(args.keys)
	# If -k/--keys was provided with a string, not a file path, the IO error is used as the indicator
	# to switch to expecting a string only, rather than a file.
	except IOError:
	keys = args.keys
	else:
	print("Couldn't determine a key from your provided file or string. Double check your file, or ensure your string is quoted correctly.")
	if args.verbose is not False:
	if args.invert is False:
	print('Fetching the following keys from: ' + inFile)
	for key in keys:
	print(key)
	else:
	print('Ignoring the following keys, and retreiving everything else from: ' + inFile)
	for key in keys:
	print(key)
	# Parse in the multifasta and assign an iterable variable:
	seqIter = SeqIO.parse(inFile, 'fasta')
	# For each sequence in the multifasta, check if it's in the keys[] tuple. If so, print it out:
	for seq in seqIter:
	if args.invert is False:
	if seq.id in keys:
	print(seq.format("fasta"))
	if args.outfile is not None:
	SeqIO.write(seq, outFile, "fasta")
	else:
	# If the --invert/-i flag is given, print all fastas NOT listed in the keyfile
	if seq.id not in keys:
	print(seq.format("fasta"))
	if args.outfile is not None:
	SeqIO.write(seq, outFile, "fasta")
	if __name__ == "__main__":
	main()

view raw fastafetcher.py hosted with ❤ by GitHub

Answer 5

[shenwei356]

I see. firstly, getting the mapping relationship between acc and gene name and saving them to acc2gene.tsv

$ grep '>' seqs.fa  | \
    perl -ne 'next unless />(\S+).+\|(.+?)\|/;   print "$1\t$2\n";'
XP_016469325.1  C3H

secondly, joining acc.txt and acc2gene.tsv using csvtk (http://bioinf.shenwei.me/csvtk/download/) or other tools：

$ csvtk join -H -t -k acc.txt acc2gene.tsv | csvtk cut -H -t -f 2

[jvire1]

Thank you so much! This worked and I was impressed how fast it processed everything. I am still going to see if I can get the solution in python, because I will likely need to be able to do similar tasks and it would be nice to know how to do it myself.

Here is the output:

NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
NF-YA
CAMTA
CAMTA
CAMTA