CODEWITHSUNDEEP
pythonpython-2.7pysam
Loki
Loki

Reputation: 99

Changing an orientation of a reverse sequence in the fasta file doesnt work

I am trying to get the reverse sequences orientated correctly in a file. This is the code:

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

When I print out the reverse sequences, the code works. But when in the file it is printed out as a reverse sequence. Can anyone help me to find out what is going wrong? Here is the link to my infile: https://www.dropbox.com/sh/68ui8l7nh5fxatm/AABUr82l01qT1nL8I_XgJaeTa?dl=0

Upvotes: 0

Views: 387

Answers (1)

E.Serra
E.Serra

Reputation: 1574

Note the ugly counter is just to print 10000 sequences, not more.

comparing one without ever reversing with one that reverses if needed Here's the output on a couple of seqs, feel free to test it, I think your issue is that yield returns an iterator but you are not iterating it, unless I am missunderstanding what you are doing:

Original:

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

Becomes:

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

And if not reverse:

Original:

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

Becomes:

SOLEXA-1GA-2:2:12:96:1547#0 ACACACAAACACACACACACACACACACACACCCCC Here's my code:

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

Upvotes: 1

Related Questions