15 votes
Accepted

Meaning of BWA-MEM MAPQ scores

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway. Strictly speaking, you have ...
  • 6,121
12 votes

Difference between BWA-backtrack and BWA-MEM

TL;DR: BWA-backtrack is based on backtracking. This approach is appropriate only when the dissimilarity between the reads and the reference is low, or when you want to find all best hits or enumerate ...
12 votes
Accepted

Definition of "seed" in sequence alignment

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the ...
  • 19.4k
11 votes
Accepted

Library for computing BWT-based alignments

First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as ...
10 votes
Accepted

How to count the number of mapped read in 100-bp window from a BAM/SAM file

one-liner Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ...
  • 3,111
9 votes
Accepted

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

1. Adapter Trimming One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of ...
  • 12.5k
9 votes

Library for computing BWT-based alignments

BWA-MEM can be used as a library. File bwa/example.c shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header bwa/bwamem.h contains ...
  • 6,121
9 votes

Difference between BWA-backtrack and BWA-MEM

To quote the Introduction to BWA on sourceforge: BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the ...
9 votes

What is the difference between SAM mapping quality and Blast E-value?

The E-value and the mapping qualities are two very different things. The E-value is "a parameter that describes the number of hits one can 'expect' to see by chance when searching a database of a ...
  • 952
7 votes
Accepted

Why is bwa-mem the standard algorithm when using bwa?

bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads. From the bwa man page (presumably in ...
  • 12.5k
7 votes
Accepted

What is local realignment and what is the problem it solves?

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ...
  • 19.4k
6 votes

Filter Trinity transcriptome based on RNASeq reads

I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows: Run BUSCO in short/transcript mode on the ...
  • 12.5k
6 votes

Meaning of BWA-MEM MAPQ scores

Those numbers are not arbitrarily picked (well... maybe 255/60/40 is arbitrarily picked). To convert from log10 Q values like these (also used for error rates in FASTQ files) to probabilities, divide ...
  • 12.5k
6 votes
Accepted

What is mate rescue in bwa mem?

Using -SP is equivalent to running bwa mem on each of the two mates as if they are single-end reads, but it formats the output ...
  • 76
6 votes

Definition of "seed" in sequence alignment

Devon's answer gives a good, concise definition. But it's also helpful to consider why seed-and-extend is used and what benefits it provides. Finding approximate string matches requires operations ...
5 votes

what percentage of the human genome is MAPQ=0?

As has been said before, mappability to the 'human genome' depends on a number of factors, among these the reference version and type of reads, for which you are interested in GRCh38 and 2x150bp reads....
  • 705
5 votes

what percentage of the human genome is MAPQ=0?

Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper ...
  • 6,121
5 votes

Visualisation of long read RNA-Seq splicing

DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots. Edit: not sure why this answer is being downvoted, unless it's because the ...
5 votes

Which reference to use for read mapping for popular model organisms

I will answer one of the points – E.coli. TL;DR Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple ...
5 votes

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

The GATK is probably the most widely used program for this. There are several workflows you can use depending on what types of variants you're interested in (germline or somatic, short variants, copy ...
5 votes
Accepted

RapMap: reference transcriptome for simulated reads

Thanks for your interest in RapMap. At that time we were using flux simulator for simulating read sequence data. We used the genome and gtf file together as an input to flux. I dug into the scripts ...
5 votes

Why total RNA-seq usually yields low mapping rate?

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple ...
  • 61
5 votes
Accepted

What are "split reads" and "intron clusters?"

split reads - These are read that have two or more alignments to the reference from unique region of the read. In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' ...
  • 2,574
5 votes

Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools. But regarding your actual question, STAR can accept SAM/BAM as input ...
  • 19.4k
5 votes

What do the symbols mean in minimap2's gap cost equation?

Here is the translation: $l$ = gap length $q$ = initial gap penalty for shorter gaps $e$ = short gap extension penalty per every additional gap $\tilde{q}$ = initial gap penalty for longer gaps $\...
5 votes

Contamination on genome assembly

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things: Run something like FRC_align. Check if there's a clear signal between regions flagged as ...
5 votes
Accepted

Understanding PAF format

As @user172818 mentioned in a comment to your question, it looks like you might have used the query as the reference; they should be switched round. A minimap2 mapping command is typically represented ...
5 votes
Accepted

Kallisto error: index input file could not be opened!

Usually errors like these occur when an input file doesn't exist or is corrupt. But you might also see these errors if you don't have permissions to access the file for whatever reason. Although ...
  • 2,156
4 votes

Library for computing BWT-based alignments

SeqAn supports BWT tables for use with their parametrisable alignment algorithms. To use it, follow the general outline for building a SeqAn short-read aligner, and use the ...

Only top scored, non community-wiki answers of a minimum length are eligible