reporting the best alignment(s) only in pairwise2
1
2
Entering edit mode
9.5 years ago

I have a fastq file of reads, say "reads.fastq". I want to align the sequences to a string saved as a fasta file ref.faa. I am using the following code for this

reads_array = []

for x in Bio.SeqIO.parse("reads.fastq","fastq"):
    reads_array.append(x)

for x in Bio.SeqIO.parse("ref.faa","fasta"):
    refseq = x

result = open("alignments_G10_final","w")

aligned_reads = []

for x in reads_array:
    alignments =pairwise2.align.globalms(str(refseq.seq).upper(),str(x.seq),2,-1,-5,-0.05)
    for a in alignments:
        result.write(format_alignment(*a))
        aligned_reads.append(x)

But I want to report only the best alignment for each read. How can I choose this alignment from the scores in a[2]. I want to report the best alignment(s) with the highest value of score

alignment next-gen-sequencing • 2.7k views
ADD COMMENT
1
Entering edit mode
6.3 years ago

just add at the end: one_alignment_only=True

ADD COMMENT

Login before adding your answer.

Traffic: 1732 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6