SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)
-   -   How to extract the longest sequence from a fasta file (http://seqanswers.com/forums/showthread.php?t=77978)

clarissaboschi 09-07-2017 03:12 PM

How to extract the longest sequence from a fasta file
 
I have a fasta file with many sequences, an I would like to get just the longest one (does not matter the id). How can I do that?

I would like to have as output a fasta file with only the longest sequence.

thanks!

Brian Bushnell 09-08-2017 02:07 PM

You can make a read length histogram with BBMap's readlength.sh (or reformat.sh lhist), then run:

Code:

reformat.sh in=file.fa out=longest.fa minlength=X

pmiguel 09-14-2017 01:37 PM

Must. Not. Write. Perl. One. Liner.

--
Phillip

kmcarr 09-15-2017 06:08 AM

Quote:

Originally Posted by clarissaboschi (Post 210766)
I have a fasta file with many sequences, an I would like to get just the longest one (does not matter the id). How can I do that?

What about the case where there is more than one sequence which match the longest length? Save all? Only the first or last?

Brian Bushnell 09-15-2017 10:39 AM

Quote:

Originally Posted by pmiguel (Post 210941)
Must. Not. Write. Perl. One. Liner.

--
Phillip

Go right ahead :)

Quote:

Originally Posted by kmcarr
What about the case where there is more than one sequence which match the longest length? Save all? Only the first or last?

These are interesting questions that make me wonder whether the original problem is worth solving. I'm guessing the OP wants the main chromosome from a perfect bacterial assembly, but perhaps the Clarissa should chime in on what she's trying to do.

Biocomputronics 09-20-2017 02:54 PM

This sort of task is nearly trivial in Python, it is well worth your while to learn it and the important packages such as BioPython which make jokes like seen in this thread all the more funny.

So here is a quick and dirty way. A little thought could make it more elegant and/or generalizable. This would just parse through your multifasta file and spit back the longest sequence.

Code:

from Bio import SeqIO

myList = []

for seq_record in SeqIO.parse("test.fa", "fasta"):
    myList.append([seq_record.id, str(seq_record.seq), len(seq_record)])

myList.sort(key=lambda x: x[2])

print("the longest sequence is:")
print(">", myList[-1][0], sep='')
print(myList[-1][1])


clarissaboschi 09-21-2017 07:04 AM

Thanks for the suggestions.
I have a list of protein ids, and I need to retrieve the respective nucleotide sequence. I am using eutils to retrieve them based on these ids, but for some ids, I have many sequences for one id. So for this specific case, I am selecting the longest sequence.

Richard Finney 09-21-2017 07:36 AM

/*
how to compile: gcc -Wall -O2 -o longestseq longestseq.c
example: ./longestseq example.fa

*/

#include <stdlib.h>
#include <string.h>
#include <stdio.h>

char s[100000]; // input buffer
int main(int argc, char *argv[])
{
FILE *fp;
int i;
size_t cur_spot;
size_t last_spot;
size_t best_spot;
size_t cursize;
size_t bestsofarsize;

if (argc != 2) { fprintf(stderr,"Error: need inputfilefasta. This program finds first longest fasta sequence(just one, no same length "ties"). Usage: ./longestseq example.fa\n"); exit(0); }
fp = fopen(argv[1],"r");
best_spot = cur_spot = last_spot = ftell(fp);
cursize = bestsofarsize = 0;
while (fgets(s,99999,fp))
{
for (i=0;s[i];i++) { if ((s[i]=='\r')||(s[i]=='\n')) {s[i] = (char)0;break;}} //handle dos or unix
if (s[0]=='>')
{
last_spot = cur_spot;
cursize = 0;
continue;
}
cursize += strlen(s);
if (cursize > bestsofarsize)
{
bestsofarsize = cursize;
best_spot = last_spot;
}
cur_spot = ftell(fp);
}

fseek(fp,best_spot,SEEK_SET);
fgets(s,99999,fp);
printf("%s",s);
while (fgets(s,99999,fp))
{
if (s[0] == '>') break;
printf("%s",s);
}
fclose(fp);
return 0;
}


[FONT="Courier New"][SIZE="1"]____


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

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.