16
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,485
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 ...
- 1,899
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.6k
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 ...
- 1,899
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,141
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 ...
- 13.8k
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,485
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 ...
- 663
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 ...
- 13.8k
7
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 ...
- 13.8k
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.6k
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 ...
- 13.8k
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,060
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....
- 827
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,485
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
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 ...
- 1,899
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 ...
- 76
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.6k
5
votes
Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?
With bwa-mem or minimap2, the recommended way is
samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p -
This passes the BC...
- 6,485
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
$\...
- 155
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 ...
- 3,059
Only top scored, non community-wiki answers of a minimum length are eligible