I am trying to extract alignment scores after running a bunch of fastq Illumina reads against the human pre-miRNAs data set. (converted to DNA U-to-T)
I have ran SSAHA2 with the default parameters . I am interested in this exact output. However, it is not working for some reason.
-> No alignments are returned even though the searches I have done using Biostrings returned a certain amount of matches. ( number of exact matches is relatively low compared to the query database)
I have tried to align the collapsed sequences as well without success.
./ssaha2 hairpins_to_dna.fa NP.fastq >out.txt
Matches For Query 498 (79 bases): hsa-mir-874
===================================================
Score Q_Name S_Name Q_Start Q_End S_Start S_End Direction
===================================================
Does anyone know why is this happening?
It is not working for the hashtable generated from the ref sequences either....
The reads are single end and around 22 nucleotides long. I need to align them to one/more pre-miRNAs and get the alignment scores.
pre-miRNA : TGTCGGGTAGCTTATCAGACTGATGTTGACTGTTGAATCTCATGGCAACACCAGTCGATGGGCTGTCTGACA
reads:
TAGCTTATCAGACTGATGTTGA
TAGCTTATCAGACTGATGTTAT
TAGCTTATCAGACTGATGATC
TAGCTTATCAGACTGATGTATC
The data set contains around 240 000 reads.
•
written
by
Agatha •
1 answer
You got very short query sequences and short database sequences. I would try to set the seed-length to a low value using -kmer option. Try 8, to start with. I would also try BWA, Mosaik, BFAST and others, lower seed length if possible. In principle to get it working you have to try more or less any tool, until you succede.
To get an impression what you are looking for get the FASTA tools http://fasta.bioch.virginia.edu/fasta_www2/fasta_down.shtml and install them locally, extract your query to fasta format. Then run ssearch36 on the first few hundred sequences against your database. This is same algorithm as in R-function but much faster, then look at the output and make a collection of sequences that do align, then with this, procede testing the faster short read aligners.
Log in to answer this question.
-
Speeding Up Psipred (Or Finding A Better Locally-Run Alternative)
written by Dave
-
Bowtie Inaccuracy Limit
written by Arpssss •
-
Isomir Discovery Or Sequencing Contaminants?
written by Agatha •
-
Illumina Rna-Seq Ambiguous Base Calls Identification
written by Agatha •
-
Differing Bwa Mapq Scores For A Read Depending On How I Build The Transcriptome
written by David Quigley
-
Detect Regions In A List That Are Within 1 Megabase Of Each Other
written by Rubal
Do you mean that it is not working because there are not alignments found? Possibly there are no alignments, or none ssaha2 can find with standard settings. Try using matchPattern (if using R) for example on a few sequences from your input to see if you can match some reads.
@Michael Dondrup-I have actually used match pattern and I did find alignments..I am actually running the pairwise alignment function in R..That is all I need actually but it's been taking forever...
@Michael Dondrup -And I am only trying to align it to one pre-miRNA sequence..it has been running for 26 hours now..I need the start and end positions plus the alignment score . I am using one loop and I am calling the pairwise alignment function once..
I meant for every read ..I don;t know what I can do to speed up the process it is not even half way done..
Please provide a little bit more information about how many reads you have, how long they are, are they paired end. Maybe post a few reads which align well, so one could test it. It is clear that you cannot use pairwiseAlign on thousands of reads. ssaha2 is optimized for aligning many short sequences against a few long ones, not against many very short sequences, maybe just try a different tool.
@Michael Dondrup- I have added extra info..what other tool would you recommend? (in order to get start and end positions and alignment scores?)