Как определить количество контигов, которые выше, благодаря указанному размеру c

Question

Надеюсь, у вас все хорошо

пожалуйста, у меня есть файл фаста, например

>contig1
sequence
>contig2
sequence
>contig3
>sequence

, у каждой последовательности есть длина, я хочу определить число контигов, превышающее 9000 ( поэтому длина последовательности превышает 9000)

Спасибо

Answer 1

Вы пометили grep, поэтому grep -c '.\{9000\}' your_fasta.fa, вероятно, самый простой метод.

Более «биоинформатический» подход заключается в использовании seqkit (https://bioinf.shenwei.me/seqkit/): seqkit seq -m 9000 your_fasta.fa > newfile.txt для извлечения последовательностей более 9000 оснований в 'newfile.txt' и grep -c ">" newfile.txt для подсчета количества последовательностей длиной> 9000.

Кроме того, вот несколько решений awk / perl / bioawk, которые вы можно адаптировать: https://www.biostars.org/p/79202/

Answer 2

Вы можете выполнить эту задачу, если установите Bio Perl модуль Bio :: SeqIO . Затем вы можете сохранить приведенный ниже скрипт как count_contigs.pl в том же каталоге, что и файл с contigs с именем " contigs.fasta ", и запустить скрипт с perl count_contigs.pl. Он будет считать контиги длиной более 9000 б.п. из входного файла и печатать результат.

#!/usr/bin/perl
use strict;
use warnings;    
use Bio::SeqIO;

# Setting minimum length to be more than 9000
my $min_len = 9000;

# Reading the input fasta file
my $seqio_in = Bio::SeqIO->new(-file => "contigs.fasta", 
                                     -format => "fasta" );
# Setting the counter
my $counter = 0; 

# Counting sequences if length > min_len     
while ( my $seq = $seqio_in->next_seq ) {
    if ( $seq->length  >  $min_len ) {
        $counter++;
    }
}

# Print the result
print "There are '$counter' sequences that are longer than $min_len\n";