Hot answers tagged

14

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 two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the ...


13

To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields: samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam Explanation follows... I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM ...


12

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 all possible alignments up to a specified number of errors. In all other situations, BWA-MEM is preferable as it can, thanks to its sophisticated strategy ...


11

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 a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned ...


9

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 basic documentation. Paired-end mapping is doable, but not well exposed. Several teams, including GATK and MS genomics, have been using bwa-mem as a library.


9

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 human genome. It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ...


8

The SAM format originally had only M, I, D, N, S, H, and P CIGAR operators. See the original SAM specification (if you can view Apple Pages documents) and Table 1 in The Sequence Alignment/Map format and SAMtools (Li et al, 2009). This was in line with previous tools using CIGAR strings, notably exonerate which introduced them with just the M, I, and D ...


7

Assuming the read name lines only contain things like @JQ714243-920-1 and nothing else, then sed -i '1~4 s/-[12]$//' file.fastq (N.B., example updated a bit for safety, thanks to @terdon) will remove the -1 and -2 suffixes. However, you should try to figure out why the read names were munged into that form to begin with (at least Illumina machines won't ...


7

bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads. From the bwa man page (presumably in Heng Li's own words): BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM. The first ...


6

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


6

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 as a proper paired-end output, i.e. with all pair-related flags added properly. Without -SP, by default bwa mem forces an alignment of a poorly aligned read if its mate is aligned somewhere. -SP turns off the forced alignment. We use ...


5

GATK best practices are explicably meant to consume BWA MEM generated BAM. Whilst SNAP may be faster, the Broad will not have tested it for compatibility with GATK as such you can't guaranty using it won't have unexpected consequences. As such you'd be better off using BWA MEM because I assume accurately called variation is always better than fast and ...


5

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 the number by 10, negate it, then raise 10 to the power of the result. Another way of looking at it is to consider the decade to mean the number of 9s in the [...


5

I don't believe that is an option. Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.


5

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to: bwa mem (...options) | samtools view -o out.bam samtools fastq -f 4 out.bam > unmatched.fastq


5

This is an error from bwa saying that it can't find your input files: Here's a checklist: Check that you have no whitespace after your \ newline breaks. I.e. \ is the very last character Check that you have no whitespace in any path names Quote your variables: $REF_PATH -> "$REF_PATH" Use set -x at the top of your bash script to debug. In ...


4

A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference. If you are concerned about space and have a reference sequence available, you can store mapped reads ...


4

You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c. ...


4

Papers The space-efficient BWT construction algorithm is in bwt_gen.c. One of the first things you see there is the copyright line: Copyright (C) 2004, Wong Chi Kwong. With that name, you can find the relevant papers: Lam et al: Compressed indexing and local alignment of DNA (Bioinformatics, 2008). This paper describes the BWT-SW alignment algorithm. ...


4

No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the ...


4

For much faster and stable$ implementation, sambamba can handle this using following one-liner: sambamba view -t 12 -h -f bam -F "mapping_quality >= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam See manpage and sambamba synatx for filters for details on further usage. $: Besides ...


4

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 FMIndex specialisation instead of — as in the example — the IndexQGram.


4

See here, in particular slide #10 from the tutorial: bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.


4

Usually it doesn't actually matter, but if you want it to be very correct: RGID=SomeLibraryID \ RGLB=SomeLibraryID \ RGPL=illumina \ RGPU=Some_Sequencer_ID \ RGSM=SomeSampleID That assumes you have one library per sample and aren't splitting the output BAM files by lane.


3

Yes, there bwa-mem was published as a preprint BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference position x with the best extension score achieved at query position y. ... Secondly, while extending a seed, BWA-MEM tries to keep track of the best ...


3

In the bwa source code, file bwt.{c,h} implement basic operations on .bwt and .sa files. They depend on utils.{c,h}, but you can create a fake utils.h to cut this dependency as is shown below. Then you only need three files: kvec.h, bwt.h and bwt.c to access .bwt and .sa. The following example demonstrates how to find super-maximal exact matches (SMEMs). ...


3

bwa mem is using exactly the same meaning of the 0x200 flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for ...


3

Minimap2 doesn't work well with 32bp reads by design. Older mappers like bowtie1 and bwa-aln will be better.


3

[W::sam_read1] Parse error at line 1 samtools sort: truncated file. Aborting SAMtools sort has been unable to parse its input, which it thought was SAM (mostly because it couldn't be recognised as another format e.g. BAM). This is because sed 's/^/LP1-/' is putting LP1- at the front of every line. If the output of samtools fixmate is SAM, then this LP1 is ...


3

You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment. From bwa docs: If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair. BWA will not search through FASTQ files ...


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