SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
KaKs calculator help (merging fasta files) Shorash Bioinformatics 4 07-28-2013 08:19 AM
Tools to generate VCF from two FASTA, or mutant FASTA from Ref FASTA and VCF? jeffseq Bioinformatics 3 05-28-2013 10:59 AM
merging fastq running bwa gives different results from running bwa then merging sam. bw. Bioinformatics 1 06-25-2012 03:08 AM
Make Multi-fasta with annotations detq182 Bioinformatics 4 04-15-2012 05:50 AM
merging a tab and a fasta file arg General 2 10-21-2010 10:53 AM

Reply
 
Thread Tools
Old 07-22-2014, 06:59 AM   #1
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Lightbulb Merging annotations and fasta headings

I have a fasta file like this:
Quote:
>comp0_c0_seq1 len=204 path=[4976:0-203]
ATTGTACTTAATCTA.....
>comp0_c1_seq1 len=222 path=[8835:0-221]
GTACAATCACGGCACATT......
and an annotation file like this:
Quote:
comp1558_c0_seq1 repressor protein
comp8142_c0_seq1 aspartate aminotransferase
comp8357_c0_seq1 l-xylulose reductase
comp8387_c0_seq1 succinyl- synthetase alpha
comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase
And I want to replace the fasta files` heading by the annotation to produce:
Quote:
>comp0_c0_seq1 xxx protein
ATTGTACTTAATCTA.....
>comp0_c1_seq1 xxx reductase
GTACAATCACGGCACATT......
Here is what I wrote with python trying to do this:
Quote:
f = open("test")
ano = open("annotation.txt")
output = open("merged.fasta",'w')
for line in ano:
x = line.spilt('\t')
for line in f:
str = line.split(' ')
if x[0] == str[0]:
output.write(x)
else:
output.write(line)
But it does not work, says "AttributeError: 'str' object has no attribute 'spilt'
"
... any better ideas to do this or found any way to fix the codes?

Thanks
Simonli is offline   Reply With Quote
Old 07-22-2014, 07:45 AM   #2
Wallysb01
Senior Member
 
Location: San Francisco, CA

Join Date: Feb 2011
Posts: 286
Default

Hum, looks like you’ve type “spilt” when you mean “split”.

You could also convert you fasta file into tab format then just use a join command. Ie it could be:

File1:
Code:
someheader1	GATC….
someheader2	ATCG….
File2:
Code:
someheader1	othername1
someheader2	othername2
So then it would be:
Code:
sort -k1,1 File1 > File1.sort
sort -k1,1 File2 > File2.sort
join -j1 -o1.1,2.2,1.2 -e “NA” -a1 File1.sort File2.sort | sed ’s/ /_/‘ | see ’s/ /\t/‘ > Newfile
That will give you something that would look like this that could be transformed back in fasta format:

Code:
someheader1_othername1	GATC...
someheader2_othername2	ATCG...
And careful with “\t” in that code above. In case you don’t know, you’ll probably need to replace it with Ctrl-V-tab, to get it Unix to recognize it.

Also you can you the fastx toolkit to convert to and from tab format.
Wallysb01 is offline   Reply With Quote
Old 07-22-2014, 07:48 AM   #3
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

This is a syntax error. You need to use the proper function name. You have two similarly named functions -- 'spilt' and 'split'. Only one of them is giving you an error.

[edit: I see that Wallysb01 beat me to the answer ... and in a much nicer fashion :-) ]
westerman is offline   Reply With Quote
Old 07-22-2014, 10:15 AM   #4
crazyhottommy
Senior Member
 
Location: Gainesville

Join Date: Apr 2012
Posts: 140
Default

you may use awk:
http://stackoverflow.com/questions/1...n-a-dictionary

Usage: awk -f foo.awk dict.dat user.dat
http://www.gnu.org/software/gawk/man...Functions.html
http://www.gnu.org/software/gawk/man...de/Arrays.html

NR == FNR {
rep[$1] = $2
next
}

{
for (key in rep) {
gsub(key, rep[key])
}
print
}
crazyhottommy is offline   Reply With Quote
Old 07-22-2014, 10:42 AM   #5
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Default

Thank you Wallysb01,
But the thing is, it might need some levels of "header matching",
the headers in annotation file are not completed in the fasta file..

Quote:
someheader1 othername1
someheader3 othername2


[QUOTE]
someheader1
ATTTTCTG...
someheader2
ATCGGG...
someheader3
ATCTTG...
someheader4
ATCGAG...[QUOTE]

by the way, my code still not work after fix the syntax error...
then I wrote a R code trying to do the same thing,,no error message comes out, but the output all containing the " "..

Quote:
ano <- read.table("annotation.txt",sep="\t")
fasta <- readLines("TrinityS_condensed.fasta")
i<-1
j<-1
while (i<5000){
ano1 <- paste(">",ano[i,1],sep="")
while (j<40604){
str <- strsplit(fasta[j]," ")
item <- str[[1]]
if (ano1 == item[1]) {
fasta[j] <- paste(fasta[j],ano[i,2])
}
j <- j+2
}
i <- i+1
}
dput(fasta,file="mergeoutput")
Simonli is offline   Reply With Quote
Old 07-22-2014, 11:23 AM   #6
crazyhottommy
Senior Member
 
Location: Gainesville

Join Date: Apr 2012
Posts: 140
Default

I am sorry, the awk does not work in this case, awk use $1 for first word and $2 for the second word.

I just wrote a python script below, and I tested it, it worked fine.
myfasta file

>comp1558_c0_seq1 len=204 path=[4976:0-203]
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 len=222 path=[8835:0-221]
GTACAATCACGGCACATT......
>comp8570_c0_seq1 len=232 path=[3344:0-232]
GCATCGCTTATCGGTTTA......

annotation file:

comp1558_c0_seq1 repressor protein
comp8142_c0_seq1 aspartate aminotransferase
comp8357_c0_seq1 l-xylulose reductase
comp8387_c0_seq1 succinyl- synthetase alpha
comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase

script:
Code:
with open ("C:/Users/Tang Ming/Desktop/anotation.txt", "r") as annotation:
    anotation_dict = {}
    for line in annotation:
        line = line.split()
        if line: #test whether it is an empty line
            anotation_dict[line[0]]=line[1:]
        else:
            continue

# really should not parse the fasta file by myself. there are
# many parsers already there. you can use module from Biopython 
ofile = open ("C:/Users/Tang Ming/Desktop/output.txt", "w")

with open ("C:/Users/Tang Ming/Desktop/my_fasta.txt", "r") as myfasta:
    for line in myfasta:
        if line.startswith (">"):
            line = line[1:] # skip the ">" character
            line = line.split()
            if line[0] in anotation_dict:
                new_line = ">" + str(line[0]) + " " + " ".join(anotation_dict[line[0]])
                ofile.write ( new_line + "\n")
        else:
            ofile.write(line)
                
                
ofile.close() # always remember to close the file.
output:
>comp1558_c0_seq1 repressor protein
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 aspartate aminotransferase
GTACAATCACGGCACATT......
>comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase
GCATCGCTTATCGGTTTA......
crazyhottommy is offline   Reply With Quote
Old 07-22-2014, 01:16 PM   #7
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Default

thank you crazyhottommy for you codes,
it works but somehow the lines were not right..
in fasta file it should looks like 1,3,5,7,9..... lines are the headings and 2,4,6,8,10.... lines are ATCGs...

and in the output there are ATCGs take one line and some take multiple lines.

I`m new to python so not sure what was wrong...
Simonli is offline   Reply With Quote
Old 07-22-2014, 01:42 PM   #8
crazyhottommy
Senior Member
 
Location: Gainesville

Join Date: Apr 2012
Posts: 140
Default

I see, just changed the second last line code to:

ofile.write(line + "\n")

let me know if it works or not.

https://gist.github.com/crazyhottomm...2b48e6ad3f4630
crazyhottommy is offline   Reply With Quote
Old 07-28-2014, 10:40 AM   #9
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Default

it seems it still doesn`t work, the multi-lines issue still exists...
It is probably caused by the maximum line length limit while read in the data---now sure how to fix it.
Simonli is offline   Reply With Quote
Old 07-31-2014, 07:50 AM   #10
crazyhottommy
Senior Member
 
Location: Gainesville

Join Date: Apr 2012
Posts: 140
Default

Quote:
Originally Posted by Simonli View Post
it seems it still doesn`t work, the multi-lines issue still exists...
It is probably caused by the maximum line length limit while read in the data---now sure how to fix it.
see my updated code. you need to add "\n" at the last line of the code.
https://gist.github.com/crazyhottomm...2b48e6ad3f4630
crazyhottommy is offline   Reply With Quote
Old 08-01-2014, 01:56 AM   #11
syfo
Just a member
 
Location: Southern EU

Join Date: Nov 2012
Posts: 103
Default

That simple awk command works fine:
Code:
awk 'NR==FNR{des[">"$1]=$0;next}/^>/ && des[$1]{$0=">"des[$1]}1' annotation_file fasta_file
cat annotation_file
Code:
comp1558_c0_seq1 repressor protein
comp8142_c0_seq1 aspartate aminotransferase
comp8357_c0_seq1 l-xylulose reductase
comp8387_c0_seq1 succinyl- synthetase alpha
comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase 
comp0_c1_seq1 bli
comp0_c0_seq1 blahblah blah blah
cat fasta_file
Code:
>comp0_c0_seq1 len=204 path=[4976:0-203]
ATTGTACTTAATCTA.....
>comp0_c1_seq1 len=222 path=[8835:0-221]
GTACAATCACGGCACATATCTATGCA
CACAGTCAGCTAC
resulting output:
Code:
>comp0_c0_seq1 blahblah blah blah
ATTGTACTTAATCTA.....
>comp0_c1_seq1 bli
GTACAATCACGGCACATATCTATGCA
CACAGTCAGCTAC
syfo is offline   Reply With Quote
Old 08-11-2014, 06:09 AM   #12
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Default

Thank you syfo,
but I`m not sure how to save the output into a txt file... awk command seems just print the result.
Simonli is offline   Reply With Quote
Old 08-11-2014, 06:28 AM   #13
syfo
Just a member
 
Location: Southern EU

Join Date: Nov 2012
Posts: 103
Default

Quote:
Originally Posted by Simonli View Post
Thank you syfo,
but I`m not sure how to save the output into a txt file... awk command seems just print the result.
Just redirect the default output into a file by adding at the end of the command something like

> name-of-your-output-file.txt
syfo is offline   Reply With Quote
Old 08-11-2014, 06:30 AM   #14
syfo
Just a member
 
Location: Southern EU

Join Date: Nov 2012
Posts: 103
Default

Quote:
Originally Posted by syfo View Post
Just redirect the default output into a file by adding at the end of the command something like

> name-of-your-output-file.txt
which gives you more explicitly:

Code:
awk 'NR==FNR{des[">"$1]=$0;next}/^>/ && des[$1]{$0=">"des[$1]}1' annotation_file fasta_file > output-file.txt

Last edited by syfo; 08-11-2014 at 06:31 AM. Reason: typo
syfo is offline   Reply With Quote
Old 08-11-2014, 06:32 AM   #15
Simonli
Member
 
Location: Clemson, SC

Join Date: Oct 2013
Posts: 14
Default

You are awesome! Thank you so much!
Simonli is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:25 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO