I have a fasta file with multiple sequences comprising 38 sequences. The length of the sequences are around 11000 bp.
How can i get changes in the genomes based in a reference genome? (aa subtitutions or deletions)
Bioinformatics Stack Exchange is a question and answer site for researchers, developers, students, teachers, and end users interested in bioinformatics. It only takes a minute to sign up.
Sign up to join this communityI have a fasta file with multiple sequences comprising 38 sequences. The length of the sequences are around 11000 bp.
How can i get changes in the genomes based in a reference genome? (aa subtitutions or deletions)
Python3 via the difflib
library/package. AlignIO
can be leveraged (Biopython) to construct the align
list. To be honest if your coding isn't strong you could just cut and paste the alignment into this code. Nucleotides not amino acids ... same idea just easier to type. It will work for aa and gaps. "Ref" (output) refers to the first alignment position only.
import difflib
align=['AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGGTGATATGA', 'AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGAAGATATGA', 'AGCTAGAGGATATGA', 'AGCTAGATCATATGT']
ref= 'AGCTAGATGATATGA'
for query in align:
print('{} (ref) => {} (query)'.format(ref,query))
for loc,pos in enumerate(difflib.ndiff(ref, query)):
if pos[0]==' ':
continue
elif pos[0]=='-':
print(u'Ref "{}"'.format(pos[-1]))
elif pos[0]=='+':
print(u'Query has "{}" at position {}'.format(pos[-1],loc))
print ('')
Output,
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)
AGCTAGATGATATGA (ref) => AGCTAGGTGATATGA (query)
Ref "A"
Query has "G" at position 7
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)
AGCTAGATGATATGA (ref) => AGCTAGAAGATATGA (query)
Ref "T"
Query has "A" at position 8
AGCTAGATGATATGA (ref) => AGCTAGAGGATATGA (query)
Ref "T"
Query has "G" at position 8
AGCTAGATGATATGA (ref) => AGCTAGATCATATGA (query)
Ref "G"
Query has "C" at position 9
Query has "T" at position 16
Sequences must be aligned, i.e. ref sequence within the global alignment and separated. Hope that makes sense. It will work with multiple mutations between ref and query as the last example points out (two mutations per example query was just too much to type).
A reproducible example in future please.