Извлечь чтение из файла BAM/SAM указанной длины

Question

Я немного новичок в Perl и хочу использовать его для извлечения чтений определенной длины из моего файла BAM (выравнивания).

Файл BAM содержит чтения, длина которых составляет от 19 до 29 нт. Вот пример первых двух прочтений:

YT:Z:UUA00182:193:HG2NLDMXX:1:1101:29884:1078   0   3R  6234066 42  22M *   0   0   TCACTGGGCTTTGTTTATCTCA  FF:FFFF,FFFFFFFF:FFFFF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:22   

YT:Z:UUA00182:193:HG2NLDMXX:1:1101:1777:1094    16  4   1313373 1   24M *   0   0   TCGCATTCTTATTGATTTTCCTTT    FFFFFFF,FFFFFFFFFFFFFFFF    AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:24   

Я хочу извлечь только те, которые, скажем, 21 нт в длину.

Я пытаюсь сделать это с помощью следующего кода:

my $string = <STDIN>;    
$length = samtools view ./file.bam | head | perl -F'\t'  -lane'length @F[10]';    
if ($length == 21){    
        print($string)    
}        

Однако программа не дает никакого результата... Может ли кто-нибудь предложить правильный способ сделать это?

Answer 1

Обратите внимание, что 10-е поле в вашей выборке имеет длину 22 или 24. Кроме того, синтаксис, который вы используете, является неправильным. Вот Perl-однострочник, чтобы сопоставить поле с длиной =22.

$ cat pkom.txt
YT:Z:UUA00182:193:HG2NLDMXX:1:1101:29884:1078   0   3R  6234066 42  22M *   0   0   TCACTGGGCTTTGTTTATCTCA  FF:FFFF,FFFFFFFF:FFFFF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:22

YT:Z:UUA00182:193:HG2NLDMXX:1:1101:1777:1094    16  4   1313373 1   24M *   0   0   TCGCATTCTTATTGATTTTCCTTT    FFFFFFF,FFFFFFFFFFFFFFFF    AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:24

$ perl -lane ' print if length($F[9])==22 ' pkom.txt
YT:Z:UUA00182:193:HG2NLDMXX:1:1101:29884:1078   0   3R  6234066 42  22M *   0   0   TCACTGGGCTTTGTTTATCTCA  FF:FFFF,FFFFFFFF:FFFFF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:22

$

Answer 2

Ваш вопрос немного сбивает с толку. Должен ли фрагмент кода быть сценарием Perl или сценарием оболочки, который вызывает Perl-однострочник?

Предполагая, что вы намеревались написать сценарий Perl, в который вы передаете вывод samtools view чтобы:

#!/usr/bin/perl
use strict;
use warnings;

while (<STDIN>) {
    my @fields = split("\t", $_);

    # debugging, just to see what field is extracted...
    print "'$fields[10]' ", length($fields[10]), "\n";

    if (length($fields[10]) eq 21) {
        print $_;
    }
}

exit 0;

С вашими данными испытаний в dummy.txt Я получил:

# this would be "samtools view ./file.bam | head | perl dummy.pl" in your case?
$  cat dummy.txt | perl dummy.pl
'FF:FFFF,FFFFFFFF:FFFFF' 22
'FFFFFFF,FFFFFFFFFFFFFFFF' 24

Ваши тестовые данные не содержат образец с длиной 21, поэтому if пункт никогда не выполняется.