Hi all, I have two sets of data (let's say A & B ) each includes 1000,000 sequences. I would like to see how many of sequences of one file (e.g. A) has at lease 90% matches to the sequences in the other file (e.g. B). I have already wrote a script which take one sequence from A and do pairwise alignment to all sequences in B and store the score for each matches. (I have used Biopython module) My problem is, the script is very slow and still it did not finish after 4 days run on 21 cores. Do you guys have any idea that how can I improve the script? or is there any better way to do this massive pariwise alignment?
Thanks
3 answers
I am answering my question, since It may help people with similar problem. The fast massive alignment is possible by applying Bowie. When you do alignment using Bowtie2 it will tell you how many sequences are repeated. And you can delete the repetitive sequences by extracting it from sam file (Bowtie output)
"massive pairwise alignment" ? makes me think of blast ! :-)
I thus suggest to do a blast search and parse the resulting output file.
Yes it might be, you could benefit from tweaking the blast parameters a little (eg. do ungapped blast).
I do wonder however how different the result are (can be?) between these two approaches, despite that I'm very well aware from the fact blast is not actually an alignment tool. It will nonetheless be able to answer your question if that is to look for >90% identity.
Are we talking DNA or protein here btw?
Log in to answer this question.
-
Cut down CD-HIT-EST time
written by Neeraja •
-
Pairwise alignment of multi-FASTA file sequences
written by aurora •
-
Protein semi-global alignment tools? Blast for glocal alignments?
written by chaperoneline •
-
Get pairwise scores from alignment without re-aligning?
written by mpd •
-
How to calculate all vs all sequence divergence from an alignment?
written by ThulasiS
-
How important is --frag-bias-correct in running Cuffdiff (it makes Cuffdiff very slow in our case)?
written by tunl
-
Efficient implementation of Fuzzy set operations on sets of sequences?
written by m.aljaff90 •
-
crossover operator for evolutionary algorithm
written by Moses
-
MSA building algorithm
written by ga29qal •
-
Pairwise Sequence Alignment With R
written by Lakshmi •
Instead of re-inventing the wheel use one of the known MSA programs (clustal, muscle, t-coffee, mafftt). How long are these sequences? A MSA of a million sequences may be tough to view/decipher. Can you de-dupe your data before doing MSA?
Those softwares are deigned for multiple sequence alignment. right? But I dont want to do MSA. seq size ranges btw 100-300 bp.
My bad. Perhaps using an NGS aligner would still be much faster. I suggest looking at bbmap.sh from BBMap suite. Let me try a small test.
Actually these sequences are not related to RNA-seq. I am trying to find unique sequences among a dataset with a million sequence then create a distance matrix based on that.
That is ok. NGS aligners are going to be inherently fast compared to other alignment programs, as long as they can handle fasta format sequences (which bbmap can).
This part is very easily done with
clumpify.sh. It will be insanely fast. You can allow for errors (subs=option). See this thread for usage information: A: Introducing Clumpify: Create 30% Smaller, Faster Gzipped Fastq Files You can useclumpifywith fasta sequences. Where ever the commands sayfq.gzyou can replace those with.fafor your files.that I am not sure how to do right away.
I'm also puzzled by the 'distance matric' part. Is it correct to assume you mean 'distance matrix' as in phylogenetic context?
You want to create it based on the unique sequences in your dataset? What do you consider ''unique' , not otherwise present in the dataset or less than 90% identical/similar?