Fasta Length

Question

[k.nirmalraman]

I was wondering if there are any pre built packages in R or Python/ scripts, where one can extract FASTA sequences that are longer than N bases or shorter or exactly N bases from a multi fasta file.

[Biomonika (Noolean)]

I believe biopieces can do that as well https://code.google.com/p/biopieces/wiki/HowTo#Howto_filter_sequences_on_length

Answer 1

[Ashutosh Pandey]

I think this post should help you.

How to Filter Multi fasta by length??

Answer 2

[brentp]

https://github.com/lh3/bioawk

bioawk -c fastx '(length($seq) > 60){ print ">"$name"\n"$seq }' input.fasta > output.fasta

Answer 3

[Neilfws]

One option for R is the seqinr package.

Use read.fasta() to read in a multi-sequence fasta file:

library(seqinr)
f <- read.fasta("myfile.fa")  # default is DNA; add ", seqtype = "AA" for protein

Use summary() for a range of sequence statistics, including length:

summary(f[[1]])  # first sequence
summary(f[[1]])$length  # length of first sequence

Use lapply() to get sequences longer than N bases:

f1 <- f[which(lapply(f, function(x) summary(x)$length) > 60)]  # N > 60
# ignore warnings; they relate to sequence composition

There is also sure to be a R/Bioconductor solution using Biostrings.

[Manu Prestat]

Shortened version using the getLength() dedicated function:

f1 <- f[which(getLength(f)>60)]

[Neilfws]

Nice; I always forget the get*() functions in seqinR.