This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Script To Print Number Of Occurences Of Genes From Multifasta File

Hi, i need help to write a script that will take for input a multi fasta file and output the gene names and the number of times the gene is found in the file in two columns.

fasta scripting
• 7,371 views • link
updated by Rm Rm • written by Syed Imtiyaz • Syed Imtiyaz

It will be useful if you indicate whether the answers are helpful. When we get this type of "please write my code for me" question, I often wonder whether the answer even means anything to the questioner. If your problem is that you know nothing about scripting, my advice is to go away and learn some.

• 0 views • link
written by Neilfws Neilfws

2 answers

simplest will be: Asuming fasta header as ">gene_name description..."

grep "^>" multi_fasta.txt | sed 's/>//' | awk '{print $1}' | sort | uniq -c | awk '{print $2 "\t" $1}' >gene_count.txt

Edit: using the link provided

curl http://dl.dropbox.com/u/43445136/examplefasta.fa |  grep "^>" | sed 's/>//' | awk -F"|" '{print $1}' | sort | uniq -c | awk '{print $2 "\t" $1}'

Output

ENSTGUG00000000002      1
ENSTGUG00000000010      1
ENSTGUG00000000018      1
ENSTGUG00000000021      1
ENSTGUG00000000026      1
ENSTGUG00000000027      1
ENSTGUG00000000029      1
ENSTGUG00000000037      1
ENSTGUG00000000043      1
....
• 1 views • link
written by Rm Rm

could also be a single awk command something like: awk '(/^>/){ a[substr($0, 2)]++ }END { for(header in a){ print a[header], header }}'

• 0 views • link
written by brentp brentp

thanks alot for spending ur precious time on my problem RM. But there is some problem with the commands i guess coz it's giving the gene name but it is not doing the counting part as expected could you plz check it out thank you

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

to test it, take off last part of the awk to see if uniq -c is giving out put or not...: BTW can you paste the sample file by editing your Question...

• 0 views • link
written by Rm Rm

This code works for me using a simple test fasta file. If it doesn't work for you, the problem might lie with your data.

• 0 views • link
written by Neilfws Neilfws

ENSTGUG00000012287|ENSTGUT00000012814|1475 ACCGGTGCCAGGGGCCGCGGTTGGCTGCGAAGCGGCGGCTCCCGCCCCCTGCGGAATCAGCCCCAGGTCCGGGGCGGCTCTACCTGCCGGCACGATGAACCTCACCGCCGAGAGCCACCGCATTCCGCTGAGCGACGGCAACAGCATCCCGCTCTTGGGGCTGGGCACCTACGCCGACCCGCAGAAAACTCCCAAAGGTTCCTGTCTGGAGGCGGTGAAGATTGCCATCGATGCTGGTTACCGCCACATCGACGGTGCCTTTGTCTACTTCAATGAGCATGAAGTGGGACAAGCCATCCGGGAGAAGATTGCTGAAGGGAAGATCAAGAGAGAAGACATATTTTACTGTGGCAAGCTGTGGAATACCTGCCACCCCCCAGAGCTGGTGCGTCCCACACTGGAGAAAACCCTGAAGATCCTGCAGCTGGACTACGTTGACCTCTACATTATTGAGCTGCCAATGGCTTTCAAGCCTGGAGATGCACTCTACCCAAAAGATGAAAATGGAAAATTTATCTACCATGAGACAGACTTATGTGCCACTTGGGAGGCTCTG

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

@syed: its working for me with the sequence you provided.: output ENSTGUG00000012287|ENSTGUT00000012814|1475 1

• 0 views • link
written by Rm Rm

thanks for checking. But for the count part it is not giving the expected it is always giving 1 as the count for each and every gene.

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

could u give me ur email id r some thing by taht i can send u the original file

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

check the headers for all the sequences if they are consisitant; Try replacing awk '{print $1}' with awk -F"|" '{print $1}'

• 0 views • link
written by Rm Rm

ENSTGUG00000012287|ENSTGUT00000012814|1475 ACCGGTGCCAGGGGCCGCGGTTGGCTGCGAAGCGGCGGCTCCCGCCCCCTGCGGAATCAGCCCCAGGTCCGGGGCGGCTCTACCTGCCGGCACGATGAACCTCACCGCCGAGAGCCACCGCATTCCGCTGAGCGACGGCAACAGCATCCCGCTCTTGGGGCTGGGCACCTACGCCGACCCGCAGAAAACTCCCAAAGGTTCCTGTCTGGAGGCGGTGAAGATTGCCATCGATGCTGGTTACCGCCACATCGACGGTGCCTTTGTCTACTTCAATGAGCATGAAGTGGGACAAGCCATCCGGGAGAAGATTGCTGAAGGGAAGATCAAGAGAGAAGACATATTTTACTGTGGCAAGCTGTGGAATACCTGCCACCCCCCAGAGCTGGTGCGTCCCACACTGGAGAAAACCCTGAAGATCCTGCAGCTGGACTACGTTGACCTCTACATTATTGAGCTGCCAATGGCTTTCAAGCCTGGAGATGCACTCTACCCAAAAGATGAAAATGGAAAATTTATCTACCATGAGACAGACTTATGTGCCACTTGGGAGGCTCTG

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

upload the file and share the link

• 0 views • link
written by Rm Rm

same problem downloading it..its blocked by websense

• 0 views • link
written by Rm Rm

Content blocked by your organization

• 0 views • link
written by Rm Rm

upload it to dropbox or other more safe site...

• 0 views • link
written by Rm Rm

i ran the script its working fine, all your sequences are present only once...Iam editing my answer to include part of the result

• 0 views • link
written by Rm Rm

This can be done using Biopieces www.biopieces.org) like this:

read_fasta -i test_big.fna -n 10 |
count_vals -k SEQ_NAME |
uniq_vals -k SEQ_NAME |
write_tab -ck SEQ_NAME_COUNT,SEQ_NAME -x

Cheers,

Martin

• 1 views • link
written by Martin A Hansen Martin A Hansen

Thanks for reply, actually i am not aware of 'Biopieces'. Could you please suggest how to use it

Thank you

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

The above does exactly what you wanted, but do have a look at the website and the documentation there.

• 0 views • link
written by Martin A Hansen Martin A Hansen

thanks for ur suggestion maasha, i checked it is not installed on the server which i am working. could u plz try it in perl or awk

THank you

• 0 views • link
written by Syed Imtiyaz • Syed Imtiyaz

Log in to answer this question.

More posts like this