0

我正在尝试在文件中正确定位反向序列。这是代码：

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:])`

当我打印出反向序列时，代码有效。但是当在文件中时，它会以相反的顺序打印出来。谁能帮我找出问题所在？这是我的 infile 的链接： https ://www.dropbox.com/sh/68ui8l7nh5fxatm/AABUr82l01qT1nL8I_XgJaeTa?dl=0

4

1 回答 1

1

请注意，丑陋的计数器只是打印 10000 个序列，而不是更多。

比较一个没有反转的和一个反转的如果需要这是几个seq的输出，随意测试它，我认为你的问题是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 ACACACAAACACACACACACACACACACACACCCCC

变成：

SOLEXA-1GA-2:2:12:96:1547#0 ACACACAAACACACACACACACACACACACACCCCC 这是我的代码：

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:])
于 2018-08-01T12:31:13.873 回答