15
The obvious answer is that different people wrote them. It's fairly common in bioinformatics for people with a computer science background to get frustrated with existing tools and create their own alternative tool (rather than improving an existing tool). Over time, tools with similar initial aims will have popular functionality implemented in them (and ...
14
samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.
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 ...
14
Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details.
Both Kallisto and Salmon ...
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
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 ...
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
Whole genome aliment can be done using Progressive Mauve, LAST or Mummer. For bacteria I used Mauve since it has also very nice visualisation engine. A very new tool is Minimap2, a super fast mapper that supposed beside read mapping be able to handle reference vs reference. However, I do not know how performance of it compares to the tools mentioned above.
...
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
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 | \
...
8
Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency. So samtools merge (use the most up-to-date version, as there have been improvements in merge header handling in the 1.3.x and 1.4.x versions), picard MergeSamFiles, etc.
These tools need to hold all the input BAM files ...
8
The quick way to get the number of alignments on each reference is
samtools idxstats my_bam.bam
Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us ...
8
I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", ...
7
There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal_region
It is therefore normal to have some female-derived reads mapping in Y chromosome.
You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.
7
The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy. For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.
It may be worthwhile to consider aligning ...
7
Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted in silico before alignment to minimize mapping bias.
I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the --score-min option to allow ...
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 ...
7
Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.
But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.
6
Qualimap will do this for you.
Go to qualimap.bioinfo.cipf.es
Run qualimap (default params are fine) on each BAM file
Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.
One thing to watch out for (you don't mention it in your ...
6
I've never tried this myself, so I don't know how easy this is...
One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You ...
6
You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad ...
6
For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.
In studies where they compare ...
6
If I understand your question correctly, then I think for case of pairwise alignment, there is a simple explanation.
I believe the key insight is that: a mismatch should always score better than a gap.*
This follows biologically since the insertion/deletion (indel) rate is roughly 1/10th that of the substitution rate (i.e. the occurrence of single ...
6
Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match.
zcat My_Hiseq_Data.fq.gz | \
grep -Eo '.{0,3}GATCGATC.*' | \
sed -En 's/.*/ \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
grep --color=always GATCGATC
Here are some ...
6
BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.
EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based ...
6
Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.
5
Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.
Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a ...
5
As per my answer to @_julien_roux on twitter:
Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA:
http://pasapipeline.github.io/#A_ComprehensiveTranscriptome
We did something similar in much simpler organisms in our recent paper:...
Only top voted, non community-wiki answers of a minimum length are eligible