Hey,
I would like to map reads from my metatranscriptome to contigs of my reference metagenome. I used bowtie2 with the commands:
./bowtie2 -x [reference file] -f [sequence file] -S [SAM output] --local
I ended up with a 2GB .sam file, which I converted to .bam using samtools. Unfortunately, samtools didnt manage to sort the .bam on its own so I used sorter (https://github.com/brendanofallon/StreamSorter). After creating an index I used the command samtools flagstat [both sorted and unsorted .bam file] for the following output:
6270380 + 0 in total (QC-passed reads + QC-failed reads)0 + 0 duplicates
4564799 + 0 mapped (72.80%:nan%)
I indexed the .bam file and used samtools idxstats [sorted or unsorted .bam file] > [output file]
I end up with a tab separated file containing the sequence identifier, some number followed by 0 and 0
eg. Ref_gDNA_contig_397172|m.2379519 1363 0 0
According to the manual the second tab refers to the length, the third to the number of mapped reads and the fourth to the number of unmapped reads. The length cant be true since my original contig was way shorter. And since no sequences appear to be mapped this doesn't match the output of 'flagstat'. Furthermore, I know that the 'flagstat' command gives the right output since I crosschecked with CLC Genomics Workbench mapping (same mapping fraction).
Anyone who can help me here?
You're better off using samtools sort to sort your BAM files, so something like
samtools view -Sbu foo.bam | samtools sort -@ 4 - output_prefix. The likelihood of samtools idxstats getting something like that wrong if the file was properly sorted is quite low, so just rerun things after sorting with samtools.Also, the contig lengths are reported by bowtie2, so in the unlikely even that that value is wrong then bowtie2 is to blame.