These are the results for
samtools idxstats input.sorted.bam
X1 3580 9483 0
X2 2836 8110 0
X3 2188 8220 0
X4 1468 6698 0
* 0 0 0
These are the results of samtools flagstat
samtools flagstat input.sorted.bam
32511 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
32511 + 0 mapped (100.00%:-nan%)
32511 + 0 paired in sequencing
16715 + 0 read1
15796 + 0 read2
31500 + 0 properly paired (96.89%:-nan%)
31522 + 0 with itself and mate mapped
989 + 0 singletons (3.04%:-nan%)
10 + 0 with mate mapped to a different chr
6 + 0 with mate mapped to a different chr (mapQ>=5)
How can I calculate mapping rate say for X1 from this data? Will it be
(number of reads mapped to X1/ Total mapped reads)*100
(9483/32511)*100 = 29.17%
Please let me know if this is correct.
However you calculate it, your mapping rate is really low, indicating problems with your data. Is this just an example, or are these numbers from your data?
Thank you for replying. These are the numbers from my data. The multi fasts file reference which I used is a synthetic construct. The total mRNA data was mapped to this reference. The data is Miseq with PE 300. What can I do to increase the mapping rate.? I used bowtie2 with -I 0 and -X 1000 with -fr option.
Could you add further details? Is this a complete genome, with a synthetic construct added? Or your reference genome consists of a synthetic construct only? What is the organism you are studying?
The reference is sequence of 4 transcripts. I performed a BLAST and came to know it mapped closely to immunoglobulin light and heavy chain. This is all the information which I could gather.
You can't map a RNAseq dataset to a reference of just four transcripts, there probably will be a lot of spurious mappings due to relatively similar transcripts mapping to your "reference", as the correct genes are not present in your reference.
Why are you doing this? What exactly do you expect to accomplish?
I wish to know variants in my data with respect to this reference. But even though the alignment rate is low I was interested to know if mapping rates for each sequence in a multi-fasta reference file can be calculated using SAM or BAM file?
You first need to address your very low mapping rate as addressed by d-cameron
You also should align to the full genome using a spliced aligner such as HISAT2 or STAR.
Thank you for the suggestion. I will use these tools and align input data to a complete genome.
With a <1% mapping rate you have a serious issue with your data. Try the following:
Run fastqc on your input files. Are there any errors or warning?
Align to the reference genome for your species with default bowtie2 setting. Are you now getting most of your reads mapped? If so, your synthetic construct is likely to be garbage. Work out what when wrong and make a new one.
If you still have 99% of your reads not mapping to your reference genome then it's likely to be an issue with library preparation/sequencing. The most likely explanation is that your library preparation failed (e.g did you accidentally wash all your DNA away in one of your steps?) and there's nothing you can do bioinformatically to recover.
Try BLASTing a few hundred of the unmapped reads. Do they actually come from your species?
When you say overall alignment rate (OVR), you need to consider a lot of factors
All of these will contribute to the overall alignment rate which I have tried explaining here
Calculating OVR for all the sequences individually will be losing a lot of that information (for example, exon skipping, probable duplications, translocations etc.). If you want to know the just the number of raw reads, than that's a separate stuff.