I'm using BWA MEM to align paired-end fastq files to two reference sequences (50bp and 80bp):
>ref1
TCGTAACGCAAGTTGGATACTCTCGA******************************GGATGTTGCCGTCCTCCTTGAAGT
>ref2
TCGTAACGCAAGTTGGATACTCTCGATTGCAAGTAGTCGATTGCATTGTCAATCTAGGATGTTGCCGTCCTCCTTGAAGT
These two sequences are idential except the middle bit. I then used samtools to filter out properly paired reads from each sequence based on identifiers but surprisingly it showed overwhelming majority of reads have the number of matches that correspond to **ref2** and 0 properly paired reads for **ref1**
I verified this result by repeating the process for each sequence individually and got the same output:
```
$ samtools flagstat mysam.sam
1055596 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
20960 + 0 mapped (1.99% : N/A)
1055596 + 0 paired in sequencing
527798 + 0 read1
527798 + 0 read2
0 + 0 properly paired (0.00% : N/A)
20934 + 0 with itself and mate mapped
26 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
```
I then tried setting high values for gap opens and extension penalties in order to limit the chance that hits from alignments of **ref1** reads result in origin of **ref2**. By doing this, I now can see some properly paired reads from **ref1**. However, one thing I don't understand is that the total number of reads as well as properly paired reads increased while these should be reduced because of gap penalties.
Can someone explain to me:
1. Is it normal to see 0 properly paired reads in this case?
2. Why does disallowing reads with open and extension gaps increase the total of reads?
Many thanks
>ref1
TCGTAACGCAAGTTGGATACTCTCGA******************************GGATGTTGCCGTCCTCCTTGAAGT
>ref2
TCGTAACGCAAGTTGGATACTCTCGATTGCAAGTAGTCGATTGCATTGTCAATCTAGGATGTTGCCGTCCTCCTTGAAGT
These two sequences are idential except the middle bit. I then used samtools to filter out properly paired reads from each sequence based on identifiers but surprisingly it showed overwhelming majority of reads have the number of matches that correspond to **ref2** and 0 properly paired reads for **ref1**
I verified this result by repeating the process for each sequence individually and got the same output:
```
$ samtools flagstat mysam.sam
1055596 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
20960 + 0 mapped (1.99% : N/A)
1055596 + 0 paired in sequencing
527798 + 0 read1
527798 + 0 read2
0 + 0 properly paired (0.00% : N/A)
20934 + 0 with itself and mate mapped
26 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
```
I then tried setting high values for gap opens and extension penalties in order to limit the chance that hits from alignments of **ref1** reads result in origin of **ref2**. By doing this, I now can see some properly paired reads from **ref1**. However, one thing I don't understand is that the total number of reads as well as properly paired reads increased while these should be reduced because of gap penalties.
Can someone explain to me:
1. Is it normal to see 0 properly paired reads in this case?
2. Why does disallowing reads with open and extension gaps increase the total of reads?
Many thanks