从指定长度的 BAM/SAM 文件中提取读取答案

【问题标题】：Extract reads from a BAM/SAM file of a designated length从指定长度的 BAM/SAM 文件中提取读取
【发布时间】：2019-02-04 14:56:59
【问题描述】：

我对 Perl 有点陌生，并希望使用它来从我的 BAM（对齐）文件中提取特定长度的读取。

BAM 文件包含长度为 19 到 29 nt 的片段。以下是前 2 次读取的示例：

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 nt 的内容。

我尝试使用以下代码来做到这一点：

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

但是，程序没有给出任何结果... 谁能建议这样做的正确方法？

【问题讨论】：

标签： perl samtools

【解决方案1】：

你的问题有点混乱。代码 sn-p 应该是 Perl 脚本还是调用 Perl one-liner 的 shell 脚本？

假设您打算编写一个 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 子句。

【讨论】：

根据@stack0114106 的回答判断，您实际上可能正在寻找第 10 个字段。在 Perl 数组中，从索引 0 开始，因此您需要将上面的代码从 10 更新为 9。

【解决方案2】：

请注意，示例输入中的第 10 个字段的长度为 22 或 24。此外，您使用的语法是错误的。这是匹配长度为 22 的字段的 Perl 单行代码。

$ 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

$

【讨论】：