Entering edit mode
21 months ago
EarthShaker-EAC
•
0
I used a genome A and an A+B genome to construct two A.db and AB.db with bwa respectively. The reads can be alignment with A alone, but only the B genome is alignment in the results of AB. I used the -a parameter of bwa mem during the alignment, but why not all the alignment results are output in the result.
Input_A:
bwa mem -a A.ref.fna query.fa |samtools view -F 4 -bS|grep "NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76"
OutPut_A:
NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76 0 ref_A 1569869 60 76M * 0 0 ACTTAATTGTGTTAGAAGCTCTGAATTTTCAATCCAACCTAATACAGCTGGACCGATAACGATACCGACAATTAAT * NM:i:1 MD:Z:9A66 AS:i:71 XS:i:0
Input_AB:
bwa mem -a AB.ref.fna query.fa |samtools view -F 4 -bS |grep "NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76"
OutPut_AB:
NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76 0 ref_B 119652595 0 76M * 0 0 ACTTAATTGTGTTAGAAGCTCTGAATTTTCAATCCAACCTAATACAGCTGGACCGATAACGATACCGACAATTAAT * NM:i:0 MD:Z:76 AS:i:76 XS:i:76
NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76 256 ref_B 33868779 0 76M * 0 0 * * NM:i:0 MD:Z:76 AS:i:76
NDX550718_RUO:5:HF5H7BGXK:1:11101:8944:11737|len=76 256 ref_B 96156965 0 76M * 0 0 * * NM:i:0 MD:Z:76 AS:i:76
genome size(bp):
ref_A 5427183
ref_B 273598709
I suggest you look into the NM:i values. In ref_A is 1, in ref_B is 0. By looking at older posts it looks like the value is the edit distance between the sequence and the reference. So apparently you map exactly to ref_B and with a SNP to ref_A. This is my understanding, you may want to investigate in more detail.