Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • BWA - read identity

    Hi,

    I'm using BWA to map my short sequencing reads (Illumina MiSeq) on a reference sequence. Is this possible to set the minimum % of identity with BWA? For exemple, if I set 90%, then all reads under this threshold will be discarded. I know the parameter "-n" can do something like this but the manual is not really clear (for me). Is "-n" the minimal absolute mismatch in term of base pairs?

    Thank you!

    Antony

  • #2
    Since nobody else has answered this, I will mention that I do not know of a way to generically do that in BWA. The manual says:
    -n NUM Maximum edit distance if the value is INT, or the fraction of missing alignments given 2% uniform base error rate if FLOAT. In the latter case, the maximum edit distance is automatically chosen for different read lengths. [0.04]
    ...which is clear in the case of int but unclear in the case of float. If all your reads are the same length you can calculate the edit distance that would correspond to that identity, though: e.g. 150bp * (1-0.90) = 15.

    Alternately, you can use BBMap, which has an "idfilter" flag - e.g. "idfilter=0.90", which will discard reads that don't map with at least 90% identity, regardless of their length. It also has an "idtag" flag, which will annotate all the alignments in the sam file with the percent identity, so you can manually filter them if you want.

    Comment


    • #3
      [This is a hack. I'm putting it here more for curiosity than anything else]

      On the aligned bam file you could capture the read length and the NM tag ("Edit distance to the reference, including ambiguous bases but excluding clipping") and print the alignment if L/NM < threshold. With samtools/awk:

      Code:
      samtools view -h aln.bam \
      | awk -v maxpct=0.05 -v OFS='\t' -v FS='\t' '{
          if($0 ~ /^@/){
              print $0
          } else {
              xnm=gensub(/.*\tNM:i:/, "", "g", $0); 
              nm=gensub(/\t.*/, "", "g", xnm); 
              pct=nm/length($10);
              if(pct < maxpct) {
                  print $0, "XP:f:"pct
              }
          }
      }' | less
      (The L/NM goes in tag XP)

      Comment

      Latest Articles

      Collapse

      • seqadmin
        Choosing Between NGS and qPCR
        by seqadmin



        Next-generation sequencing (NGS) and quantitative polymerase chain reaction (qPCR) are essential techniques for investigating the genome, transcriptome, and epigenome. In many cases, choosing the appropriate technique is straightforward, but in others, it can be more challenging to determine the most effective option. A simple distinction is that smaller, more focused projects are typically better suited for qPCR, while larger, more complex datasets benefit from NGS. However,...
        10-18-2024, 07:11 AM
      • seqadmin
        Non-Coding RNA Research and Technologies
        by seqadmin




        Non-coding RNAs (ncRNAs) do not code for proteins but play important roles in numerous cellular processes including gene silencing, developmental pathways, and more. There are numerous types including microRNA (miRNA), long ncRNA (lncRNA), circular RNA (circRNA), and more. In this article, we discuss innovative ncRNA research and explore recent technological advancements that improve the study of ncRNAs.

        Nobel Prize for MicroRNA Discovery
        This week,...
        10-07-2024, 08:07 AM

      ad_right_rmr

      Collapse

      News

      Collapse

      Topics Statistics Last Post
      Started by seqadmin, 11-01-2024, 06:09 AM
      0 responses
      18 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 10-30-2024, 05:31 AM
      0 responses
      18 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 10-24-2024, 06:58 AM
      0 responses
      24 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 10-23-2024, 08:43 AM
      0 responses
      53 views
      0 likes
      Last Post seqadmin  
      Working...
      X