36
Let’s start with what they have in common: All three formats store
sequence data, and
sequence metadata.
Furthermore, all three formats are text-based.
However, beyond that all three formats are different and serve different purposes.
Let’s start with the simplest format:
FASTA
FASTA stores a variable number of sequence records, and for each record it ...
21
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 ...
17
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 "^>" | ...
16
The maternal and paternal copies of a chromosome are called haplotypes. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states.
Response to Q1
Your question, in other words, is: Why do .bam files not differentiate between haplotypes?
Your question ...
14
In a nutshell,
FASTA file format is a DNA sequence format for specifying or representing DNA sequences and was first described by Pearson (Pearson,W.R. and Lipman,D.J. (1988) Improved tools for biological sequence comparison. Proc. Natl Acad. Sci. USA, 85, 2444–2448)
FASTQ is another DNA sequence file format that extends the FASTA format with the ability ...
13
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 ...
12
I attended a talk by Josh Quick at PoreCampAU 2017, in which he discussed some common barriers to getting both good sequencing yield and long read length. It mostly boils down to being more careful with the sample preparation. Bear in mind that the MinION will still sequence a dirty sample, it will just be at a reduced yield. Here are my notes from that talk:...
9
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
awk '{if(NR%4==1) print $1, $5}' file.fastq | sed -e "s/ start_time=/, /" -e "s/^@//"
The awk command gets the first of every 4 lines, printing the first and fifth "word". sed is then used to strip the initial @ and replace start_time= with ,. The output on your example file is:
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-...
9
Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ. For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.
If the FASTQ is malformed, then you should really ask yourself ...
8
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.
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
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
Incidentally, the first part of your question is something you could have looked up yourself as the first hits on Google of "NAME format" point you to primers on Wikipedia, no less. In future, please do that before asking a question.
FASTA
FASTQ
SAM
FASTA (officially) just stores the name of a sequence and the sequence, inofficially people also add comment ...
7
Script to run everything in a loop for extracting tar.gz files into fastq and to bam with alignment?
This should work for you:
for i in *.tar.gz
do
dir2="/destination/directory"
tar xvzf $i -C $TMPDIR
files=($TMPDIR/*)
module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done
(If this is the ...
7
You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads. The library contains an example program to demonstrate how to use it with multiple threads. It's worth noting that kseq is already very fast, ...
6
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 ...
6
Since the string start_time will only appear on the header line, or else you don't have a valid fastq file, you can simply do:
$ perl -ne '/^@(\S+).*start_time=(.*)/ && print "$1, $2\n"' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720,2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c,2017-07-04T06:56:41Z
Alternatively, since you mentioned ...
6
As @AaronBerlin mentioned, you didn't remove reads that were completely trimmed. Next time use the --minimum-length option and set it to something reasonable, like 20. Alternatively, use "Trim Galore!", which is a wrapper around cutadapt that has more reasonable defaults.
6
NOTICE:
I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.
This is something I have been meaning to get around to for a while, so thanks for the prompt.
I have ...
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 ...
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
If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command:
sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'
The -i.bak makes sed modify the original file and create a backup copy with the same name and the extension .bak. The command just means "delete the line matching the pattern and ...
6
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 be as fast.
// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
// Compile with: gcc -O2 -o myprog this-file.c -lz
#include <zlib.h> ...
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
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 ...
5
I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:
A fastq file (optionally compressed)
A bam file
The ...
5
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 ...
5
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 -c fastx makes ...
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: ...
Only top voted, non community-wiki answers of a minimum length are eligible
