Извлечение операций чтения из файла BAM/SAM заданной длины

Я немного новичок в 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)    
}

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


perl samtools
person pkom    schedule 03.02.2019    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Ваш вопрос немного сбивает с толку. Фрагмент кода должен быть сценарием 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 никогда не выполняется.

person Stefan Becker    schedule 03.02.2019
comment
Судя по ответам @stack0114106, вы, возможно, ищете 10-е поле. В Perl массив начинается с индекса 0, поэтому вам нужно будет обновить приведенный выше код с 10 на 9. - person Stefan Becker; 04.02.2019

arrow_upward
1
arrow_downward

Обратите внимание, что 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

$
person stack0114106    schedule 04.02.2019