Entering edit mode
9.5 years ago
sumithrasank75
▴
140
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