user172818
  • Member for 4 years, 8 months
  • Last seen this week
What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
Accepted answer
24 votes

For FASTQ: seqtk fqchk in.fq | head -2 It gives you percentage of "N" bases, not the exact count, though. For FASTA: seqtk comp in.fa | awk '{x+=$9}END{print x}' This command line also works with ...

View answer
What's the most efficient file format for the storage of DNA sequences?
20 votes

The standard and the most common sequence format is FASTA for sure. You can compress it with a compressor. For the ~3GB human genome, gzip reduces the size to ~900MB, depending on the option in use. ...

View answer
Meaning of BWA-MEM MAPQ scores
Accepted answer
14 votes

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

View answer
The state, limitations and comparisons of large variant stores
Accepted answer
13 votes

An epic question. Unfortunately, the short answer is: no, there are no widely used solutions. For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the ...

View answer
What is indel calling and what is its purpose?
10 votes

There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently. <50bp short indels are called from ...

View answer
Single-sample vs. joint genotyping
Accepted answer
10 votes

Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you ...

View answer
What are the de facto required fields in a SAM/BAM read group?
Accepted answer
9 votes

The sample tag (i.e. SM) was a mandatory tag in the initial SAM spec (see the .pages file; you need a mac to open it). When transitioned to Latex, this requirement was mysteriously dropped. Picard is ...

View answer
Library for computing BWT-based alignments
9 votes

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

View answer
Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value?
9 votes

Let f and r be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base c in a proper 2-bit encoding, we update the two integers as follows: f = (f<&...

View answer
Are variant calling files personally identifiable information?
8 votes

Yes, definitely identifiable. The combination of ~80 unlinked common SNPs can be fairly unique in the entire human population, let alone the whole VCF file. EDIT: 30 in the original answer is an ...

View answer
How to get fasta alignment file from SAM/BAM file?
8 votes

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

View answer
Why Ti/Tv ratio?
Accepted answer
8 votes

I more like to use "ts/tv" for transition-to-transversion ratio. This abbreviation had been used in phylogenetics. When NGS came along, some important developers started to use "ti/tv", but I am still ...

View answer
How do kmer counters determine which kmer is 'canonical'?
Accepted answer
8 votes

When a k-mer is identical to its reverse complement, both are canonical. Note that a canonical k-mer is a sequence, irrelevant to its position(s) in the input string. More precisely, give a string $s$,...

View answer
What Ensembl genome version should I use for alignments? (e.g. toplevel.fa vs. primary_assembly.fa)
Accepted answer
8 votes

Generally, you should use the soft-masked or unmasked primary assembly. Cross-species whole-genome aligners, especially older ones, do need to know soft-masked regions; otherwise they can be ...

View answer
Accuracy of the original human DNA datasets sequenced by Human Genome Project?
8 votes

While the quality of the reference human assembly keeps improving, there are still misassemblies in it. A common problem is recent segmental duplications are occasionally collapsed into one sequence ...

View answer
How should the popular press compare similarity of genomes?
Accepted answer
7 votes

1–4% is from an evolution point of view. 99.7% is from a sequence point of view. Because they are measuring different things, they can be compatible with each other. The correct interpretation is: 1–4%...

View answer
How to write a hash function for canonical kmers?
Accepted answer
7 votes

You can convert any string hash function to a "canonical" DNA string hash function. Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$...

View answer
How do PCR duplicates arise and why is it important to remove them for NGS analysis?
7 votes

PCR polymerases introduce errors. When an error arises in the first few cycles of amplifications, it will appear in a reasonably high fraction of DNA fragments in the library. After sequencing, you ...

View answer
Why sequence the human genome at 30x coverage?
7 votes

Solexa Inc. sequenced NA12878 chrX to ~30x in early 2007, which later became part of Bentley (2008). This, I believe, was the first time that 30x showed up. I don't recall they had a particular reason ...

View answer
Rule Extraction from nnet results
7 votes

R's nnet package only supports fully connected neural networks with one hidden layer. This is the most primitive type of network. I doubt it will work well for promotors finding. In addition, such ...

View answer
Fast processing of fastq data
6 votes

Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose. C++, Rust, Go, Nim and Julia can ...

View answer
Can blat use more than one core/CPU to speed up the alignment?
Accepted answer
6 votes

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

View answer
What is the most compact data structure for canonical k-mers with the fastest lookup time?
Accepted answer
6 votes

The question and the accepted answer are not about k-mer data structure at all, which I will explain in detail below. I will first answer the actual question OP intends to ask. The simplest way to ...

View answer
Are these standard species abbreviations and how to look up others?
6 votes

These are not. The closest to "standard" is the 5-character abbreviation by Swissprot. It names 25,886 species and has been used for decades. It is also easy to remember for some common species such ...

View answer
Remove/delete sequences by ID from multifasta
6 votes

Suppose you keep sequence names in ids.txt and sequences in seq.fa: awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

View answer
Counting repeated kmers sequences that match at least x % of reads sequence
6 votes

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

View answer
Why does cutadapt remove low quality bases from the ends of reads only?
Accepted answer
6 votes

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

View answer
How to calculate the memory usage of storing kmers in RAM
6 votes

I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM? Hash table based k-mer counters only keep k-mers seen in data. For $16<k\le32$, you need 64-bit ...

View answer
Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?
6 votes

I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, ...

View answer
Do variant calls change when you call from CRAM?
Accepted answer
6 votes

By default, a CRAM you create with samtools is lossless. It typically halves the input BAM in terms of file size. If you want to compress more, you can let samtools convert most read names to integers....

View answer