Extract mutations from fasta sequences

Question

[dasilvauirajan]

I have a large amount of align protein sequences in the .fasta forma, and a reference sequence, every of that has the same length. I would like to extract only the amino acid mutations from these sequences, so that, in the end, I want to have a list that looks something like this: I456L, W675T, etc . Is there a software or any way to do this? Thankful

[GenoMax]

Pierre has a complete solution but in case that does not work you could use blastp with -outfmt 3 which will identify the difference and output it so.

Query_1    181  VAATMCIGPEGDLHGVPPGECAVRLVKAGASIIGVNCHFDPTISLKTVKLMKEGLEAARL  240
Subject_1  181  ..............................V.............................  240

Biopython blast parser may be able to help finish the rest.

Answer 1

[Pierre Lindenbaum]

Using bioalcidaejdk; http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html and a fasta file where the very first sequence is the reference:

java -jar dist/bioalcidaejdk.jar -e 'FastaSequence ref=null; while(iter.hasNext()) { final FastaSequence seq = iter.next(); if(ref==null) { ref=seq; } else { for(int i=0;i< seq.length() && i< ref.length();i++) { char aa1 = ref.charAt(i); char aa2 = seq.charAt(i); if(aa1!=aa2) println(seq.getName()+"\t"+aa1+(i+1)+aa2); } } }'  input.fasta

[dasilvauirajan]

Pierre Lindenbaum I cannot find the fold dist, and bioalcidaejdk.jar too. There is a file named bioalcidaejdk.java in the bioalcidae folder.

[GenoMax]

Did you follow install instructions:

Requirements / Dependencies

java compiler SDK 11. Please check that this java is in the ${PATH}. Setting JAVA_HOME is not enough : (e.g: https://github.com/lindenb/jvarkit/issues/23 )

Download and Compile

$ git clone "https://github.com/lindenb/jvarkit.git"
$ cd jvarkit
$ ./gradlew bioalcidaejdk

[dasilvauirajan]

Thank you both genomax and Pierre Lindenbaum. I had a problem with my JDK compiler, just solved it, i get running and everything went well!!!