Изменение ориентации обратной последовательности в файле fasta не работает

Question

Я пытаюсь получить обратные последовательности, правильно ориентированные в файле. Это код:

import os
import sys import pysam
from Bio import SeqIO, Seq, SeqRecord

def main(in_file):
    out_file = "%s.fa" % os.path.splitext(in_file)[0]
    with open(out_file, "w") as out_handle:
        # Write records from the BAM file one at a time to the output file.
        # Works lazily as BAM sequences are read so will handle large files.
        SeqIO.write(bam_to_rec(in_file), out_handle, "fasta")

def bam_to_rec(in_file):
    """Generator to convert BAM files into Biopython SeqRecords.
    """
bam_file = pysam.Samfile(in_file, "rb")
for read in bam_file:
    seq = Seq.Seq(read.seq)
    if read.is_reverse:
        seq = seq.reverse_complement()
    rec = SeqRecord.SeqRecord(seq, read.qname, "", "")
    yield rec

if __name__ == "__main__":
    main(*sys.argv[1:])`

Когда я распечатываю обратные последовательности, код работает. Но когда в файле это распечатано как обратная последовательность. Может ли кто-нибудь помочь мне выяснить, что происходит не так? Вот ссылка на мой инфил: https://https//www.dropbox.com/sh/68ui8l7nh5fxatm/AABUr82l01qT1nL8I_XgJaeTa?dl=0

Answer 1

Обратите внимание, что уродливый счетчик просто печатает 10000 последовательностей, а не больше.

сравнивая одно без обращения к тому, которое обращает вспять при необходимости. Вот вывод на пару последовательностей, не стесняйтесь тестировать его, я думаю, ваша проблема в том, что yield возвращает итератор, но вы его не повторяете, если только я неправильно понимаю, кто вы делать:

Оригинал:

SOLEXA-1GA-2: 2: 93: 1281: 961 # 0 GGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAG

становится:

SOLEXA-1GA-2: 2: 93: 1281: 961 # 0 CTAACCCTAACCCTAACCCTAACCCTAACCTAACCC

И если не поменять

Оригинал:

SOLEXA-1GA-2: 2: 12: 96: 1547 # 0 ACACACAAACACACACACACACACACACACACACCCCC

становится:

SOLEXA-1GA-2: 2: 12: 96: 1547 # 0 ACACACAAACACACACACACACACACACACACACCCCC Вот мой код:

import os
import sys 
import pysam
from Bio import SeqIO, Seq, SeqRecord

def main(in_file):
    out_file = "%s.fa" % os.path.splitext(in_file)[0]
    with open('test_non_reverse.txt', 'w') as non_reverse:
        with open(out_file, "w") as out_handle:
            # Write records from the BAM file one at a time to the output file.
            # Works lazily as BAM sequences are read so will handle large files.
            i = 0
            for s in bam_to_rec(in_file):
                if i == 10000:
                   break
                i +=1 
                SeqIO.write(s, out_handle, "fasta")
            i = 0
            for s in convert_to_seq(in_file):
                if i == 10000:
                   break
                i +=1

                SeqIO.write(s, non_reverse, 'fasta')

def convert_to_seq(in_file):
    bam_file = pysam.Samfile(in_file, "rb")
    for read in bam_file:
        seq = Seq.Seq(read.seq)
        rec = SeqRecord.SeqRecord(seq, read.qname, "", "")
        yield rec


def bam_to_rec(in_file):
    """Generator to convert BAM files into Biopython SeqRecords.
    """
    bam_file = pysam.Samfile(in_file, "rb")
    for read in bam_file:
        seq = Seq.Seq(read.seq)
        if read.is_reverse:
            seq = seq.reverse_complement()
        rec = SeqRecord.SeqRecord(seq, read.qname, "", "")
        yield rec

if __name__ == "__main__":
    main(*sys.argv[1:])