Hi everybody,
Sorry if this question is not relevant to this forum, however, I don't know anywhere to ask it. I have a fasta file containing thousands of sequences that their headers are like below:
Gene.1019::c44525_g4_i4::g.1019::m.1019 Gene.1019::c44525_g4_i4::g.1019 ORF type:complete len:339 (+) c44525_g4_i4:48-1064(+)
and in the another text file, there are two columns with hundred rows like below: (In this text file, just the part, not total, of headers from the fasta file exist)
Gene.1019 negative
I want to have the text file with the below information:
Gene.1019::c44525_g4_i4 negative
I'm familiar with Linux, but such a task for me as a biologist isn't easy. Could you please advise me with your helpful commands or tools.
Thanks in advance
3 answers
Here is then the final answer :-D
#Change identifiers of a fasta file (argument 1) based on additional file (argument 2)
#wdecoster
import sys
from Bio import SeqIO
with open(sys.argv[2]) as identifiersf:
#Create a dictionary with key the first field of the file and value the second field, remove empty lines (dictionary comprehension)
iddict = { line.split('\t')[0] : line.strip().split('\t')[1] for line in identifiersf.readlines() if not line == "" }
for seq_record in SeqIO.parse(sys.argv[1], "fasta"):
id2 = '::'.join(seq_record.id.split('::')[:2]) #This is the part we want to keep
seq_record.description = "" #The description of the fasta identifier is removed
try: #use try-except because not all keys are necessarily present in the dictionary
id1 = iddict[seq_record.id.split('::')[0]] #The key is the first part by splitting on '::'
seq_record.id = id2 + ' ' + id1 #Add the two strings, separated by space
except KeyError: #when not present, only the identifier is used
seq_record.id = id2
print(seq_record.format("fasta").strip()) #Avoid that newlines are printed by using .strip() on format("fasta")
Late to the party..
join -1 1 -2 1 -t $'\t' \
<(awk 'BEGIN{FS="::";OFS="\t"}{if(/^>/){print substr($1,2),$2}}' file.fasta | sort -k1,1 -t $'\t') \
<(sort -k1,1 -t $'\t' file.txt) \
| awk 'BEGIN{OFS=FS="\t"}{print $1"::"$2,$3}' \
> output
For big files one should add to the sort commands (if supported and enough cores): use up to 90% of RAM and 8 threads. This is for GNU sort. I'm not sure if BSD sort (e.g. in OS X) supports this. man sort to check.
-S 90% --parallel=8
awk learning curve isn't really that steep. I highly recommend it to anyone who has to deal with large text files..
Wow, looks fancy. Perhaps your solution is faster but I'd prefer python every day of the week over this :D then I don't need to be awake will coding.
Thanks, but it's not that fancy. I wrote it as a one-liner while thinking through the problem. @genomax2, maybe you M-w C-y'd only to the pipe on the second line leaving the process substitution open?
This should do the trick:
Code removed, see other answer. Answer was incomplete.
Save e.g. as connectingidentifiers.py, and execute as python input.fa identifiersfile.txt > results.fa
Note:
- Requires Biopython to be installed (pip install biopython)
- Probably contains python (>=2.7) synthax (which can be changed if really can't be upgraded on your system)
Let me know if it doesn't work as it should!
Edit: this produces an error (explicitly) when the identifier is not present in the dictionary. But this wasn't the desired performance.
Thank you, WouterDeCoster. I run it, but it gave me the error of "Aborted: Could not find Gene.4 in file of identifiers!". As I mentioned in my post, part of (not total) the headers of fasta file exist in the text file and it caused the error. Sorry, I tried to modify the code to solve the error, but I couldn't.
Right, didn't catch that. I just wanted to make sure that unexpected results were gained by explicitly mentioning that error in the script. But okay :p I'll make sure it doesn't produce the error then and just writes the identifier. In addition, I've added excessive comments to my code to explain the lines.
code see other answer
Will do as soon as confirmed.
Right, my bad. Didn't test my update. I made a rookie mistake here ;)
#Change identifiers of a fasta file (argument 1) based on additional file (argument 2)
#wdecoster
import sys
from Bio import SeqIO
with open(sys.argv[2]) as identifiersf:
#Create a dictionary with key the first field of the file and value the second field, remove empty lines (dictionary comprehension)
iddict = { line.split('\t')[0] : line.strip().split('\t')[1] for line in identifiersf.readlines() if not line == "" }
for seq_record in SeqIO.parse(sys.argv[1], "fasta"):
id2 = '::'.join(seq_record.id.split('::')[:2]) #This is the part we want to keep
seq_record.description = "" #The description of the fasta identifier is removed
try: #use try-except because not all keys are necessarily present in the dictionary
id1 = iddict[seq_record.id.split('::')[0]] #The key is the first part by splitting on '::'
seq_record.id = id2 + ' ' + id1 #Add the two strings, separated by space
except KeyError: #when not present, only the identifier is used
seq_record.id = id2
print(seq_record.format("fasta").strip()) #Avoid that newlines are printed by using .strip() on format("fasta")
Okay, I've submitted this as a final answer to this thread. Good luck with the rest of your analysis and my apologies for the mistakes which kept you waiting!
Log in to answer this question.
-
Help with Biopython for Beginner
written by cput •
-
Using transdecoder to filter a fasta file with preserving whole transcripts
written by poecile.pal •
-
Extracting fasta sequence from CDS matching with sequence ID
written by metallicasster •
-
sed command to extract sequence
written by mxlsherry1992 •
-
extract headers of a file using a second file containing a list of IDs
written by paraskevopou •
-
Extracting fasta sequences with identical header
written by seta
-
script for collecting selected IDs prediction result from Transdecoder .pep file
written by Farbod
-
Re-naming just part of header in fasta file
written by seta
-
Fasta sequence replacement based on header name
written by seta
-
FASTA sequence extraction using list of sequence name that are not identical with the headers?
written by seta
You can use any scripting language like python here. Step1 : Open fasta file Step2 : read fasta sequence and then store Gene.XXX as a key and fasta sequence as a value of that gene in a dictionary. Step3 : Open 2nd file and read it. Step4 : Open the dictionary and see any key of that fasta sequence match with any row of 1st column, if yes then grab 2nd column and print what you want. Even you can do without making dictionary by just reading two files same time.
I would suggest using a dictionary, definitely most straightforward.
While pseduocode is great @seta is looking for a solution that can be used right away :-)
Do you have any experience with (Bio)Python? Did I understand correctly that you want to keep the part of the identifier up to the second "::" and the part before the first "::" needs to match the field in the second file?
Sorry, I'm learning to work in this field and have no previous experience. In the example, I just have
Gene.1019::c44525_g4_i4as the identifier. Your straightforward help would be highly appreciated.Could you give me a few more examples? I'll write the script tomorrow morning (CEST).
Sure. Here is an example of the fasta file:
and an another text file is something like below:
I would like to have identifiers as following:
Thank you for your help in advance.