SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
gap penalty in BWA Victory Bioinformatics 1 01-14-2021 06:19 AM
How to Set Clipping Penalty in BWA ALN? svos Bioinformatics 4 03-17-2015 11:36 PM
mapping reads using BWA, no MM, no GAP nejma Bioinformatics 1 07-01-2014 02:50 AM
BWA change mismatch penalty between N-mismatch and base-mismatch andreac Bioinformatics 0 12-19-2012 08:05 AM

Reply
 
Thread Tools
Old 01-14-2021, 11:07 AM   #1
annguyen
Junior Member
 
Location: Cambridge

Join Date: Nov 2020
Posts: 3
Default High gap penalty in BWA

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 goes up while this number should go down because of gap penalties.

Here are the result (I combined ref1 and ref2 in one file called refx.fasta):

```
bwa mem -O 20 -E 20 refx.fasta R1.fastq.gz R2.fastq.gz > mysam.sam

samtools flagstat mysam.sam

1117216 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
61620 + 0 supplementary
0 + 0 duplicates
1116527 + 0 mapped (99.94% : N/A)
1055596 + 0 paired in sequencing
527798 + 0 read1
527798 + 0 read2
1053928 + 0 properly paired (99.84% : N/A)
1054666 + 0 with itself and mate mapped
241 + 0 singletons (0.02% : N/A)
6 + 0 with mate mapped to a different chr
5 + 0 with mate mapped to a different chr (mapQ>=5)```

```
bwa mem -O 30 -E 30 refx.fasta R1.fastq.gz R2.fastq.gz > mysam.sam

samtools flagstat mysam.sam

1121698 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
66102 + 0 supplementary
0 + 0 duplicates
1121007 + 0 mapped (99.94% : N/A)
1055596 + 0 paired in sequencing
527798 + 0 read1
527798 + 0 read2
1053912 + 0 properly paired (99.84% : N/A)
1054662 + 0 with itself and mate mapped
243 + 0 singletons (0.02% : N/A)
6 + 0 with mate mapped to a different chr
5 + 0 with mate mapped to a different chr (mapQ>=5)
```

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

Last edited by annguyen; 01-15-2021 at 01:23 AM.
annguyen is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:35 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO