Hot answers tagged

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

Arbitrary record access in constant time To get a random record in constant time, it is sufficient to get an arbitrary record in constant time. I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix. Thanks to ...


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

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality ...


9

It is still slow but grep has a -f option to take in a file samtools view inbam.bam | grep -f read_names.txt > read_locs.txt


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

As wkretzsch suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the FASTQ. Index it As much as I typically hesitate to jump to a solution that requires a script or framework (as opposed to just unix command line tools), there is sadly no samtools fqidx (perhaps there should be), and existing answers suggest a ...


7

If you have your counts in a data.frame called counts, something like this might work: filtered.counts <- counts[rowSums(counts==0)<10, ] For example lets assume the following data frame. > A <- c(0,0,0,0,0) > B <- c(0,1,0,0,0) > C <- c(0,2,0,2,0) > D <- c(0,5,1,1,2) > > counts <- data.frame(A=A, B=B, C=C, D=D) > ...


6

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the --split RG argument. Here's some example output from a BAM file with combined read groups: $ samtools stats --...


6

I am commenting on this part: The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. According to its user guide, cutadap is designed this way: it trims off bases from the 3'-end until it sees a base with quality higher than a threshold. This is not a ...


6

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution. EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k*2 or longer is counted. For example in sequence ...


5

You might have to manually strip those auxiliary tags off: samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam - The sed bit is searching for the aux tag and removing everything up to the next tab.


5

You could shuffle the FASTQ once and then read sequences off the top of the file as you need them: gzip -dc input.fastq.gz | paste - - - - | shuf | tr '\t' '\n'| gzip -c > output.fastq.gz I would recommend pigz as a replacement for gzip in the compression step if you have it available. The downside of this approach is that you only get n reads before ...


5

Not sure if my explanation is any good but let's try... First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2). Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550: n = 10000 (n * 0....


4

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses. If you did your variant calling on samples separately and then merged ...


4

using samjdk and invoking the function getReadGroup() getReadGroup() returns The SAMReadGroupRecord from the SAMFileHeader for this SAMRecord, or null if 1) this record has no RG tag, or 2) the header doesn't contain the read group with the given ID.or 3) this record has no SAMFileHeader java -jar dist/samjdk.jar -e 'return record.getReadGroup()==null;' ...


4

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it m: sel = order(apply(m, 1, var), decreasing=TRUE)[1:100] sel then contains the indices into your matrix (or the original dataframe). BTW, I hope you're not filtering by variance ...


4

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be. They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases. This would include N's, but a small deletion or insertion would also cause the last 6 bases to ...


4

Typically I use samtools for operations like this. Specifically I use samtools view with either -r or -R flag depending on the use case. -r STR Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release. -R FILE Output alignments in read groups ...


3

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message.


3

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same ...


3

While this answers explains how to do it I want to address when and why and which thresholds to do it. Filtering the genes with low counts is usually done because the counts are not reliable it would be noise, specially when there are low number of samples these genes disturb power of the analysis. However with 1200 samples having 10 or more samples without ...


3

I've just been generating data like this, so can tell you about why/how missing calls are created in my dataset. There are two main reasons: 1. Sequencing failure When reads don't map across the variant region, then it's impossible to accurately determine a genotype for that region. This will commonly happen just outside the borders of the selected regions ...


3

One possibility is to: reformat the data such that each record is a single line containing the read description, bases, and quality scores pad out each record to a maximum length in each field such that every record in the file is the same number of bytes the total number of records can now be calculated as file size / record size choose a random record ...


3

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.


3

I actually disagree with that.. I guess I write this down as a discussion. The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read. However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads: ...


2

I wrote a tool called sample that you can use to do random sampling without reading the entire file into memory. It can be used where GNU shuf fails for lack of sufficient memory. It requires two passes through the file to do a random sample, but the second pass is generally fast(er) as it uses mmap routines to do cached reads. If you do repeated ...


2

One of the most thorough treatments of this question (or a similar question: grabbing a random subset of reads) was given by Jared Simpson in a blog post a few years ago. http://simpsonlab.github.io/2015/05/19/io-performance/ If you just want to grab a single random read, Jared's benchmarks suggest that seeking to a random position in the file and then ...


2

For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic. Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. ...


2

A complete solution to validate a SAM/BAM file would be picard's ValidateSamFile. This will find many other issues with BAM files. http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile


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