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 ...


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 ...


12

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 alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.


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 ...


9

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 particular size". More details can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect The mapping ...


9

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 ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence. Note, this one-liner does not double-count reads because of uniq. samtools view file_sorted.bam | \ ...


7

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 built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it ...


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 ...


7

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 ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get: ACTGACTGACTGTac (lower case denotes soft-clipping) ACTGACTGACTGTTAAGAACGACTACGAC or ACTGACTGACTGTAC (N.B., you've created some ...


6

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 Trinity-generated sequences Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome Merge the BUSCO full ...


6

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' end of exon2 and bases 76-150 aligning to the 5' end of exon3. This would be a split read because it have two alignments (exon2 and exon3) and those alignment ...


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 ...


6

If you're lucky the information is stored in the BAM header in the lines starting with @PG which stores the command lines run on the BAM, use samtools view -H your.bam to display the header.


5

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. Although I am not aware of numbers accounting for these particular reference and type of Illumina reads, the 1000 genomes project has provided the community ...


5

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 is not necessarily MAPQ=0 to another. That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a ...


5

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 assemblies for individual monophyletic clusters. Long answer: Whereas a single reference sequence can make sense for human or mice (to a certain extent), ...


5

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 software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of. Here is a video demonstrating the use of Sashimi plots ...


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

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 and got hold of the gtf link ftp://ftp.ensembl.org/pub/release-80/gtf/homo_sapiens/Homo_sapiens.GRCh38.80.gtf.gz, and the genome file link ftp://ftp.ensembl....


5

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 genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (...


5

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 with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as ...


5

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 that are expensive computationally. On the other hand, finding exact string matches can be performed using much less expensive operations—that is, it can be done ...


5

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 $\tilde{e}$ = long gap extension penalty per every additional gap To understand the first of the two expressions being minimized, see the Wikipedia section on ...


5

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 suspicious by it and your coverage graph. Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage ...


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

Given the high level of multimapping in this region, you'll need to use the -M --primary options if you want to keep many of the alignments. I would be very hesitant to use these numbers as input for DESeq2 or similar programs, since it's fairly questionable whether one should fully trust the "randomness" of the aligner's assignments. I'm more comfortable ...


4

First the Nextera adapters and the custom barcode adapters overlap each other Nextera 1 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG Barcode 1 GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT Nextera 2 GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG Barcode 2 ...


4

Don't restrict the portion of the genome you're mapping to, you'll just increase false-positive alignments by doing so. If you're worried about mapping time, then (A) use more cores or (B) use a different mapper (BSMAP was the quickest out there last I looked, though in my opinion it's fine to wait a bit longer for better quality results (I normally suggest ...


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