I am doing data testing on haplotype. I generate two genome: a raw hg38 genome and a mask hg38 genome after mask haplotype SNP position.
example region sequence (Italic base indicates SNP):
echo -e "chr1\t107904139\t107904229" | bedtools getfasta -fi hg38.fa -bed -
chr1:107904139-107904229
TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATT_GA_GTTAAAATG
echo -e "chr1\t107904139\t107904229" | bedtools getfasta -fi hg38.maskallele.fa -bed -
index file genome.maskallele.fa.fai not found, generating...
chr1:107904139-107904229
TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATT_NN_GTTAAAATG
example for reads alignment to raw genome, it is correct:
SRR1609148.25 83 chr1 107904140 60 90M = 107903755 -475 TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATTATGTTAAAATG 3DDDDDECCEEEFDECHHHFJIGGGAIICHFEJJIIIJIJIHIJJJJJJIJJIIHGFCJJJJIIHHFEHHGJIIJJJGHHHHDFFFFCCC NM:i:2 MD:Z:79G0A9 AS:i:80 XS:i:19
example for reads alignment to mask genome, it may have error for MD tag (the MD tag is 79TG9, not 79NN9):
SRR1609148.25 83 chr1 107904140 60 90M = 107903755 -475 TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATTATGTTAAAATG 3DDDDDECCEEEFDECHHHFJIGGGAIICHFEJJIIIJIJIHIJJJJJJIJJIIHGFCJJJJIIHHFEHHGJIIJJJGHHHHDFFFFCCC NM:i:2 MD:Z:79T0G9 AS:i:80 XS:i:19
I am doing data testing on haplotype. I generate two genome: a raw hg38 genome and a mask hg38 genome after mask haplotype SNP position.
example region sequence (Italic base indicates SNP): echo -e "chr1\t107904139\t107904229" | bedtools getfasta -fi hg38.fa -bed -
example for reads alignment to raw genome, it is correct: SRR1609148.25 83 chr1 107904140 60 90M = 107903755 -475 TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATTATGTTAAAATG 3DDDDDECCEEEFDECHHHFJIGGGAIICHFEJJIIIJIJIHIJJJJJJIJJIIHGFCJJJJIIHHFEHHGJIIJJJGHHHHDFFFFCCC NM:i:2 MD:Z:79G0A9 AS:i:80 XS:i:19 example for reads alignment to mask genome, it may have error for MD tag (the MD tag is 79TG9, not 79NN9): SRR1609148.25 83 chr1 107904140 60 90M = 107903755 -475 TCAGGGCCACTACGTAGACACTGGCCAAATGACCCTTTGGTATCAATCAGAAAAGGCACCCTTCCTTCCTTCCTCAATTATGTTAAAATG 3DDDDDECCEEEFDECHHHFJIGGGAIICHFEJJIIIJIJIHIJJJJJJIJJIIHGFCJJJJIIHHFEHHGJIIJJJGHHHHDFFFFCCC NM:i:2 MD:Z:79T0G9 AS:i:80 XS:i:19
I don't know if bwa is intentional or a bug.