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

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

0 голосов
/ 03 февраля 2019

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

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

Ответы [ 2 ]

0 голосов
/ 04 февраля 2019

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

$
0 голосов
/ 03 февраля 2019

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

...