A .fq or fastq file is an ASCII encoded text file that stores DNA or RNA sequences and their corresponding IDs and quality scores.

It uses unix newlines and consists of 4 lines per sequence unless wrapping occurs due to sequence length.

The first line begins with an "@" followed by an identifier (ID) which acts as a label for the read/sequence.

The second line represents a DNA or RNA sequence, and should consist only of standard bases, and IUPAC ambiguity codes (ACTGNURYSWKMBDHV). This line must be wrapped with newlines if the read is longer than 80nt.

The third line must be a single "+" which signifies the end of the sequence.

The fourth line is a quality score string showing the quality of each base in the prior sequence, represented as the ASCII character corresponding to the quality Phred score + 33. Phred scores must be 0 and 60 (ASCII chars 33 aka "!" to 93 aka "]").

The predominant quality score here is "#" which has an ASCII value of 35.

A .fq or fastq file is an ASCII encoded text file that stores DNA or RNA sequences and their corresponding IDs and quality scores.

It uses unix newlines and consists of 4 lines per sequence unless wrapping occurs due to sequence length.

The first line begins with an "@" followed by an identifier (ID) which acts as a label for the read/sequence.

The second line represents a DNA or RNA sequence, and should consist only of standard bases, and IUPAC ambiguity codes (ACTGNURYSWKMBDHV). This line must be wrapped with newlines if the read is longer than 80nt.

The third line must be a single "+" which signifies the end of the sequence.

The fourth line is a quality score string showing the quality of each base in the prior sequence, represented as the ASCII character corresponding to the quality Phred score + 33. Phred scores must be 0 and 60 (ASCII chars 33 aka "!" to 93 aka "]").

The predominant quality score here is "#" which has an ASCII value of 35.

Hi Friend!

Some extra info for you:

the first string after the char @ is unique for the sequencing run. This is then followed by the x and y coordinates of the particular read on the sequencing plate, so they indentify unique your read. The second part is well explained by Ms. Nagraj.

Here I share you my simple and general short-read alignment consensus calling pipeline, nevertheless i use this only for small genomes like viruses or bacteria. It also includes the mapping itself, but you can modify it, it isn't needed. You will need some experience with shell scripting, but if vcfutils, bwa, seqtk and samtools are in the path, it's no big deal.

#!/bin/bash

ref=$1

reads=$2

bam=${reads//\.fastq/}.on.${ref//\.fasta}

bwa index $ref

bwa mem -t 8 $ref $reads | samtools view -bS -F 4 - | samtools sort - $bam

samtools index $bam.bam

samtools mpileup -uf $ref $bam.bam | bcftools view -cg - | vcfutils.pl vcf2fq | seqtk seq -A -l 70 - | sed -r "s/(>)(.*)/\1$bam.conensus/g" > $bam.consensus.fasta

1) What is the .bam file you have? Normally .bam files contain reads mapped to some reference, but you say you sequenced the genome?

2) You mention .fastq files. Are these the reads you got from the sequencing company, or did you extract them from the .bam file?

3) Normally raw sequencing reads from .fastq need to be assembled to get the draft genome sequence (in the form of contigs and scaffolds/super-contigs). Was assembly performed on the reads you have?

Lowercase and uppercase a/t/g/c are sometimes used as a proxy for base calling OR mapping/assembly quality. If this is the case, then upper-case means higher quality. .fastq files normally only use uppercase letters.

I'm not sure I understand your question about stretches of 200bps... Do you want to close the "gaps" in scaffolds? (Where "gap" means a stretch of Ns, e.g. 200bp long) If so, then you may try GapCloser (part of the SoapDenovo tools) and other similar tools. However, for IonTorrent these tools most probably will not help at all, as you have few long reads instead of many short reads.

You may want to add a good-coverage short-reads library (e.g. Illumina paired-end), then either perform hybrid assembly, or use shorter reads for gap closing.

Dear Daniel,

Thank you very much for your help and support.

I obtained the consensus sequence and I would like to extract sequences >200 bp that are present between runs of n's. For example I have nnnnnnnnnATCGGTCGGTAACCC....>200bpnnnnnnnnnnn and I want only the DNA seq. Is there any way (perl script or other) that do this work?

I need to obtain conting >200 bp to submit to genbank database.

Thank you again for your help.

Since the sequencing was done with Ion PGM, it's impossible to see any lowercase basecalls in the fastq file itself. It's logic, because the basecall uncertainities are already marked with the quality scores. I don't say its impossible to get lowercase bases in the fastq files, but i have seen some PGM runs, and i don't remember such a phenomen. However, if the basecalling is made with some custom flowgram handling tool, i have no idea.

He talks about a bam file, so one way or another but he got a mapping against some reference. I use the bwa aligner for mapping purposes, and now i checked some bam files, and not even in them are lowercase bases.

My coworker often (or mostly...) want to have the consensus of these mappings againts the reference (i dont think its a good idea, but they are biologists, so i do it :) ), and i make that with the pipeline in my previous comment. That gives you a fasta file, and this fasta includes lowercase letters indeed, and in this case, they stand for low coverage and/or low confidence. The 'n'-s are then on the sites, where the coverage is too low (0), but because of the reference, you suppose that there is a nucleotide, you just don't know, which one. That is, what N or n stands for.

I don't know, what do you mean wit that "extrapolate" term, it's a very unproper term on my opinion. I think, your problems are the unmapped regions on the reference. That's not so easy as it seems. In the most common cases, if you assume a sequencing error there, the simpliest solution is to create a primer and let sequence that gap with the good old sanger technology. That should work for a few hundred bases. If there are too many of those gaps, you might have to resequence the thing with an other technology, because there are technology biases of every one and they "don't like" some specific regions of the dna, but at least they are not common. GC bias is for example a famous thing that infulences the sequencing depth.

But the more distant are the species you want to sequence, the greater is the probability, that the missing piece in your consensus is really missing in your sample, so thats a smaller-greater structural variation. For this reason, you have to prove the coverage of the mapping (Tablet is a great tool for working with bams), because coverage differences can also mean copy number variatios, what you don't see from just a fasta file.

Thats why we use a de novo assembler, that assembles the genome without a reference, but thats a very long story to explain here. If you expect more than just snps in your experiment, you should use de novo assembly, for Ion (and 454) data mostly Mira or Newbler, for Illumina some de bruijn graph based one such as abyss or velvet.

The "cyrillic characters" can also be found in the fastq and in the bam files, but the meaning is here already well explained.

You should read at least this:

http://en.wikibooks.org/wiki/Next_Generation_Sequencing_%28NGS%29

book to get the basics, because without brain, you only waste your resources.

Dear Bogdan,

I sequenced a fungal genome using Ion Torrent machine and I assembled the reads against a reference genome (closely related to the sequenced fungus). Then I used the .bam file to extract consensus sequences. Now I am interested to obtain sequences >200 bp to submit to genbank.

I am new in this field and I would like to assemble my data to obtain contigs. Unfortunately I do know to perform a de novo assembly. Could you help me? are you interested to collaborate with me?

Seqtk. If you work often with NGS data, you should keep it in your mind.

for example:

cat input.fasta | seqtk cutN -n 5 - | seqtk seq -L 200 - > output.fasta

where the sequences are cutted if there are at least 5 n-s in between and then the shorter than 200nt ones are dropped. Welcome :)

Dear Daniel,

thank you again for your help.

As I told you in my last reply I am new in this field and I'm looking to obtain assembled contigs from my fungus because curently there are no genomic data on this species. I agree that the best way is de novo assembling but I am not able to perform it. Therefore I decided to map my reads against a reference genome and extract sequences from .bam file according to your suggestions. Here is an exemple of the obtained file:

gtgtcnnnnnnnnnnnnnnnnnnnnatgacaatgcctggccnnnnnnnnnnnnnn

nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnncgctagc

tagctccatagaaaacccacgGGAAGGCGAGTGCAGGGTGATGCGGTTGCCCAGGTTGCGCATAATGTTGTAGCCGTAAGTCCAGAGACCGATGACCAGGGCAGCGGCACCAAAGACCAGAATCCAGATGGGCACATCGGCCTTGCCCTTCTCGGGGATGGCACCGTCGTTCCAGATCTGGAAAATGGTGGCATAGGGGCCAATGGCGTTGGCAATATCGTTGGCTCCGTGCGTGAACGAGGCCGTGGCGGCCGTCATGACCTGCAGGAAGGTGAACATGTACTCGCTCTTGTTGTCGTAGTGGGCAGCACGGGCGTGCATGGCCTCCAGGTCGCCAGTAAGCATGTCGCGCTTGTCCTTTTGCGAGTTCACAACATCGACCTCAACGCCGCGGAAGATGGCGTGCTTGAAGTACcaaaagaggaccgagccggaataccacttgccattgggcttGgggccgacaatggacttctgnnnnnnnnnnnncttctcgttttggatnnnnnnnnnnnnnnnnnnnnggccacgtcgtgccgggttgttcacnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

nnnctcactgacaacggcgcgcgcacgccgtgtatccAGCTGCTCCTGGGTCAGATGGCC

CTGGTAGAAGTTGCGAGTGCCACGGTAGCCCGGAGGAGGCGGGCCAGGTTCGGGACGA

As you can see I have consensus sequences between n's that I want "extract" (I used the term extrapolate before...sorry). I will re-align these contigs against my reference genome with MAUVE to get a better view of my assembled contigs.

Then I will try to close gaps using sanger technology.

However now I am interested to submit my data to genbank and only contigs >200 bp can be deposited in the database.

Given that you are very experienced in this field could you give me more suggestions. If you want we can also collaborate.

Many thanks

cat input.fasta | seqtk cutN -n 5 - | seqtk seq -L 200 - > output.fasta

This should do the work for you, from the linux cmd of course. Write me, and I look what i can do for you, but notice, I never worked with fungi.

bartha dot daniel at agrar dot mta dot hu

Hi Orazio,

given the fact that the sequence (and maybe annotation?) of a closely related fungus is already available, I suggest you to try Mauve (http://gel.ahabs.wisc.edu/mauve/).

Hi friend!

People have already given you nice hints about your file.

A .fq file is basically a FASTA and Quality file merged together.

This file includes details about quality of each base, the ASCII characters that you may have seen.

To process the data, use PRINSEQ.

It has a labyrinth of tools for quality control.

You need to provide it the FASTA or FASTQ (.fq) file and it will give you multiple options to control the quality of your final data. Also, it has an option to convert all letters to UPPERCASE or LOWERCASE. Output formats are also very broad. FASTA, FASTQ, FASTA + QUAL.

If you are comfortable with Linux then install the web version of PRINSEQ locally on your system and use it in your own browser. Or you can also use it completely in the commandline. If you can't do any of these then use their server. You will need to upload the whole data to the server in that case.

Here's the link to the software:

http://prinseq.sourceforge.net/

Here's the link to the online PRINSEQ server:

http://edwards.sdsu.edu/cgi-bin/prinseq/prinseq.cgi

I myself use PRINSEQ for all Quality Control tasks and it performs quite flawlessly.

Good Luck!

one liner:

samtools mpileup -uf $INFA $INBAM | bcftools view -cg - | vcfutils.pl vcf2fq | awk '/^@'${CHROM_NAME}'.*$/,/^+$/' | perl -pe "s/@/>/ ; s/\+//"

Hi Orazio,

how did you get the consensus from the file bam? my bam file come out from a dna pool of 30 individuals. Thanks

Marco

I have 6 wheat samples sequenced using RNA-seq. I received forward and reverse fastq files and I generated bam files by using hisat2 tool which are aligned with the reference wheat genome. I have been asked to build multiple sequence alignment for 3 genes from this sequenced rna seq data. I believe I need to select a gene sequence from all the samples and do a multiple sequence alignment. But I am struck in fetching the gene sequence from bam files. How do I select one gene sequence for all the 6 samples? Any suggestions? kidly send me commands for fetching the sequence of genes from this RNA seq data in fastq file or aligned.bam file.