0

I have two files of each with 3 sequences (200 length each) with a header like this:

#File1

>1  
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTT  
>2  
TGGGAAAACAGCATTCAGCGGTGGCTTATTCCTGCTAAGGATGTTGGCCGCATTCATGCTGAGCACAACCTCGACGGCCTGCTGAGGGGCGATTCGGCATCCCGCGCTGCCTTTATGAAGGCAATGGGAGAGGCAGGGCTACGCACCATCAACGAGATGCGACGAACGGACAACCTCCCGCCATTGCCGGGTGGCGATGT  
>3  
GAAATGGGAACCGCGAACATGCCTGCACATCCGTTTGTGCGACCCGCTTACGATACTCGCGAGGAAGAGGCCGCCAGCGTCGCCATTGCCAGGATGAATCAGGCTATTGATGAGGTATTGAGCAAGTGAATGAAGATAATATCTACGCCTTGCTTTCTCCCCTGGCAGAAGGACGGGTATATCCCTATGTTGCGCCATTA
#File2

>1  
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCGGGAAGGTGCAGCATAATGTGCTTTGGATT  
>2  
TGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCATCCATATACTGTATAGACCATCTGAGCGTT  
>3  
TACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTGCAGGAATGAGTTTTTGTACGACCGTATTCGCGCAGCTTTACTTCAATTTTGTGCTGTTTGCTCAGCTTCGTGAAAGAGGCCTGACTTTTTAAAGCATCAATTGCTGGCTGCACAAGATGTATCACCCTGTCGGTTCCTGCCTGGGTTTTCGGCAGGGTG  

I would like to compare each sequence from file1 vs file2 (without considering the headers) (File1: 1,2,3 vs File2: 1,2,3) and if exactly 30 continuous characters were identical from both sequences, I would like to save the headers of the sequences having a match (only those with a match) in an outputfile.

For example the 30 characters string:

TGATTGCATAACCACTTAACATCTTGTTTT is present in seq1 from file1 and seq1 from file2.
TCCATATACTGTATAGACCATCTGAGCGTT is present in seq1 from file1 and seq2 from file2.

so I will end up with an output file like:

[1,1]  
[1,2]  
...
3

1 Answer 1

4
$ cat tst.awk
BEGIN { wid = 30 }
sub(/^>/,"") { hdr=$1; next }
NR == FNR { a[hdr]=$0; next }
{
    for ( hdrA in a ) {
        strA  = a[hdrA]
        lgthA = length(strA)
        for ( idxA=1; idxA<=(lgthA - wid + 1); idxA++ ) {
            substrA = substr(strA,idxA,wid)
            if ( index($0, substrA) ) {
                printf "[%s,%s]\n", hdrA, hdr
                break
            }
        }
    }
}

$ awk -f tst.awk file1 file2
[1,1]
[1,2]
10
  • this was so helpfull but would this also work for big files? lets say 2 files with 10,000 or 50,000 sequences each? Thank you again.
    – GSQ
    Dec 16, 2020 at 16:32
  • @GSQ yes, that shouldn't be a problem.
    – Ed Morton
    Dec 16, 2020 at 17:19
  • Dear @EdMorton I was "playing" with the script you kindly provided. However, I have a question that maybe you could help me to answer. When I run your script ``` awk -f file1 file2 > output_file ``` for some reason it does not recognize a match between this two sequences (next comment) but I know there is match "CACGACGGACGTTAGCACCCGCCGTGTGTC" I already check the input files and they seem to be ok. Would it be possible that when running a lot of sequences, there could be false negatives in the output? In this case I ran 10000 sequences from file1 and 10000 from file2 Thank you!
    – GSQ
    Dec 28, 2020 at 0:26
  • (This comes from file1) >8603 TAGCCCCGTTACATCTTCCGCGCAGGCCGACTCGACCAGTGAGCTATTACGCTTTCTTTAAATGATGGCTGCTTCTAAGCCAACATCCTGGCTGTCTGAGCCTTCCCACATCGTTTCCCACTTAGTAGTAATTTGGGACCTTAGCTGGCGGTCTGGGTTGTTTCCCTCTCCACGACGGACGTTAGCACCCGCCGTGTGTC (This comes from file2) >9977 CTTGGGGACCTTAGCTGTCGGTCAGGGTTGTTTCCCTCTTCACGACGGACGTTAGCACCCGCCGTGTGTCTGCCGACTAGTACTCCCGGGTATTCGGAGTTTGGTTAGGTTTGGTAAATCGGTGAGATCCCCTAGCCCATCCAGTGCTCTACCCCCCGGGGTATTCGGTCGACGCTCTACCTAAATAGATTTCGCGGAGA
    – GSQ
    Dec 28, 2020 at 0:27
  • 1
    Thank you @Ed Morton! Now its working perfectly! :)
    – GSQ
    Dec 29, 2020 at 4:15

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .