count reads mapped to multiple locations only once

Question

[snishtala03]

Hello,

I have paired end data which I am mapping (using bwa - mem) to 17 viral references which are ~85% similar to one another (i.e I have one reference file that I give to bwa mem which has 17 different fastqs). So, obviously, I have most of my reads mapping to multiple references. Which is needed for my analysis.

However, I also want to get a count of uniquely mapped reads, i.e, if a read maps to 6/17 references, I want to count it only once. I have tried to filer my bam file based on quality or using grep on XA and XS tags but this removes all the 6 mapping.

samtools view -h sorted.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique.bam

samtools view -h sorted.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique.bam

Is there a way I can keep one instance of those 6 mappings?

[GenoMax]

Not answering your question directly but providing a suggestion to get what you need.

If you align with bbmap.sh you get the option of treating multi-mapping reads by placing them in only one of those 6 locations by using ambig=random option.

[GenoMax]

You could also grab the alignments via

samtools view -h sorted.bam | grep  -e 'XA:Z:' -e 'SA:Z:'

make then non-redundant (leave one copy you want) and then add them back to your unique.bam.

[snishtala03]

Thanks for your response!

So while I was trying to this, I came across this example below. It shows that a read has secondary alignments and supplementary alignments, but not all of them are seen in the bam file. Is this how secondary and supplementary alignments usually appear or does this mean that secondary ans supplementary alignments are marked but they are not shown as outputs?

M02091:48:000000000-CBW5K:1:2105:5397:6372      81      HBV_genotype_A_X02763   3271    0       16S112M23S      HBVpcRNA1_Genotype_C_3405   11      0       GCTCTTCCGAGAGAGGGGTCTGCGAAAAAGCACCATGCAACTTTTTCAAAGCTGCCTAAGAATATCTTGTACATGTACAACTGTTCAAGAAGACAAGCTGTGAATTGGGTGGCTTTGGGGCATGGACATTGACCCTTATAAAGAATTTGGA     HHHGGGGGFGF?FCGHHHGGGGGFHHGFHGFDHGHHHHHHCHHHHHEHFHGGHHHHHHGGDGE3GHHHHHHHHHGF5GEHEGHHFHHHGHHGGGHGHHHHGFF2HGGHHGHHHHGGGGHHHHHHHHHHHGGGGGGDGGGFFFFFFFBBBBB     NM:i:17 MD:Z:8C1C0C20C0C0T8T0C2C12C1C10C0C0T0C9C0C24        AS:i:78 XS:i:76 SA:Z:HBV_genotype_G_AB064313,3403,-,112S39M,0,4;
M02091:48:000000000-CBW5K:1:2105:5397:6372      2129    HBV_genotype_G_AB064313 3403    0       112H39M HBVpcRNA1_Genotype_C_3405   11      0       CTTTGGGGCATGGACATTGACCCTTATAAAGAATTTGGA HHGGGGHHHHHHHHHHHGGGGGGDGGGFFFFFFFBBBBB NM:i:4  MD:Z:4T0T3T1A27     AS:i:31 XS:i:31 SA:Z:HBV_genotype_A_X02763,3271,-,16S112M23S,0,17;      XA:Z:HBV_genotype_G_AB064313,-146,112S39M,4;HBV_genotype_G_AF160501,-146,112S39M,4;HBV_genotype_G_AF160501,-3403,112S39M,4;HBV_genotype_F_X75658,-159,125S26M,0;
M02091:48:000000000-CBW5K:1:2105:5397:6372      161     HBVpcRNA1_Genotype_C_3405       11      0       114M24S HBV_genotype_A_X02763       3271    0       ATTGGTCTGCGCACCAGCCCCCTTCAACTTTTTCCCCTCTGCCTAATCCTCTCTTGTACCTTTCCCACTGTTCCCTCCTCCAAGCTGTTCCTTTGGTGGCTTTTGGGCATGGACATTGACCCTTATAAAGAATTTGGA  >>>AAFFFFFBBGGGGGGFGFGGCGHHHHHHHHHHHHHGFFHHHHGHHHHHFHHHHHHHHH5FHHGHHHHHHHHFGFFGHHHEHGHHFDFGGG@FFGFAHHFH3FEGHHHHHHHHHHHHHHHHGGHGFFHHHHHHHHF  NM:i:17 MD:Z:9T0T7A2A1G24A4A3T1A1G3T7A0A0G12G4G9G10     AS:i:80     XS:i:80 SA:Z:HBV_genotype_G_AF160501,146,+,99S39M,0,3;  XA:Z:HBVpcRNA1_Genotype_C_3405,+3268,114M24S,17;
M02091:48:000000000-CBW5K:1:2105:5397:6372      2209    HBV_genotype_G_AF160501 146     0       99H39M  HBV_genotype_A_X02763       3271    0       CTTTTGGGCATGGACATTGACCCTTATAAAGAATTTGGA HHFH3FEGHHHHHHHHHHHHHHHHGGHGFFHHHHHHHHF NM:i:3  MD:Z:5T3T1A27       AS:i:33 XS:i:33 SA:Z:HBVpcRNA1_Genotype_C_3405,11,+,114M24S,0,17;       XA:Z:HBV_genotype_G_AF160501,+3403,99S39M,3;HBV_genotype_G_AB064313,+3403,99S39M,3;HBV_genotype_G_AB064313,+146,99S39M,3;

To get secondary and supplementary reads in different lines, I also tried using -a and -M flags when running bwa-mem, but then these secondary and supplementary reads appear as * . A little confused here, am I missing something? Showing a few examples below -

M02091:48:000000000-CBW5K:1:2105:5397:6372      81      HBV_genotype_A_X02763   3271    0       16S112M23S      HBVpcRNA1_Genotype_C_3405       11      0       GCTCTTCCGAGAGAGGGGTCTGCGAAAAAGCACCATGCAACTTTTTCAAAGCTGCCTAAGAATATCTTGTACATGTACAACTGTTCAAGAAGACAAGCTGTGAATTGGGTGGCTTTGGGGCATGGACATTGACCCTTATAAAGAATTTGGA HHHGGGGGFGF?FCGHHHGGGGGFHHGFHGFDHGHHHHHHCHHHHHEHFHGGHHHHHHGGDGE3GHHHHHHHHHGF5GEHEGHHFHHHGHHGGGHGHHHHGFF2HGGHHGHHHHGGGGHHHHHHHHHHHGGGGGGDGGGFFFFFFFBBBBB NM:i:17 MD:Z:8C1C0C20C0C0T8T0C2C12C1C10C0C0T0C9C0C24    AS:i:78 XS:i:76 SA:Z:HBV_genotype_G_AB064313,3403,-,112S39M,0,4;
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_genotype_D_AF121240 14      0       16H112M23H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:18 MD:Z:8C1C0C20C0C0T8T0C2C6T5C1T10C0C0T0C9C0C24   AS:i:76
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_Genotype_A_EuroHEP_AJ012207 14      0       16H110M25H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:17 MD:Z:8C1C0C20C0C0T8T0C2C12C1C10C0C0T0C9C0C22    AS:i:76
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_Genotype_A_EuroHEP_AJ012207 3271    0       16H110M25H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:17 MD:Z:8C1C0C20C0C0T8T0C2C12C1C10C0C0T0C9C0C22    AS:i:76
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_genotype_A_X02763   14      0       16H110M25H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:17 MD:Z:8C1C0C20C0C0T8T0C2C12C1C10C0C0T0C9C0C22    AS:i:76
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_genotype_E_AB032431 14      0       16H111M24H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:18 MD:Z:8C1C0C20C0C0T8T0C2C6T5C1T10C0C0T0C9C0C23   AS:i:75
M02091:48:000000000-CBW5K:1:2105:5397:6372      337     HBV_genotype_D_AF121240 3271    0       16H117M18H      HBVpcRNA1_Genotype_C_3405       11      0       *       *       NM:i:21 MD:Z:8C1C0C20C0C0T8T0C2C6T5C1T10C0C0T0C9C0C22A0C2A2     AS:i:75