perl 从 2 个数组中提取公共元素（fastq 文件中的公共序列）

Question

我有 2 个文件（paired_end），读取来自同一生物体的 fastq 格式（File_R1.fastq 和 File_R2.fastq）。我想使用 bwa 和 samtools （生物信息学）确定深度覆盖，为此我需要两个具有相同读取次数、相同名称的文件（file_1：@XX00341:4450:6341 1:N:0:AACGTTAA+ TTGCAATT and file_2: @XX00341:4450:6341 2:N:0:AACGTTAA+TTGCAATT)，在两个文件中的顺序相同，问题是两个文件都有空读取（每个文件对应的读取，可能在一个，但不是另一个！！！！），如果我尝试用 bwa 运行它，它会由于空序列而失败。

所以我正在尝试制作一个 perl 脚本来提取两个文件中具有相同名称和顺序的 NO 空读取。

这是fastq文件的格式（每次读取有4行：名称（@XX00341:19:000H27K25:1:11101:4450:6341 1:N....），序列（GAGGTGCGTGGTTGTCACCCC...）， Q_head (+) 和 Quality (FFFFFFFFFFFFFFF....)，按此顺序。

file_r1.fastq（以 = .... 1:N:0:AACGTTAA+TTGCAATT 为特征）

@XX00341:4450:6341 1:N:0:AACGTTAA+TTGCAATT
CGCCCATCATGAGGTGCGTGGTTGTCACCCCCGCAAACGCGGAGGTGTAAACAGGCTCACCTGTGGGTTGTTGGTAGTCGTTATTGTGCTTGCGCTGTTCATCTGGATTACTGGATTCGACGCGTTGAGC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@XX00341:14420:6341 1:N:0:AACGTTAA+TTGCAATT
GCTTTGTCGTCGTCGGTTTTAAAGTGAACCGCTTTACCTGTTTCTGCTTGAACTTGTTCTGCTTGAGGTGCTGCTTCTGGTTTTGTTTCTTCTTTCTGACAACCTACCGCTAGCATTACTGTTGCAGCAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF/FFFFAAFFFFFFAFFFFAFFFFFFFFFFAFAFFFFFAFFFFFFFFFFFFFFFFAFFFFFFFFFFFFFAFFFFFFFAFFFFFFF/FFFFFFFFFFFF
@XX00341:10259:6342 1:N:0:AACGTTAA+TTGCAATT

+

@XX00341:6685:6342 1:N:0:AACGTTAA+TTGCAATT
CATTGGAAGGCAAGCCAGAACAAGGCAAGAATATTCCAGATGACATCGTGCGCGTTCGCATCGATCGCAACAGCGGTCTGCTGACTCATAAAGTGGATAGCACGTCCATGTTTGAATACTTTGAGAAAGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF6FFFFAFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFAFFFFFFAFFFFFFFFFFFFF
@XX00341:4051:6343 1:N:0:AACGTTAA+TTGCAATT
TCACCATGATCGGATTTATGAATGGTTTAGTGGACAGCATGATCAAAAATGCGATTGCTTGGCAAACCAGCCATTTGCAGATTCATCAAAGTGCTTATCTCGTTAATCCTGAACTGAAAGACACCATACC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@XX00341:12307:6344 1:N:0:AACGTTAA+TTGCAATT
TCCAAATTGAGAGTTAATGTGAAAAAAGAATGACTTTTTGCTTGTCATGCAGCGGATTTGTGTGATACCACTAATGACGCAAACATTTTCGTAGAACATTACACAATACTATTTAACGAAAAAAGAACGA
+
FAAFFAFAAFFAFFF/FF/FFAFFFFFFFFFAAFFF=F=AFFFFFFFFFFFFAFFFFFFFFFFAFFFFFFFFFFAFFFFFAAFFFFFFFFFF/FAFFFFAFFFFFFFFFFFFAAAFFAF///FAFFFFFF
@XX00341:24250:6345 1:N:0:AACGTTAA+TTGCAATT
CGAACATAGAGCAAGCTCTGGTAAGTCCCGATGGGAGCTTAAAGACACAAGTAGACCAAAAGCTCGCAGAACATGGCTTGAGCCGTAAAGTCACCGTGGCGTCGCGCAACTTTCTCACCATTCGGCATCT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

file_r2.fastq（以 = .... 2:N:0:AACGTTAA+TTGCAATT 为特征）

@XX00341:4450:6341 2:N:0:AACGTTAA+TTGCAATT
GCGGCTTTGGTAGCGAAGCGCGTCTACCAACCGCAGCCATGAAACAACTGGCGTTTGAAGTGGAAAAAACGGCAGCGGGCAGTATTCCGGTACTGATTGAAGCCATCGAACAAGTGGCGGTGCCTACCGC
+
AFFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFF/FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFAFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@XX00341:14420:6341 2:N:0:AACGTTAA+TTGCAATT
GCGGTTCGAGTGGCCAACGTTGAACTTCATGCTCGGTAAAAAAGCAACCATTTAACGTGGTGATGACAATTAAATATAGGAATAAATGAGAAATTCTTTGCATCACATCTAATAATCCGGTTTGTGTCTT
+
AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF/FFFAFAFFFFFFFF/FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFFF
@XX00341:10259:6342 2:N:0:AACGTTAA+TTGCAATT
TCTTTTTTCGTTAAATAGTATTGTGTAAT
+
FFFFFFFFFFFFFAFFFFFFFFFFFFFFF
@XX00341:6685:6342 2:N:0:AACGTTAA+TTGCAATT
GAATAAATGCTGTCTTCGAGGCTGTTACCAACGTATTCGGTTGGCTCCGTGCCTTTCTCAAAGTATTCAAACATGGACGTGCTATCCACTTTATGAGTCAGCAGACCGCTGTTGCGATCGATGCGAACGC
+
AFFFFFFFFAFFFFAFFFFFFFFFFFFFFFFFFFFFFFF/FF/FFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFF=FF=FFFFFFFAFFFFFFF
@XX00341:4051:6343 2:N:0:AACGTTAA+TTGCAATT
GAAGCTATCATGCCATCGGCGAGAAACCGCTCCGATACGGCTTTCACGCTTTGATGCTTGTCTAACGTTGTCACGATGCTTTGCGAGTCAGGTATGGTGTCTTTCAGTTCAGGATTAACGAGATAAGCAC
+
AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@XX00341:12307:6344 2:N:0:AACGTTAA+TTGCAATT
GTTCATCTGTCTGTCGTTCTTTTTTCGTTAAATAGTATTGTGTAATGTTCTACGAAAATGTTTGCGTCATTAGTGGTATCACACAAATCCGCTGCATGACAAGCAAAAAGTCATTCTTTTTTCACATTAA
+
AFFFFFFFFFAFFFFFFFFAFFFFFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFFFFF=FFFFFFFAFFFFFFFAFFFFFFFFFF=FAFFFFFFFFAFFFAFFAFFFFFFFFA=AFF6AFFAFAFF
@XX00341:24250:6345 2:N:0:AACGTTAA+TTGCAATT
CCATAGCATGGCAATATCAAAATCCGCAACGGCGATCGGCGGCTTAACAGTAATGAGATCGTCTGAAAAGCCCTCCGCCAACGCCATTTTTTCTGGCACGATGGCAATGAGATCGCGCTTTTTCAATAGA
+
AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF/FFFFFFFFFFFFFFFFFFFFFFFFFFF

在这种情况下，file_1.fastq (@XX00341:10259:6342 1:N....) 中的一次读取是空的，第二和第四行是空的（序列和质量），但不是在 file_2 .fastq (@XX00341:10259:6342 2:N....);在此示例中，必须省略两个序列（每个文件中的 4 行）！！！

这是我要完成的代码：

#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;

my ($fastQ_R1, $fastQ_R2);
      GetOptions (
            'R1|r1=s'   =>\$fastQ_R1,
            'R2|r2=s'   =>\$fastQ_R2,

      );

    sub extract_list {
        my ($file_in) = @_;
        open FILE, '<', $file_in or die "cant open the $file_in\n";

        my (@elements, @list, @seq);
        while(
            defined(my $head    =   <FILE>) &&   # 1 line
            defined(my $seq     =   <FILE>) &&   # 2 line
            defined(my $qhead   =   <FILE>) &&   # 3 line
            defined(my $quality =   <FILE>)      # 4 line
        ){
            if ($seq=~ m/^$/g){
                next;
            }
            else {  push @seq, $head;   }
         }
        close FILE;
        foreach (@seq){
            chomp;
            @elements = split '\s', $_;
                push @list, $elements[0]; # split to eliminate 1:N... (file_1) and 2:N... (file_2)
        }
        return @list;
    }

    my @list_R1 = extract_list ($fastQ_R1);

    my @list_R2 = extract_list  ($fastQ_R2);

到目前为止，我在两个文件中都有没有空序列（@list_R1 和 @list_R2）的列表，两个数组的比较将给出两个文件（@ common_elements) 从两个数组的比较中获得。

@list_R2=

@XX00341:4450:6341   
@XX00341:14420:6341   
@XX00341:6685:6342   
@XX00341:4051:6343   
@XX00341:12307:6344   
@XX00341:24250:6345

@list_R2=

@XX00341:4450:6341   
@XX00341:14420:6341   
@XX00341:10259:6342
@XX00341:6685:6342   
@XX00341:4051:6343   
@XX00341:12307:6344   
@XX00341:24250:6345

@common_elemnts=

@XX00341:4450:6341   
@XX00341:14420:6341   
@XX00341:6685:6342   
@XX00341:4051:6343   
@XX00341:12307:6344   
@XX00341:24250:6345

因此，新数组 (@common_elemnts) 将用于搜索和提取每个文件的常见读取（4 行），并将提供两个输出文件：Files_R1_common.fastq（从 File_R1.fastq 获得）和 File_R2_common.fastq（从 File_R2.fastq 获得）。

任何建议，将不胜感激！！！！ 非常感谢

Answer 1

您最好使用哈希映射来检查列表中是否存在项目。
构造公共哈希图：

my %list1_map = map { $_, 1 } @list_R1;
my %common_map = map { $list1_map{$_} ? ($_, 1) : () )} @list_R2;

然后处理您的列表：

for my $item(@list_R1) {
if (defined ($common_map{$item})) {
    # ... process ...       
}

Answer 2

在以下解决方案中，我已将 FASTQ 解析移至单独的类，该类将从文件中读取 4 行并返回代表一条记录的 id、hdr、seq 和 qual 的 hashref。我把这门课写给你作为练习。这样做使下面的逻辑易于理解。可读性很重要。

此代码仅保存来自 file1 的空序列的 ID（即 $hdr=~/^@(\S+)/ 中的 $1）。然后它读取file2并输出完整的记录，并保存空记录的ID。最后，您再次读取 file1 并通过删除由 %filtered 哈希指示的记录来输出完整记录。

此解决方案效率低下，因为它读取 file1 两次。原始问题中没有提到 FASTQ 文件通常包含数百万条记录，因此将 file1 中的所有未过滤记录保存在 %filtered 哈希中将需要大量 RAM，但如果您有大量 RAM，则可以更改以适应，但就我个人而言，我不会担心两次读取文件。

# SAVE IDS OF EMPTY RECORDS FROM FILE1 IN HASH
my %filter;
my $fq = new Fastq($file1);
while (my $rec = nextSeq($fq))
{
    $rec->{seq} or $filter{$rec->{id}} = undef; # not using value
}

# ADD IDS OF EMPTY RECORDS FROM FILE2 TO HASH AND OUTPUT FILTERED_FILE2
open(my $out, '>', "$file2.filtered") or die($!);
$fq = new Fastq($file2);
while (my $rec = nextSeq($fq) )
{
    if ( $rec->{seq} )
    {
        # READ2 FOR THIS PAIR IS NOT EMPTY
        if ( ! exists($filter{$rec->{id}}) )
        {
             # READ1 FOR THIS PAIR IS ALSO NOT EMPTY
             print $out join("\n", $rec->{hdr}, $rec->{seq}, '+', $rec->{qual}), "\n";
        }
    } else
    {
         # READ2 IS EMPTY, ADD TO FILTERED HASH
         $filtered{$rec->{id}} = undef;
    }
}
close($out);

# OUTPUT FILTERED_FILE1
open($out, '>', "$file1.filtered") or die($!);
$fq = new Fastq($file1);
while (my $rec = nextSeq($fq) )
{
    if ( $rec->{seq} and ! exists($filter{$rec->{id}}) )
    {
         print $out join("\n", $rec->{hdr}, $rec->{seq}, '+', $rec->{qual}), "\n";
    }
}  
close($out);