# Tag Info

24

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 FASTQ, but it will be slower as awk is slow. EDIT: ok, based on @BaCH's reminder, here we go (you need kseq.h to compile): // to compile: gcc -O2 -o count-N ... 23 I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions: import string tab = string.maketrans("ACTG", "TGAC") def reverse_complement_table(seq): return seq.translate(tab)[::-1] The thing with hashing is that it adds a good bit of overhead for a replacement set this small. For what it's worth, I ... 18 5 hours and no benchmarks posted? I am sorely disappointed. I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are: R with the ShortRead package (even if not the fastest, certainly a super convenient method). A pipeline of grep -v "^>" | tr -cd A | wc -c A pipeline of grep -v "^>" | ... 14 It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer. People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C. I ... 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 I think that this should be pretty fast: FASTA: grep -v "^>" seqs.fa | tr -cd N | wc -c FASTQ: sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c See this answer on SO about how to count characters in BASH using different approaches. 9 Here's a Cython approach that might suggest a generic approach to speeding up Python work. If you're manipulating (ASCII) character strings and performance is a design consideration, then C or Perl are probably preferred options to Python. In any case, this Cython test uses Python 3.6.3:$ python --version Python 3.6.3 :: Anaconda custom (64-bit) To ...

8

FASTQ As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is: sed -n '2~4p' seqs.fastq | grep -io N | wc -l sed -n '2~4p' will print every fourth line grep -o N will output a line with N for every matching symbol wc -l will count the lines I suspect this python approach ...

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

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve. Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl ...

7

The most reliable and simplest way is probably using Biopython: from Bio.Seq import Seq def my_reverse_complement(seq): return Seq(seq).reverse_complement() print(my_reverse_complement('ATGCGTA')) # TACGCAT As Devon has already said here using Biopython isn't as fast as the naive Python solution, and I also tested that shown here with ipython. ...

6

Decompression of gzipped FASTQ is the main issue If we take real world gzipped FASTQ (which as the OP suggested would be beneficial) rather than trivial FASTA as the starting point then the real issue is actually decompressing the file not counting the Ns and in this case the C program count-N is no longer the fastest solution. Additionally it would be ...

6

Using bioawk With bioawk (counting "A" in the C. elegans genome, because there seem to be no "N" in this file), on my computer: $time bioawk -c fastx '{n+=gsub(/A/, "",$seq)} END {print n}' genome.fa 32371810 real 0m1.645s user 0m1.548s sys 0m0.088s bioawk is an extension of awk with convenient parsing options. For instance ...

6

The following is more than twice as fast; however, wc counts newline characters as well. We thus need to subtract the line count from the base count: fix_base_count() { local counts=($(cat)) echo "${counts[0]} $((${counts[1]} - ${counts[0]}))" } gunzip -c "$file" \ | awk 'NR % 4 == 2' \ | wc -cl \ | fix_base_count However, the ...

6

pigz | awk | wc is the fastest method First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer. I've chosen this file: ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz as my test file, the correct answers being: Number of reads: 67051220 Number of bases in ...

6

If you want to learn how to do things with command-line tools, you can linearize the FASTA with awk, pipe to grep to filter for items of interest named in patterns.txt, then pipe to tr to delinearize: $awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s",$0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff ... 6 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 6 Just today I wrote a script to do exactly this using Biopython. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script filter_fasta_by_list_of_headers.py is : #!/usr/bin/env python3 from Bio import SeqIO import sys ffile = SeqIO.parse(sys.argv[1], "fasta") header_set = set(line.strip() for ... 5 I get fairly quick results with my fastx-length.pl script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics: time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > /dev/null Total sequences: 301135 Total length: 283.902419 Mb Longest sequence: 5.601 kb Shortest sequence: ... 5 If it is raw speed you're after, then writing an own little C/C++ program is probably what you need to do. Fortunately, the worst part (a fast and reliable parser) has already been tackled: the readfq from Heng Li is probably the fastest FASTA/FASTQ parser around. And it's easy to use, the example on GitHub can easily be expanded to do what you need. Just ... 5 Honestly, the easiest way (especially for FASTQ) is probably to use a dedicated parsing library, such as R/Bioconductor: suppressMessages(library(ShortRead)) seq = readFasta(commandArgs(TRUE)[1]) # or readFastq cat(colSums(alphabetFrequency(sread(seq))[, 'N', drop = FALSE]), '\n') This may not be the fastest, but it’s pretty fast, since the relevant ... 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 Perl should be fairly fast with this when using a hash set to store the list of lines. A structure like this also works for subsetting based on a field value, where the comparison would be with the field rather than "$.": #!/usr/bin/perl use strict; use warnings; my $lines_file =$ARGV[0]; my %include_lines = (); open my $lines_fh, '<',$lines_file or ...

5

Another python extension but without cython. The source code is available at the bottom of this answer or from this gist. On Mac with Python3: string 0.38166s total, 1310067.4 strings per second 84.4% increase over baseline seqpy 0.20524s total, 2436182.7 strings per second 91.6% increase over baseline On Linux with Python2 (seqpy is the ...

4

Some related questions appear in other sites, with potentially interesting solutions, which I report here: To sample approximately 1% of the non-empty lines: awk 'BEGIN {srand()} !/^$/ { if (rand() <= .01) print$0}' input_file (from https://stackoverflow.com/a/692321/1878788) To select 1000 random lines: shuf -n 1000 input_file (from https://...

3

For FASTA files, I've implemented a relatively efficient method in pyfaidx v0.4.9.1. This post made me realize that my previous code was quite slow and easy to replace: $pip install pyfaidx$ time faidx -i nucleotide ~/Downloads/hg38.fa name start end A T C G N others chr1 1 248956422 67070277 67244164 48055043 48111528 ...

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

Use a bytearray instead of a string and then employ maketrans to translate You do not need the more advanced string encoding capabilities of string to store a string of bases, but you're still paying for it in performance. Devon Ryan's suggestion of maketrans is the huge improvement, 10x faster than your naive implementation. Using the same approach, but ...

3

Here is a revision of my original Cython answer which incorporates my suggestion to use a char lookup array: from libc.stdlib cimport malloc cdef int seq_len = 17 cdef char *seq_dest = <char *>malloc(seq_len + 1) seq_dest[seq_len] = '\0' cdef char *basemap = [ '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',...

3

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do: FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa

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