SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extract multiple fasta sequences from a fasta file based on sequenes entomology Bioinformatics 38 12-19-2015 07:28 PM
How to Extract Multiple Sequence from Multi Fasta File by ID list Anti Bioinformatics 4 02-02-2015 05:02 AM
Extract sequence from multi fasta file with PERL andreitudor Bioinformatics 25 01-31-2013 01:27 PM
Extract partial sequence from FASTA record cdlam Bioinformatics 9 10-30-2012 03:21 PM
Extract only sequence ids from fasta file with makeblastdb angeloulivieri Bioinformatics 13 07-30-2012 03:41 AM

Reply
 
Thread Tools
Old 09-07-2017, 03:12 PM   #1
clarissaboschi
Member
 
Location: US

Join Date: Apr 2010
Posts: 61
Default 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!
clarissaboschi is offline   Reply With Quote
Old 09-08-2017, 02:07 PM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

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
Brian Bushnell is offline   Reply With Quote
Old 09-14-2017, 01:37 PM   #3
pmiguel
Senior Member
 
Location: Purdue University, West Lafayette, Indiana

Join Date: Aug 2008
Posts: 2,237
Default

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

--
Phillip
pmiguel is offline   Reply With Quote
Old 09-15-2017, 06:08 AM   #4
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,135
Default

Quote:
Originally Posted by clarissaboschi View Post
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?
kmcarr is offline   Reply With Quote
Old 09-15-2017, 10:39 AM   #5
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

Quote:
Originally Posted by pmiguel View Post
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.
Brian Bushnell is offline   Reply With Quote
Old 09-20-2017, 02:54 PM   #6
Biocomputronics
Junior Member
 
Location: USA

Join Date: Sep 2017
Posts: 2
Default

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])
Biocomputronics is offline   Reply With Quote
Old 09-21-2017, 07:04 AM   #7
clarissaboschi
Member
 
Location: US

Join Date: Apr 2010
Posts: 61
Default

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.
clarissaboschi is offline   Reply With Quote
Old 09-21-2017, 07:36 AM   #8
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 688
Default

/*
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"]____
Richard Finney is offline   Reply With Quote
Reply

Tags
fasta, fasta file

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:39 PM.


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