Hello.
I have a miRBase database of mature miRNA sequences that I want to allign against my adapter-stripped small RNA fastq files.
I am using the Rsubread package in R, but I am getting 0% mapping on every run I attempt.
I have no problems when building my reference DB using the 'build' function like so:
buildindex(basename=file.path("filepath","miRBaseSeqs"),reference="miRBaseSeqs.fa")
However, when I try to run the 'align' function:
align("miRBaseSeqs", "AdapterStripped.fastq", input_format="FASTQ", output_file="testAlgn.SAM", indels=1, minFragLength=15, maxFragLength=55)
It runs fine with no errors:
//========================== subread-align setting ===========================\\
|| ||
|| Function : Read alignment ||
|| Threads : 1 ||
|| Input file : AdapterStripped.fastq ||
|| Output file : testAlgn.SAM (SAM) ||
|| Index name : miRBaseSeqs ||
|| Phred offset : 33 ||
|| ||
|| Min votes : 3 ||
|| Max indels : 1 ||
|| # of Best mapping : 1 ||
|| Unique mapping : yes ||
|| Hamming distance : yes ||
|| Quality scores : no ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
//====================== Running (13-Mar-2014 11:33:58) ======================\\
|| ||
|| The input file contains base space reads. ||
|| Load the 1-th index block... ||
|| Map reads... ||
|| 0% completed, 0 mins elapsed, total=12818k reads, rate=69.8k/s ||
|| 6% completed, 0 mins elapsed, total=12802k reads, rate=74.7k/s ||
|| 12% completed, 0 mins elapsed, total=12818k reads, rate=75.4k/s ||
... etc
But the output is always this:
//================================= Summary ==================================\\
|| ||
|| Processed : 12826628 reads ||
|| Mapped : 0 reads (0.0%) ||
|| Indels : 0 ||
|| ||
|| Running time : 5.9 minutes ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
I have tried multiple argument parameters such as:
align("miRBaseSeqs", "AdapterStripped.fastq", input_format="FASTQ", output_file="testAlgn.SAM", indels=1, minFragLength=15, maxFragLength=55, DP_GapOpenPenalty=10, DP_GapExtPenalty=8, DP_MismatchPenalty=-5, DP_MatchScore=10)
subjunc("miRBaseSeqs", "PS24h.AdapterStripped.fastq", input_format="FASTQ", output_file="testAlgn.SAM", indels=1, minFragLength=15, maxFragLength=55, DNAseq=TRUE, reportAllJunctions=TRUE)
and variation on the theme, but I always get that 0% aligned output.
Does anyone have any suggestions? Anything is appreciated.
Thanks in advance!
•
written
by
Tom A •
1 answer
Hi Tom, thank for using our Subread aligner. I noticed that you required subread not to report multi-mapping reads in the SAM file ("Unique mapping : yes"). If the reference genome is highly repetitive, you will get zero mapped read.
You can add a parameter, unique=FALSE, to the align function. This will make subread to report all mapping locations; if a read can be mapped to more than one location, a random one of these locations is reported. Or, if you want all these locations, just set the number of locations you want by using the nBestLocations parameter.
If it still does not work, please contact us on liao (at) wehi.edu.au. Thank you.
•
written
by
lavertonvic •
Log in to answer this question.
-
miRNA alignment with Bowtie2
written by Chironex •
-
ATAC seq analysis
written by aranyak111 •
-
Problem with generating BAM files with STAR
written by gokberk
-
Rsubread: ERROR: cannot finish the SAM file
written by TheCheshireCat •
-
Rsubread subjunc :Missing entry : 1733000143 for 0
written by inowkeren •
-
Interpreting Flagstat Output
written by roll
-
Analyzing Older Illumina Paired-End Data
written by danrdanny •
-
Fastq File Without Index?
written by lsvijfhuizen •
-
Bowtie Inaccuracy Limit
written by Arpssss •
-
The Meaning Of B In Illumina 1.5 Pipeline Data?
written by Martin A Hansen
Just to check, try mapping with bowtie1 and see do you still not get any reads mapped. You can use -n 2 -a -l 25 as options. Let me know
So I have used OSA to map a similar file; one that also contained miRNA seqs, but didn't have an option for all the mature seqs (3p, 5p arm) hence the need to use a different aligner (one that enables you to build your own reference library). Anyway, long story short, the adapter striped files are fine.
The problem is that Subread does not report those reads that have bases which were mapped outside of reference sequences. You can get around this by padding your miRNA sequences with N's and then build an index for them and do the alignments. For example, you may add 100 N's to the left side of each of your miRNA sequences and also add 100 N's to the right side. This will prevent 100bp reads from overhanging your reference sequences, and the mapped reads should then be reported by Subread.