Re^3: Size of sequences in fastafile

It's still not very clear precisely what you expect but here's my final stab at this.

use strict;
use warnings;

use Test::More tests => 1; # Present this SSCCE as a test

my @seq; # Sequences longer than $min are stored here
my @want = (
    'GGAGGTCTTTAGCTTTAGGGAAACCC',
); # These are the sequeneces we expect for the given data set
my $min = 15; # Minimum length of any sequence to be considered

my $this = '';
while (<DATA>) {
    chomp;
    if (/^>/) {
        push @seq, $this if $min <= length $this;
        $this = '';
    } else {
        $this .= $_;
    }
}
push @seq, $this if $min <= length $this;

is_deeply \@seq, \@want; # Check that our algorithm has worked

__DATA__
>NM_001 Homo sapiens ADA2 (CECR1)
GATCCAA
>NM_002 Homo sapiens IKBKG
GGAGGTCTTTAGCTTTAGGGAAACCC
[download]

Comment on Re^3: Size of sequences in fastafile Download Code

Replies are listed 'Best First'.
Re^4: Size of sequences in fastafile by Sofie (Acolyte) on Mar 01, 2020 at 12:03 UTC
This is my data: >NM_030643.4 Homo sapiens apolipoprotein L4 (APOL4) GAGGTGCTGGGGAGCAGCGTGTTTGCTGTGCTTGATTGTGAGCTGCTGGGAAGTTGTGACTTTCATTTTA CCTTTCGAATTCCTGGGTATATCTTGGGGGCTGGAGGACGTGTCTGGTTATTATATAGGTGCACAGCTGG AGGTGAGATCCACACAGCTCAGACCAGCTGGATCTTGCTCAGTCTCTGTCAGAGGAAGATCCCTTGGAGG AGGCCCCGCAGCGACATGGAGGGAGCTGCTTTGCTGAAAATCTTTGTCGTCTGCATCTGGAACCAAAATC >NM_001198855.1 Homo sapiens cytochrome P450 family 2 subfamily C memb +er 8 (CYP2C8) ACATGTCAAAGAGACACACACTAAATTAGCAGGGAGTGTTATAAAAACTTTGGAGTGCAAGCTCACAGCT GTCTTAATAAGAAGAGAAGGCTTCAATGGAACCTTTTGTGGTCCTGGTGCTGTGTCTCTCTTTTATGCTT CTCTTTTCACTCTGGAGACAGAGCTGTAGGAGAAGGAAGCTCCCTCCTGGCCCCACTCCTCTTCCTATTA >NR_029834.1 Homo sapiens microRNA 200a (MIR200A), microRNA CCGGGCCCCTGTGAGCATCTTACCGGACAGTGCTGGATTTCCCAGCTTGACTCTAACACTGTCTGGTAAC GATGTTCAAAGGTGACCCGC >AC067940.1 Homo sapiens clone RP11-818E9, LOW-PASS SEQUENCE SAMPLING AAATACAACTTTAAATCAAAACGGTAAAAATTCCACTCTTTCATACTAACTTCAAAAGTATTTGCTTTAA AAAAAAAGNNNNNNNNNNAAACTGAATTTCTATTAAGCATCTATTTATAGAAGAGAGTAAACACCCCGTG AATAAAAGACAGAGAATTGTAGCAGCCCGAAGTCCCTTTTCTCTCCTCCCAAGCATTTGGCTCTGGTCCA AATTCACATATCCTGCTCCGTAAAACAAAGTGCCTTGGTTAACCTAACGTTATTCCTTGAACAGTAGTTT AGTGATCAACTAGTTTTTGTTGTTGTTGTTGTTTGAGACAGAGTCTCACTCTGTCGCCCAGGCTGGAGTG CAGTGGCGAGATCTCAGCTCACTGCAACCTCTGCTGCCCAGGTTCAAGGGATTCTCCTGCCTCAGCCTCC CAAGTAGCTGGTATTACAGGCACCTGCCACCGCGCCTGGCTAATTTTTTTTTTTTTTTTTTTTTGTATTT [download] The question is, how can I find out the sequence length of each sequence and extract those of a specific length? I have tried both to put the data in an array, which puts each line in one element, and also to put it in a scalar and tried to use the > to separate each seq, but I am stuck.	[reply] [d/l]
Re^5: Size of sequences in fastafile by zubenel0 (Sexton) on Mar 01, 2020 at 12:39 UTC
This code should work but you will need to install BioPerl module Bio::SeqIO. `use Bio::SeqIO; # Setting minimum length to 250 my $min_len = 250; # Reading the input fasta file my $seqio_in = Bio::SeqIO->new(-file => "Genes.fasta", -format => "fasta" ); # Creating the output fasta file my $seqio_out = Bio::SeqIO->new(-file => ">Genes_filt_250.fasta", -format => "fasta" ); # Saving sequences to the output if length > min_len while ( my $seq = $seqio_in->next_seq ) { if ( $seq->length > $min_len ) { $seqio_out->write_seq($seq); } }` [download] The sequences with length higher than 250 should be saved in file "Genes_filt_250.fasta": >NM_030643.4 Homo sapiens apolipoprotein L4 (APOL4) GAGGTGCTGGGGAGCAGCGTGTTTGCTGTGCTTGATTGTGAGCTGCTGGGAAGTTGTGAC TTTCATTTTACCTTTCGAATTCCTGGGTATATCTTGGGGGCTGGAGGACGTGTCTGGTTA TTATATAGGTGCACAGCTGGAGGTGAGATCCACACAGCTCAGACCAGCTGGATCTTGCTC AGTCTCTGTCAGAGGAAGATCCCTTGGAGGAGGCCCCGCAGCGACATGGAGGGAGCTGCT TTGCTGAAAATCTTTGTCGTCTGCATCTGGAACCAAAATC >AC067940.1 Homo sapiens clone RP11-818E9, LOW-PASS SEQUENCE SAMPLING AAATACAACTTTAAATCAAAACGGTAAAAATTCCACTCTTTCATACTAACTTCAAAAGTA TTTGCTTTAAAAAAAAAGNNNNNNNNNNAAACTGAATTTCTATTAAGCATCTATTTATAG AAGAGAGTAAACACCCCGTGAATAAAAGACAGAGAATTGTAGCAGCCCGAAGTCCCTTTT CTCTCCTCCCAAGCATTTGGCTCTGGTCCAAATTCACATATCCTGCTCCGTAAAACAAAG TGCCTTGGTTAACCTAACGTTATTCCTTGAACAGTAGTTTAGTGATCAACTAGTTTTTGT TGTTGTTGTTGTTTGAGACAGAGTCTCACTCTGTCGCCCAGGCTGGAGTGCAGTGGCGAG ATCTCAGCTCACTGCAACCTCTGCTGCCCAGGTTCAAGGGATTCTCCTGCCTCAGCCTCC CAAGTAGCTGGTATTACAGGCACCTGCCACCGCGCCTGGCTAATTTTTTTTTTTTTTTTT TTTTGTATTT [download]	[reply] [d/l] [select]
Re^4: Size of sequences in fastafile by Sofie (Acolyte) on Mar 01, 2020 at 11:03 UTC
Thank you very much for your efforts. I am obviously not able to explain very well what I am trying to do. This seems to put all the sequences in one single element of the @seq array, and the length includes the header ( not just the seq data). I just want to extract all the sequences with a minimal length. But the sequences need to still be in fasta format. How can I use the regex for > to extract the sequence in between the headers and check the length? Since each line of sequence has a newline in between. Thanks again and sorry for not being able to explain. /Sofie	[reply]