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 is 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 numer 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 doesnt 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?