I am working on Sequence Read Archive (SRA) datasets. I have one fasta file containing the SRA identifier and the sequence (atcagtt..) and one blast output file including sequence identifier, gi number and the rest of the details. Both of them are in csv format.
I want to use the common sequence identifier to combine my files. Eventually I would like to have one file with all the sequence information as well as the sequence itself.
Attached is link for a possible solution, however the sort function does only give me the sequence (as in aaccgttc...), but no identifiers.
Dear Bilgenur,
Have you tried the UNIX join command?
http://unixhelp.ed.ac.uk/CGI/man-cgi?join
Depending on size of files, it might be handy to use the R command "merge". For example something like:
allt
If I understand properly, you are dealing with a fasta header replacement task.
If so, I can provide you a Perl script that works so:
replace_fasta_headers.pl myseqs.fasta namelist.txt
Format of namelist.txt:
----------------------
oldname1 newname1
oldname2 newname2
.
.
(Columns are separated by Tabs. White spaces are allowed in newnames, but not in
oldnames !)
In your case, would be the original SRA IDs, and the IDs+the descriptions from the BLAST output file.
Tell me, if you are interested in getting the script
Dear Joan, yes I did try it. But sorting is the first issue with my files. I get the following errors if I proceed the join function without sorting.
is not sorted: SRR......113_FCD114LACXX:1:1101:4451:2028_length=76 gi|326200233:1395-2093 100.00 89 00 1 89 79 167 5e-40 165
The manual says that I can use "join -t ''" if 'sort' has no options.
So I did run this command: join -t file1 file2 >out
However I get this error: join: multi-character tab ‘file1'
My files do not have the same headings. Since the (only common thing between the two files) identifiers are not sorted, I assume that join function may not work at all.
Hej,
I am sure we can make it work if you post sample of two input files..
You can try similar python code as attached ...
With regard,
## sample fasta file
1,ATG
2,CTA
### sample blast File
1,xyz,cct,..
2,abc,jkl,..
### sample output
1,ATG,1,xyz,cct,..
2,CTA,2,abc,jkl,..
Dear Vikas Gupta,
##sample fasta file (not separated by comma but by space)
>SRR1005707.1_FCD114LACXX:1:1101:1243:2019_length=91
CTGGGCCGCACGCGCGCTACACTGAAGGAATCAGCGTGTGCCTCCCTGGGCCGAAAGGCC
CGGGTAATCCGCTGAACCTCCTGCGTGTTG
>SRR1005707.2_FCD114LACXX:1:1101:1232:2098_length=91
GGCCCAACTGAAGAAGATGTGCCCCGGCTGCGTCTCCCGCGACGAATGAAACGACCGCAA
ACTCTTCCACCGAGCTCGACGAAAAACAAAC
## sample blast file
SRR1005707.37_FCD114LACXX:1:1101:2564:2157_length=91 gi|302747135|gb|GU194765.1 98.89 90 1 0 1 90 169 80 6e-39 161
SRR1005707.113_FCD114LACXX:1:1101:4451:2028_length=91 gi|326200233:1395-2093 100.00 89 0 0 1 89 79 167 5e-40 165
in this case, I want to obtain both the sequence and the rest of the information seen in sample blast file. However my files are not sorted, when I tried to sort them it kind of crashed my computer. I was hoping to join these files without a need for sorting.
Since the information is not separated by comma, I do think that we need to change the following in the py script that you've attached:
tokens = line.split(',')
Right?
You are right, we just need to replace the delimiter while reading in the fasta file. It does not require any sorting but if there are multiple blast hits for your fasta query sequences then you will also get multiple lines in the output. At the end, numbers of line in output file will be same as in the blast input file.
Best,
Vikas
Thanks so much for the help. I have a few question regarding to the script.
1) What this mean? hash[tokens[0]] = tokens[1]
Because this line generated the following error: Traceback (most recent call last):
File "temp.py", line 7, in
hash[tokens[0]] = tokens[1]
IndexError: list index out of range
2) Also, it is not clear to me what this line is supposed to do: if tokens[0] in hash: print tokens[0] +','+hash[tokens[0]]+','+line
First error is probably because you have empty lines (I updated the code).
hash[tokens[0]] = tokens[1] means that store all the fasta sequences using fasta_id (first column) as key.
Then for second point, it is checking the blast file if it has the fasta_id as first column. Given a successful match, final line print both fasta file information together with blast file information.
As I see, the original query name (sequence ID) is the first field of each line of the blast result file, followed by the hit ID (second field) than by score values separated by white spaces.
1. Let's say that the name of your query sequence file is , and the name of your blast result file is
2. The following command will create a namelist file from your blast result file (by dupicating the first field):
awk '{col = $1 "\t" $1; $1 = col; print }' blastres.txt > namelist.txt
3. Use the attached Perl script to add the blast information to each fasta header:
replace_fasta_headers.pl myseqs.fasta namelist.txt
A new fasta file with augmented headers will be created under the name .
No previous sorting is necessary. Fasta headers without blast hits will not be changed.
I had problems by attaching the Perl script as a file.
You find the script below. Copy the content between the xxxxxxx lines, save it as
and make it executable.
xxxxxxxxxxxxxxxx
#!/usr/bin/perl
# Found in old times on Internet. Don't remember where.
#Made little modifications on it. IN
use strict;
use warnings;
@ARGV == 2 or die "usage: $0 \n";
my ( $fasta_file, $header_file ) = @ARGV;
my $destination = $fasta_file . 'headers-replaced.fasta';
open IN2, '', $destination or die "Can't write to file $destination:
+$!\n";
while ( ) {
if ( /^>(.+)$/ && exists $head_seqs{ $1 } ) {
$_ = ">$head_seqs{ $1 }\n";
}
print OUT $_;
}
close IN1;
close OUT;
xxxxxxxxxxxxxxxxx
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics