Объединение файлов Fasta и Qual с различным порядком заголовков в FASTQ

Question

Я пытаюсь объединить файл fasta и файл qual в новый файл fastq, имея в виду случай, когда два файла могут иметь разный порядок в своих идентификаторах последовательности.Для этого я попытался на первом этапе моего сценария отсортировать последовательности, что прекрасно работает, когда я тестирую его как отдельный сценарий.То же самое с остальными, когда я запускаю отдельно ту часть, где он объединяет файлы в fastq, он работает отлично.Но теперь, когда я пытаюсь объединить два метода в одном скрипте, он не работает, и я не знаю, что еще делать!Буду признателен, если вы мне поможете.

Вот мой сценарий.Он создает новый файл fastq, но содержимое перепутано, а не то, что я хочу.Я запускаю его из терминала так:

$ perl script.pl reads.fasta reads.qual > reads.fq

Сценарий:

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

die ("Usage: script.pl reads.fasta reads.qual > reads.fq") unless  (scalar @ARGV) == 2;

open FASTA, $ARGV[0] or die "cannot open fasta: $!\n";
open QUAL, $ARGV[1] or die "cannot open qual: $!\n";

my $offset = 33; 
my $count = 0;
local($/) = "\n>";

my %id2seq = ();
my $id = '';
my %idq2seq = ();
my $idq = '';
my (@sort_q, @sort_f);

while(<FASTA>){
    chomp;
        if($_ =~ /^>(.+)/){
            $id = $1;
        }else{
            $id2seq{$id} .= $_;
        }
     }

for $id (sort keys %id2seq)
    {
     @sort_f = "$id\n$id2seq{$id}\n\n";
     print @sort_f;
    }

while(<QUAL>){
chomp;
    if($_ =~ /^>(.+)/){
        $idq = $1;
    }else{
        $idq2seq{$idq} .= $_;
    }
}

for $idq (sort keys %idq2seq)
    {
    @sort_q = "$idq\n$idq2seq{$idq}\n\n";
    print "@sort_q";
    }

while (my @sort_f) {
chomp @sort_f;
my ($fid, @seq) = split "\n", @sort_f;   
my $seq = join "", @seq; $seq =~ s/\s//g;
my $sortq = @sort_q;
chomp my @sortq;
my ($qid, @qual) = split "\n", @sortq;

@qual = split /\s+/, (join( " ", @qual));
# convert score to character code:
my @qual2 = map {chr($_+$offset)} @qual;
my $quals = join "", @qual2; `enter code here`
die "missmatch of fasta and qual: '$fid' ne '$qid'" if $fid ne $qid;
$fid =~ s/^\>//;
print STDOUT (join( "\n", "@".$fid, $seq, "+$fid", $quals), "\n");
$count++;
}
close FASTA;
close QUAL;
print STDERR "wrote $count entries\n";

Заранее спасибо

Answer 1

Прошло много времени с тех пор, как я использовал perl, но я бы подошел к этому, используя хэш пар ключ / значение как для быстрого ввода, так и для качественного ввода.Затем запишите все пары, перебирая хеш-код fastta и вытаскивая соответствующую строку качества.

Я написал в python что-то, что будет делать то, что вам нужно, вы можете увидеть это в действии здесь :

Предполагается, что ваш ввод выглядит следующим образом:
reads.fasta

>fa_0
GCAGCCTGGGACCCCTGTTGT
>fa_1
CCCACAAATCGCAGACACTGGTCGG

reads.qual

>fa_0
59 37 38 51 56 55 60 44 43 42 56 65 60 68 52 67 43 72 59 65 69
>fa_1
36 37 47 72 34 53 67 41 70 67 66 51 47 41 73 58 75 36 61 48 70 55 46 42 42

вывод

@fa_0
GCAGCCTGGGACCCCTGTTGT
+
;%&387<,+*8A<D4C+H;AE
@fa_1
CCCACAAATCGCAGACACTGGTCGG
+
$%/H"5C)FCB3/)I:K$=0F7.**
@fa_2
TCGTACAGCAGCCATTTTCATAACCGAACATGACTC
+
C?&93A@:?@F,2:'KF*20CC:I7F9J.,:E8&?F

---

import sys


# Check there are enough arguments
if len(sys.argv) < 3:
    print('Usage: {s} reads.fasta reads.qual > reads.fq'.format(s=sys.argv[0]), file=sys.stderr)
    sys.exit(1)

# Initalise dictionaries for reads and qualities
read_dict = dict()
qual_dict = dict()

fa_input = sys.argv[1]
qual_input = sys.argv[2]

# Read in fasta input
with open(fa_input, 'r') as fa:
    for line in fa:
      line = line.strip()
      if line[0] == '>':
        read_dict[line[1:]] = next(fa).strip()
      else:
        next(fa)

# Read in quality input
with open(qual_input, 'r') as qual:
    for line in qual:
      line = line.strip()
      if line[0] == '>':
        qual_dict[line[1:]] = next(qual).strip()
      else:
        next(qual)

count = 0
# Iterate over all fasta reads
for key, seq in read_dict.items():
    # Check if read header is in the qualities data
    if key in qual_dict.keys():
        # There's both sequence and quality data so write stdout
        read_str = '@{header}\n{seq}\n+\n{qual}'.format(
            header=key,
            seq=seq,
            qual=''.join([chr(int(x)) for x in qual_dict[key].split(' ')]))
        print(read_str, file=sys.stdout)
        count += 1
    else:  # not found
        # Write error to stderr
        print('Error: {k} not found in qual file'.format(k=key), file=sys.stderr)

# Print count to stderr
print('{c} reads written'.format(c=count), file=sys.stderr)

Если вам нужно использовать смещение для показателя качества, отредактируйте
qual=''.join([chr(int(x)) for x in qual_dict[key].split(' ')])) в
qual=''.join([chr(int(x) + offset) for x in qual_dict[key].split(' ')])) и определите переменную offset перед этим.