I need to rename a fasta file according to a blast output file where my fasta was a query. Not all sequences had a hit.
my.fasta:
>ASSI-1_0
TTCCTTTTTGGTTCTCGATATTATGAACAGTTTCTCATCA
>ASSI-2_1
GTGAGGGAGGAGGACGCCTCGAGCAGAGGTAGGTCTGGAG
>ASSI-3_2
ATGTTAGCAAGTATAGCCAACTATATGAAACCTATGTCTT
>ASSI-4_3
GAATATCATTAAAAATCTACATTTATTTATGAGTTAGTAC
>ASSI-5_4
AGGACACCAGAAACTTTCTCCAAAGCTGAATTTGTGTATT
blast.out (one hit per query and just first two columns):
ASSI-3_2 scaf0270669_20068.102
ASSI-4_3 scaf0189112_70083.538
ASSI-5_4 scaf0083789_70072.963
ASSI-8_7 scaf0423760_50193.589
ASSI-11_10 scaf0285971_60192.428
ASSI-12_11 scaf0409557_70062.641
What I need to get is this:
>ASSI-1_0
TTCCTTTTTGGTTCTCGATATTATGAACAGTTTCTCATCA
>ASSI-2_1
GTGAGGGAGGAGGACGCCTCGAGCAGAGGTAGGTCTGGAG
>ASSI-3_2 scaf0270669_20068.102
ATGTTAGCAAGTATAGCCAACTATATGAAACCTATGTCTT
>ASSI-4_3 scaf0189112_70083.538
GAATATCATTAAAAATCTACATTTATTTATGAGTTAGTAC
>ASSI-5_4 scaf0083789_70072.963
AGGACACCAGAAACTTTCTCCAAAGCTGAATTTGTGTATT
The way I went about this is this:
Pull out fasta headers and get rid of ">":
grep "^>" my.fasta | sed -e 's/>//g' > headers.fastaCut first two columns out of blast output and only leave one hit per query:
awk '!x[$1]++' blast.out | cut -f1,2 > headers.blastThen I merge two header files but with
tacso then I can use myawkcommand to get rid of duplicates, and then I usetacagain to restore name order:tac headers.blast headers.fasta | sort -r -V | awk '!x[$1]++' | tac > headers.clean
headers.clean:
ASSI-1_0
ASSI-2_1
ASSI-3_2 scaf0270669_20068.102
ASSI-4_3 scaf0189112_70083.538
ASSI-5_4 scaf0083789_70072.963
Then replace the original headers with new headers (only works for fasta formats where sequences are in one-line):
awk 'NR%2==0' my.fasta | paste -d'\n' headers.clean -> my_renamed.fasta
The third step is what's giving me trouble. It's a very clumsy approach and I want to see if anybody has a better way?
•
written
by
san.san •
3 answers
use the join command
$ join -t ' ' -a 2 -1 1 -2 1 <( cat input.blast | sort -t ' ' -k1,1 ) <(cat input.fasta | paste - - | cut -c2- | sort -t ' ' -k1,1) | awk -F ' ' '{printf(">%s",$1);i=2;if(NF==3) {printf(" %s",$2);i++;} printf("\n%s\n",$i);}'
the options -t (join , sort) and -F (awk) require a tabulation. (Ctrl-v tab)
>ASSI-1_0
TTCCTTTTTGGTTCTCGATATTATGAACAGTTTCTCATCA
>ASSI-2_1
GTGAGGGAGGAGGACGCCTCGAGCAGAGGTAGGTCTGGAG
>ASSI-3_2 scaf0270669_20068.102
ATGTTAGCAAGTATAGCCAACTATATGAAACCTATGTCTT
>ASSI-4_3 scaf0189112_70083.538
GAATATCATTAAAAATCTACATTTATTTATGAGTTAGTAC
>ASSI-5_4 scaf0083789_70072.963
AGGACACCAGAAACTTTCTCCAAAGCTGAATTTGTGTATT
yes that's the option -a of join
-a FILENUM
also print unpairable lines from file FILENUM, where FILENUM is
1 or 2, corresponding to FILE1 or FILE2
I'm getting this when I run it:
>STE-1_0
TCTTAACTCATTGTGTGTGTATGCC
>STE-100000_218031
scaf0293423_20071.471
>STE-10000_10255
TCTGGAAAATAAAGCGGTCCACCCC
>STE-100001_218032
scaf0391581
I wonder if it's because some blast outputs are like so:
STE-121254_274055 scaf7351783
Rather than:
STE-121254_274056 scaf0143876_10040.799
I'd like to use SeqKit, only one single command.
Just download the executable binaries for your operating system (Windows/Linux/Mac OS X) and run:
$ seqkit replace -p '^([^ ]+)' -r '$1 {kv}' -k blast.out my.fasta
[INFO] read key-value file: blast.out
[INFO] 6 pairs of key-value loaded
>ASSI-1_0
TTCCTTTTTGGTTCTCGATATTATGAACAGTTTCTCATCA
>ASSI-2_1
GTGAGGGAGGAGGACGCCTCGAGCAGAGGTAGGTCTGGAG
>ASSI-3_2 scaf0270669_20068.102
ATGTTAGCAAGTATAGCCAACTATATGAAACCTATGTCTT
>ASSI-4_3 scaf0189112_70083.538
GAATATCATTAAAAATCTACATTTATTTATGAGTTAGTAC
>ASSI-5_4 scaf0083789_70072.963
AGGACACCAGAAACTTTCTCCAAAGCTGAATTTGTGTATT
I would write a biopython script to get it done, create a dictionary out of the blast output and loop over your fasta file
Log in to answer this question.