Dear Biostars,
I'm trying to extract multi-mapped and unique reads from a SAM/BAM file after doing an alignment using bowtie aligner with single end reads.
./bowtie -p 6 -v 0 -a -f --sam -- best '/hg38' '/myfasta.fa' '/results.sam'
results.sam:
URS00001E85BB_9606 0 chr19 37279510 255 26M * 0 0 TCCAATGCAGTGAGGCCGAAAGGGTC IIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:26 NM:i:0 XM:i:15
URS00001E85BB_9606 0 chr19 37292136 255 26M * 0 0 TCCAATGCAGTGAGGCCGAAAGGGTC IIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:26 NM:i:0 XM:i:15
URS00001E85BB_9606 0 chr19 37297187 255 26M * 0 0 TCCAATGCAGTGAGGCCGAAAGGGTC IIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:26 NM:i:0
...
URS00003A9162_9606 4 * 0 0 * * 0 0 TAATGATCTGCCTCAAACCTAGTTTTCNGCC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII XM:i:0
URS00005E6300_9606 0 chrX 99575628 255 28M * 0 0 TGCAGTAAATTGGTATTGAGGTTGTGGG IIIIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:28 NM:i:0 XM:i:2
URS00003E8A0A_9606 16 chr18 405642 255 28M * 0 0 AGGTATTTCAATCACGCTGATCCATTCA IIIIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:28 NM:i:0 XM:i:2
URS00002B3810_9606 4 * 0 0 * * 0 0 TGCTTGTGAACTGCTTGTGAATGGGTGTGT IIIIIIIIIIIIIIIIIIIIIIIIIIIIII XM:i:0
URS00000E8FF1_9606 16 chr3 12851709 255 29M * 0 0 TCACCCTAGAGCTTTGCACACTGTTACCA IIIIIIIIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:29 NM:i:0 XM:i:2
After that, I convert it to bam, then extracted only mapped reads:
samtools view -u -F 4 'results.bam' > mappedResults.bam
and in order to get only unique mapped reads I did:
samtools view mappedResults.bam | awk '{print 10}' | awk 'NF > 0' | sort | uniq -u > uniqResultsSeq.txt
"Extract with awk each sequence, take into account only non-zero field lines, sort them and extract only unique sequences."
If I want to extract unique and multi-mapped reads but "collapsed", then I used:
samtools view mappedResults.bam | awk '{print 10}' | awk 'NF > 0' | sort -u > uniqResultsSeq.txt
Still, I haven't found a way to extract only the multi-mapped reads. Maybe I should use XM> 1? In this answer they use XS flag to find the uniquely mapped reads but it was used on bowtie2. Is there another straightforward way to get unique/muti-mapped reads using samtools from bowtie results?
Thank you for your time
Konstantinos
