как написать файл fastq из другого файла - PullRequest
Новая
       50

как написать файл fastq из другого файла

1 голос
/ 15 февраля 2020

Меня попросили прочитать из двух файлов (чтение слева и справа) Aip02.R1.fastq и Aip02.R2.fastq, и получить чередуемый файл fastta, используя функцию zip. Левый и правый файлы были файлами fastq, но когда я сжал их вместе, чтобы создать новый файл fastq, функция записи больше не работает. Это дает мне ошибку "SeqRecord (id =) имеет недопустимую последовательность." 

  #!/usr/bin/env python3
  # Import Seq, SeqRecord, and SeqIO
  from Bio import SeqIO
  from Bio.SeqRecord import SeqRecord
  from Bio.Seq import Seq
  leftReads = SeqIO.parse("/scratch/AiptasiaMiSeq/fastq/Aip02.R1.fastq", "fastq")
  rightReads = SeqIO.parse("/scratch/AiptasiaMiSeq/fastq/Aip02.R2.fastq","fastq")
  A= zip(leftReads,rightReads)
  SeqIO.write(SeqRecord(list(A)), "interleave.fastq", "fastq")

1 Ответ

1 голос
/ 16 февраля 2020

Ваша прямая и обратная последовательности, вероятно, имеют одинаковый идентификатор. Поэтому используйте следующий код, чтобы добавить суффикс к идентификаторам. Я использовал /1 и /2 здесь, но также используются такие вещи, как .f и .r. 

from Bio import SeqIO
import itertools

def interleave(iter1, iter2) :
    for (forward, reverse) in itertools.izip(iter1, iter2):
        assert forward.id == reverse.id
        forward.id += "/1"
        reverse.id += "/2"
        yield forward
        yield reverse

leftReads = SeqIO.parse("/scratch/AiptasiaMiSeq/fastq/Aip02.R1.fastq", "fastq")
rightReads = SeqIO.parse("/scratch/AiptasiaMiSeq/fastq/Aip02.R2.fastq","fastq")

handle = open("interleave.fastq", "w")
count = SeqIO.write(interleave(leftReads, rightReads), handle, "fastq")
handle.close()
print("{} records written to interleave.fastq".format(count))

Этот код может стать медленным для больших файлов fastq. Например, см. здесь , где сообщается, что для создания файла fastq размером 2 ГБ требуется 14 минут. Таким образом, они дают этот улучшенный способ: 

from Bio.SeqIO.QualityIO import FastqGeneralIterator
import itertools

file_f = "/scratch/AiptasiaMiSeq/fastq/Aip02.R1.fastq"
file_r = "/scratch/AiptasiaMiSeq/fastq/Aip02.R2.fastq"

handle = open("interleave.fastq", "w")
count = 0

f_iter = FastqGeneralIterator(open(file_f,"rU"))
r_iter = FastqGeneralIterator(open(file_r,"rU"))
for (f_id, f_seq, f_q), (r_id, r_seq, r_q) in itertools.izip(f_iter,r_iter):
    assert f_id == r_id
    count += 2
    #Write out both reads with "/1" and "/2" suffix on ID
    handle.write("@%s/1n%sn+n%sn@%s/2n%sn+n%sn" 
                 % (f_id, f_seq, f_q, r_id, r_seq, r_q))
handle.close()
print("{} records written to interleave.fastq".format(count)
...