Script To Print Number Of Occurences Of Genes From Multifasta File

Question

[Syed Imtiyaz]

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.

[Neilfws]

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.

Answer 1

[Rm]

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
....

[brentp]

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

[Syed Imtiyaz]

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

[Rm]

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...

[Neilfws]

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.

[Syed Imtiyaz]

ENSTGUG00000012287|ENSTGUT00000012814|1475 ACCGGTGCCAGGGGCCGCGGTTGGCTGCGAAGCGGCGGCTCCCGCCCCCTGCGGAATCAGCCCCAGGTCCGGGGCGGCTCTACCTGCCGGCACGATGAACCTCACCGCCGAGAGCCACCGCATTCCGCTGAGCGACGGCAACAGCATCCCGCTCTTGGGGCTGGGCACCTACGCCGACCCGCAGAAAACTCCCAAAGGTTCCTGTCTGGAGGCGGTGAAGATTGCCATCGATGCTGGTTACCGCCACATCGACGGTGCCTTTGTCTACTTCAATGAGCATGAAGTGGGACAAGCCATCCGGGAGAAGATTGCTGAAGGGAAGATCAAGAGAGAAGACATATTTTACTGTGGCAAGCTGTGGAATACCTGCCACCCCCCAGAGCTGGTGCGTCCCACACTGGAGAAAACCCTGAAGATCCTGCAGCTGGACTACGTTGACCTCTACATTATTGAGCTGCCAATGGCTTTCAAGCCTGGAGATGCACTCTACCCAAAAGATGAAAATGGAAAATTTATCTACCATGAGACAGACTTATGTGCCACTTGGGAGGCTCTG

[Rm]

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

[Syed Imtiyaz]

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.

[Syed Imtiyaz]

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

[Rm]

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

[Syed Imtiyaz]

Hi, Still the problem is same

[Syed Imtiyaz]

i will send u the whole file plz test it

[Syed Imtiyaz]

ENSTGUG00000012287|ENSTGUT00000012814|1475 ACCGGTGCCAGGGGCCGCGGTTGGCTGCGAAGCGGCGGCTCCCGCCCCCTGCGGAATCAGCCCCAGGTCCGGGGCGGCTCTACCTGCCGGCACGATGAACCTCACCGCCGAGAGCCACCGCATTCCGCTGAGCGACGGCAACAGCATCCCGCTCTTGGGGCTGGGCACCTACGCCGACCCGCAGAAAACTCCCAAAGGTTCCTGTCTGGAGGCGGTGAAGATTGCCATCGATGCTGGTTACCGCCACATCGACGGTGCCTTTGTCTACTTCAATGAGCATGAAGTGGGACAAGCCATCCGGGAGAAGATTGCTGAAGGGAAGATCAAGAGAGAAGACATATTTTACTGTGGCAAGCTGTGGAATACCTGCCACCCCCCAGAGCTGGTGCGTCCCACACTGGAGAAAACCCTGAAGATCCTGCAGCTGGACTACGTTGACCTCTACATTATTGAGCTGCCAATGGCTTTCAAGCCTGGAGATGCACTCTACCCAAAAGATGAAAATGGAAAATTTATCTACCATGAGACAGACTTATGTGCCACTTGGGAGGCTCTG

[Rm]

upload the file and share the link

[Syed Imtiyaz]

http://www.4shared.com/folder/GgfoEnew/_online.html here is the file plz check it out thanks

[Syed Imtiyaz]

Hi, I think there is some problem in downloading the file plz go to this link http://uploading.com/files/7bbe159c/examplefasta.fa/ THank you

[Rm]

same problem downloading it..its blocked by websense

[Syed Imtiyaz]

http://www.fileserve.com/file/kts2Z3h/examplefasta.fa this should work in think

[Syed Imtiyaz]

Hi RM, the above link is working?

[Rm]

Content blocked by your organization

[Rm]

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

[Syed Imtiyaz]

http://dl.dropbox.com/u/43445136/examplefasta.fa

[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

[Syed Imtiyaz]

THank you for all your help :)

Answer 2

[Martin A Hansen]

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

[Syed Imtiyaz]

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

Thank you

[Martin A Hansen]

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

[Syed Imtiyaz]

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