3

I have the following problem:

My data sheet looks like this

AAAGGGTTT AAAGGG
AAAGGGCCC GGGCCC

and I want to display the difference between the two sequences in the third row like this

AAAGGGTTT AAAGGG TTT
AAAGGGCCC GGGCCC AAA

what I've tried is to use diff. I extracted the individual sequences in files (f1.txt and f2.txt) and formatted them so they could be compared line by line with diff which created the problem that it only works when the start of the sequence is similar (Line 1 of datasheet).

awk '{gsub(".","&\n");printf "%s",$0}' < f1.txt >f1a.txt
awk '{gsub(".","&\n");printf "%s",$0}' < f2.txt >f2a.txt
 
diff -y f1a.txt f2a.txt 

does anyone have an idea how to achieve the desired output?

  • 1
    Does a "sequence" contain the space? Or is that second column something else that the problem you've described can just ignore? – roaima Aug 4 at 13:02
8

Is this what you are after?

awk '{$3=$1;sub($2,"",$3)}1' file
  • $3=$1 copies the 1st field to the 3rd field and

  • sub($2,"",$3) looks for the 2nd field in the 3rd field. If there is a match, sub substitutes the matched string by a null string in the 3rd field.

  • The 1 at the end prints the result. It is equivalent to a {print} statement, so you could rewrite it as {$3=$1;sub($2,"",$3);print}.

The result:

AAAGGGTTT AAAGGG TTT
AAAGGGCCC GGGCCC AAA
| improve this answer | |
  • Yes that's exactly what I was looking for. Thank you so much! Would you mind explaining your syntax a bit more for me? What does the "" and the 1 at the end do? I'm new to bash and always thankful for more insight – Daniel R Aug 5 at 13:31
  • @DanielR Glad to help, see the edit for clarification. Also, bear in mind this has only to do with awk. Bash does not do anything special here. – Quasímodo Aug 5 at 17:45
4

If you need to display pairwise alignments of sequences, use the appropriate bioinformatics tools. The sequence alignment will show differences in formats expected by the typical users. Here, for pairwise alignment of nucleotide sequence, you can use, for example, blastn from the BLAST tools (which can be installed, for example, using conda), see this command from the thread: https://www.biostars.org/p/18087/#18095 :

blastn -query querySeqSet.fasta -subject targetSeqSet.fasta
| improve this answer | |
  • Idk anything about these tools, but +1 for right tool for the job – D. Ben Knoble Aug 5 at 3:27
  • thanks for your answer but I am looking specifically for the bases that are not aligned – Daniel R Aug 5 at 13:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.