Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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, 02:23 AM.

Latest Articles

Collapse

  • seqadmin
    Exploring the Dynamics of the Tumor Microenvironment
    by seqadmin




    The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
    07-08-2024, 03:19 PM
  • seqadmin
    Exploring Human Diversity Through Large-Scale Omics
    by seqadmin


    In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
    06-25-2024, 06:43 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Yesterday, 11:09 AM
0 responses
18 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-19-2024, 07:20 AM
0 responses
148 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-16-2024, 05:49 AM
0 responses
124 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-15-2024, 06:53 AM
0 responses
111 views
0 likes
Last Post seqadmin  
Working...
X