54
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 ...
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 ...
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 "^>" | ...
17
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 ...
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
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:...
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
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
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
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
FASTA (officially) just stores the name of a sequence and the sequence, unofficially people also add comment fields after the name of the sequence. FASTQ was invented to store both sequence and associated quality values (e.g. from sequencing instruments). SAM was invented to store alignments of (small) sequences (e.g. generated from sequencing) with ...
7
I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.
It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was ...
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
This can also be done with regular awk.
awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>
The NR%4==2 count the second line out of every block of 4. length is a built-in that defaults to the length of the line, same as length($0). In this case you can inject you custom printing to the END{} block but countread and ...
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, ...
7
There are a few quick'n'dirty ways depending on the type of data.
In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC
which calcualtes the fraction of reads aligning to exon, intron, intergenic etc.
RNA-seq, if you use a standard aligner such as ...
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
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
255 is the default maximum size of a Counttable in khmer. You want to do the following:
import khmer
counts = khmer.Counttable(31, 1e7, 1) # Feel free to change the parameters
counts.set_use_bigcount(True)
The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (...
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 ...
Only top voted, non community-wiki answers of a minimum length are eligible