Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • dpryan
    replied
    No prob, it's a nice break from my own stuff.
    I'm not aware of a samtools flag that would do that, but something like:

    Code:
    samtools view -h foo.bam | awk '{if(substr($1,0,0) == "@") { print $0; } else if($5< 10) { print $0 }}' | samtools view -Sbo filtered.bam -
    should do the job for you.

    Leave a comment:


  • pguilha
    replied
    bwa aln.
    Yep I'll extract those reads and take a look. In fact, is there a way using samtools to exclude reads that have eg: mapq>10 (in the same way that you can use -F to exclude reads with certain flags?), or will I have to script that?
    Thanks again for taking the time to help me with this.

    Leave a comment:


  • dpryan
    replied
    bwa mem or bwa aln? The mem version can produce chimeric alignments, which might further complicate things. Aside from that, I don't know why that's happening, you might just have to peek at the data.

    Leave a comment:


  • pguilha
    replied
    I'm using bwa, so really don't understand what happened...

    Leave a comment:


  • dpryan
    replied
    What aligner are you using? My only guess is that it gives mates different MAPQ scores, so the filtering can then create orphans. I'm not aware of an aligner that does that, of course...

    Leave a comment:


  • Removing remaining read in pair after samtools mapq filtering

    Hi,

    I am filtering a bam file using samtools (and picard for duplicates) and am struggling to understand the results from flagstat (apologies in advance for the long post)

    This is my sorted bam file before any filtering:

    Code:
    samtools flagstat Sample.sorted.bam
    86744514 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    16772228 + 0 mapped (19.34%:nan%)
    86744514 + 0 paired in sequencing
    43372257 + 0 read1
    43372257 + 0 read2
    14664488 + 0 properly paired (16.91%:nan%)
    15934998 + 0 with itself and mate mapped
    837230 + 0 singletons (0.97%:nan%)
    833002 + 0 with mate mapped to a different chr
    317794 + 0 with mate mapped to a different chr (mapQ>=5)
    Now, if I filter with -f2 to keep only those reads mapping in proper pair, flagstsat gives this:

    Code:
    samtools flagstat Sample.sorted.mapped.bam
    14664488 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    14664488 + 0 mapped (100.00%:nan%)
    14664488 + 0 paired in sequencing
    7332244 + 0 read1
    7332244 + 0 read2
    14664488 + 0 properly paired (100.00%:nan%)
    14664488 + 0 with itself and mate mapped
    0 + 0 singletons (0.00%:nan%)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    Notice how I get the same number of read 1 and read 2, which makes sense to me.

    But if I then run -q10 to remove any reads that have mapq<10, I get this:

    Code:
    samtools flagstat Sample.sorted.mapped_qualfiltered.bam
    10803413 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    10803413 + 0 mapped (100.00%:nan%)
    10803413 + 0 paired in sequencing
    5398606 + 0 read1
    5404807 + 0 read2
    10803413 + 0 properly paired (100.00%:nan%)
    10803413 + 0 with itself and mate mapped
    0 + 0 singletons (0.00%:nan%)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    (btw, I get the exact same result if doing the flag and mapq filters in the same command)

    If I'm not doing something incredibly stupid here, my first question is this: I'm getting a different number of read1 and read2 so I'm assuming samtools is removing only the reads that fail the quality filter and not their mate? Is there a way to tell it to remove those? i can't use flags at this stage as they are still labelled as properly paired even though the mate has been removed.

Latest Articles

Collapse

  • seqadmin
    Best Practices for Single-Cell Sequencing Analysis
    by seqadmin



    While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
    06-06-2024, 07:15 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, 06-21-2024, 07:49 AM
0 responses
14 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-20-2024, 07:23 AM
0 responses
14 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-17-2024, 06:54 AM
0 responses
16 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-14-2024, 07:24 AM
0 responses
25 views
0 likes
Last Post seqadmin  
Working...
X