Others may have command-line solutions. Here's one in python, untested, that I think should do what you need.
Fasta files:
one.fasta
>miRNA1
AGTCGTCA
>miRNA2
AGTTCGTT
>miRNA3
AGTCGTGA
two.fasta
>miRNA1
TTGACT
>miRNA2
AGTCGTCA
three.fasta
>miRNA3
TTACATT
>miRNA5
TTGACT
>miRNA3
AGTCGTGA
First, make a file with all the unique sequences from each of your fasta files:
cat file1.fasta file2.fasta file3.fasta | grep -v '>' | sort | uniq > all.txt
Save this code as comp_seqs.py , and run as python comp_seqs.py > output.xls in the location that has the three fasta files, and the new file with all the sequences:
#!/usr/bin/env python
import regex
from collections import defaultdict
from itertools import izip
def open_fasta(file):
""" Open fasta file, and store in dictionary, as {'sequence:header'} """
s = {}
with open(file, 'r') as fasta:
for line in fasta:
if line.startswith('>'):
h = line.strip().split('>')[1]
n = next(fasta).strip()
s[n] = h
return s
def hamming_distance(s1, s2):
""" Count number of mismatched characters in equal length strings. """
if len(s1) != len(s2): raise ValueError('string lengths do not match')
return sum(a != b for a, b in izip(s1, s2))
def find_seqs(fh, p):
""" Find sequences from 'all sequence list' that either match as identical sequence, or mismatch of 1, store as defaultdict [{sequence:[id1, id2, id3...]}] """
seqs = open_fasta(fh)
for i in seqs:
if i in x:
match[i][p] = seqs[i]
else:
for j in x:
if len(i) == len(j):
if hamming_distance(i, j) <= 1:
match[i][p] = seqs[i] + '+1'
# Open concated fasta file, containing unique sequences
with open('all.txt', 'r') as a:
x = [line.strip() for line in a]
# empty defaultdict for reulsts
match = defaultdict(list)
for i in x:
for n in range(3):
match[i].append('')
# Open each fasta, and search for sequences in 'all sequence list'
find_seqs('one.fasta', 0)
find_seqs('two.fasta', 1)
find_seqs('three.fasta', 2)
# Print results in tab delimited format
print 'Sequence\tFile1_ID\tFile2_ID\tFile3_ID'
for m in match:
print m, '\t', '\t'.join(match[m])
Output from example above:
Sequence File1_ID File2_ID File3_ID
AGTCGTCA miRNA1 miRNA2
TTGACT miRNA1 miRNA5
AGTTCGTT miRNA2
AGTCGTGA miRNA3 miRNA3
TTACATT miRNA3
I usually like blastclust which is now a legacy software but this post recommending cd-hit might be useful How to cluster sequences present in multiple fasta files
Thx for the help. The cd-hit webserver seems to take forever for a small job and can only handle two files. So I dont think this will solve my problem
Are these long sequences? How many sequences per file? Are they in single-line format, or multi-line?
The sequences are short <30bp and the files are in single line format.
When applying the 1 mismatch between sequences, does this count for length too?
Yes, this would be perfect.
I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below: