Extraordinary Command Line: Basic Data Editing Tools for Biologists Dealing with Sequence Data

3.1. FASTA Format

FASTA is one of the most popular and one of the simplest formats in bioinformatics. FASTA is used to store DNA, RNA or protein sequence information (we use DNA information here) as plain text. Sequences can be of any length – from a few bases to whole chromosomes. A typical use of the FASTA format is to contain a whole genome assembly. Chromosomes are listed consecutively. Each chromosome is listed first with a name followed by one or more lines with sequence information. Each sequence starts on a line by itself with the greater-than character (“>”) followed by the identification number or name and description of the sequence. The next line is the first line with the actual nucleotide or amino-acid sequence, with each letter representing one nucleotide or one amino acid. Additional lines of sequence usually follow. The typical extension of a file in FASTA format is “.fasta” or “.fa” [13]. In contrast to Windows, command line applications under Linux do not rely on file extensions; however, the use of extensions is still common in order to identify the format of the file content. An example of a nucleotide sequence belonging to Ursus arctos (GenBank ID: AM411403.1) [14] is shown below (Case 1). The presented file has only one sequence, a multi-sequence FASTA file is created in Section 4.1 (Case 12).

>AM411403.1 Ursus arctos mitochondrial D-loop, isolate TAK4

ACTACTATTTTACTCCATGTCCTATTCATTTCATATATACCATTCTATGTACTGTACTATCACAGTATGT

CCTCGAATACTTTCCCCCCCCTATGTATATCGTGCATTAATGGCGTGCCCCATGCATATAAGCATGTACA

TACTGTGCTTGGTCTTACATGAGGACCTGCATTTTAGAAGTTTATCTCAGGTGTATAGTCTGCAAGCATG

TATTTCACTTAGTCCGGGAGCTTAATCACCAGGCCTCGAGAAACCATCAATCCTTGCGAGT

3.2. FASTQ Format

FASTQ [15] is the most popular format used to store high throughput sequencing data and the corresponding nucleotides qualities. This format contains information about short sequences (called reads) and their features. Each read is represented by four lines. The first line begins with a “@” character followed by information about the read (identifier, sequencing process features). The second line contains the actual nucleotide or protein sequence. The third line begins with a “+” optionally followed by information duplicated from the first line. The fourth line contains information on the phred-scaled quality score of each base in the read, encoded as a single ASCII character. The quality is expressed as the probability of incorrectly called base. The typical extension of a FASTQ file is “.fq” or “.fastq”. Moreover, due to the huge size of sequencing data, FASTQ files are often compressed using the gzip or bzip2 tools, such that the final extension is “.fq.gz”, “.fastq.gz”, “.fq.bz2”, or “.fastq.bz2”. An example of a single read in the FASTQ format (NCBI ID: SRR5078057) [16] is shown below (Case 2).

@SRR5078057.1 HWI-ST:6:1101:1149:1947/1

NGCGACCTGAACCTCTACAACAAGGAGTCCAAGCTGTCCTACTTCACCGA

+

#4=DDFFFHHGHHJJJJJJJJJJIJJJHIIJIJJJJJIJJJJJJJJJJJB

3.3. Variant Call Format (VCF)

Variant Call Format (VCF) [17] is a standard for storing information about DNA polymorphisms (e.g., SNPs, insertions, deletions and structural variants), and information on individuals’ genotypes. The actual list of polymorphisms is preceded by several lines forming a header starting with a “##” string. The header includes descriptions of the chromosomes and provides a definition for variables used to describe polymorphisms in the INFO fields. The header is followed by a single line starting with “#”, which defines the record layout for the following columns consisting of data lines, each with information on one polymorphism. These data lines contain eight tab-delimited, mandatory columns: chromosome or scaffold identifier (CHROM), the variant position expressed in base pairs (POS), variant ID (ID), a nucleotide(s) corresponding to the reference allele(s) (REF), a nucleotide representing alternative allele(s) (ALT), a quality score quantifying the polymorphism being a false positive (QUAL), filter status (FILTER), and additional information (INFO). The content of the last column is software dependent, but typically, it contains information about the type of polymorphism (e.g., SNP, SNV, indel, duplication). Missing values are coded by a dot character. An example of the VCF format (NCBI ID: GCA_000001405.27) is provided below (Case 3).

##fileformat=VCFv4.1

##fileDate=20180307

##source=ensembl;version=92;url=http://e92.ensembl.org/Homo_sapiens

##reference=ftp://ftp.ensembl.org/pub/release-92/fasta/Homo_sapiens/dna/

## [further header lines]

#CHROM POS ID REF ALT QUAL FILTER INFO

1 10177 rs367896724 A AC . . dbSNP_150;TSA=insertion;E_Freq;E_1000G;MA=C;MAF=0.425319;MAC=2130;EAS_AF=0.3363;EUR_AF=0.4056;AMR_AF=0.3602;SAS_AF=0.4949;AFR_AF=0.4909

1 10235 rs540431307 T TA . . dbSNP_150;TSA=insertion;E_Freq;E_1000G;MA=A;MAF=0.00119808;MAC=6;EAS_AF=0;EUR_AF=0;AMR_AF=0.0014;SAS_AF=0.0051;AFR_AF=0

[further body file lines]

3.4. Variant Effect Predictor Software Input and Output Formats

The Variant Effect Predictor (VEP) tool predicts the functional consequences of a polymorphism [18, 19]. The default input format is a whitespace-separated file containing the following columns: chromosome, variant start position and variant end position, both expressed in base pairs, variant polymorphism and optional information such as the strand and the name of the variant. The example below (Case 4) shows an indel variant (line 1), SNPs (lines 2-7), and a structural polymorphism (line 8). The example comes from the official VEP manual website (www.ensembl.org/info/website/upload/var.html).

1 881907 881906 -/C +

5 140532 140532 T/C +

12 1017956 1017956 T/A +

2 946507 946507 G/C +

14 19584687 19584687 C/T -

19 66520 66520 G/A + var1

8 150029 150029 A/T + var2

1 160283 471362 DUP

The default output format consists of a header starting with a “#” character. The body of the file contains 14 tab-delimited columns including: the variant polymorphism from the input, variant position, variant allele used to assign the consequence, gene ID containing or closest to the variant, transcript ID, Sequence Ontology [20] term assigned to the variant, and the predicted biochemical consequence of the variant. Missing values are indicated by '.'. The example below (Case 5) presents the first three lines of a VEP output file. It contains the header (starting with #) and genomic annotation for the structural duplication variant. Based on the presented fragment, the duplication is located close to the gene ENSG00000222623 (“upstream_gene_variant”) and overlaps with the gene ENSG00000236601 (“non_coding_transcript_exon_variant,intron_variant”). The output was generated by VEP based on the input shown in the previous example (Case 4).

#Uploaded_variation Location Allele Gene Feature Feature_type Consequence cDNA_position CDS_position Protein_position Amino_acids Codons Existing_variation Extra

1_160283_duplication 1:160282-471362 duplication ENSG00000222623 ENST00000410691 Transcript upstream_gene_variant - - - - - - IMPACT=MODIFIER;SYMBOL=RNU6-1100P;BIOTYPE=snRNA;DISTANCE=2396;STRAND=-1;SYMBOL_SOURCE=HGNC;HGNC_ID=HGNC:48063

1_160283_duplication 1:160282-471362 duplication ENSG00000236601 ENST00000412666 Transcript non_coding_transcript_exon_variant,intron_variant - - - - - - IMPACT=MODIFIER;SYMBOL=AL732372.1;BIOTYPE=lincRNA;EXON=1-2/2;INTRON=1/1;STRAND=1;SYMBOL_SOURCE=Clone_based_ensembl_gene;TSL=1

[further body lines]

4.1. Example 1: Merging Multiple Files Into One (The FASTA Format)

A reference genome is the complete sequence of an official genome of a given species. FASTA files containing a reference genome can be downloaded from publicly available databases, e.g., NCBI (www.ncbi.nlm.nih.gov) or Ensembl (www.ensembl.org). Often separate chromosomes are stored in separate files, each containing unmasked, masked and soft-masked genomic DNA sequences for a chromosome or a region not (yet) assigned to any chromosome. Although the files can be downloaded using a web browser, it is usually more convenient and faster to use the wget command followed by the web address of the file we want to download. An example of using the wget command can be found below (Case 7). The command downloads a FASTA file for mouse chromosome 1 from the Ensembl database. The output “Resolving…” shows that wget is trying to connect to the remote server and about to start the download. The progress bar keeps informing about the status of the process. The last output line indicates that the file is saved.

username@hostname:~$ wget ftp://ftp.ensembl.org /pub/release-

94/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.1.fa.gz

Resolving ftp.ensembl.org (ftp.ensembl.org)... 193.62. 193.139

Connecting to ftp.ensembl.org (ftp.ensembl.org)|193. 62.193.139|:21... connected.

Logging in as anonymous ... Logged in!

==> SYST ... done. ==> PWD ... done.

==> TYPE I ... done. ==> CWD (1) /pub/release-94/fasta/mus_musculus/dna ... done.

==> SIZE Mus_musculus.GRCm38. dna.chromosome.1.fa.gz ... 58345453

==> PASV ... done. ==> RETR Mus_musculus.GRCm38.dna.chromosome.1.fa.gz ... done.

Length: 58345453 (56M) (unauthoritative)

100%[=============================================================================================>] 58,345,453 47.9MB/s in 1.2s

2020-10-27 10:12:08 (47.9 MB/s) - ‘Mus_musculus.GRCm38.dna.chromosome.1.fa.gz’ saved [58345453]

The wget command may be modified to save more files (Case 8). The example below downloads all unmasked chromosomes, excluding other genomic elements (e.g., mitochondrial DNA). The chromosome names are provided in curly braces ({}), while [1-9XY] matches any one character between the brackets ([]). The sequence “1-9” matches any one digit. Thus, the command downloads autosomes from 1 to 9 and sex chromosomes X and Y. “?” is a wild card. It matches any single character. Therefore, all autosomes whose names begin with “1” (e.g., 13), followed by any single character, are retrieved. Here, the reference genome files in the FASTA format are named according to the following pattern: (Species).(reference_genome_version).(molecule_type).chromosome.(chromosome name).fa.gz. Thus, the file Mus_musculus.GRCm38.dna.chromosome.1.fa.gz represents the DNA sequence of the entire first chromosome of the house mouse corresponding to the GRCm38 version of the reference genome.

username@hostname:~$wget ftp://ftp.ensembl.org/pub/release-94/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.{[1-

9XY],1?}.fa.gz

The file names end with .gz. which indicates that the files are compressed. After downloading, the gunzip command must be used to uncompress the files (Case 9).

username@hostname:~$ gunzip Mus_musculus.GRCm38.dna.chromosome.*.fa.gz

The asterisk (“*”) is a wildcard, which matches zero or more occurrences of any characters. All files with names starting with “Mus_musculus.GRCm38.dna.chromosome.” and ending with “.fa.gz” will be uncompressed. Tools that allow processing and viewing of compressed files exist (e.g., the zcat command). Listing of the content of the current directory is done by the ls command (Case 10).

username@hostname:~$ ls

Mus_musculus.GRCm38.dna.chromosome.10.fa Mus_musculus.GRCm38.dna.chromosome.2.fa

Mus_musculus.GRCm38.dna.chromosome.11.fa Mus_musculus.GRCm38.dna.chromosome.3.fa

Mus_musculus.GRCm38.dna.chromosome.12.fa Mus_musculus.GRCm38.dna.chromosome.4.fa

Mus_musculus.GRCm38.dna.chromosome.13.fa Mus_musculus.GRCm38.dna.chromosome.5.fa

Mus_musculus.GRCm38.dna.chromosome.14.fa Mus_musculus.GRCm38.dna.chromosome.6.fa

Mus_musculus.GRCm38.dna.chromosome.15.fa Mus_musculus.GRCm38.dna.chromosome.7.fa

Mus_musculus.GRCm38.dna.chromosome.16.fa Mus_musculus.GRCm38.dna.chromosome.8.fa

Mus_musculus.GRCm38.dna.chromosome.17.fa Mus_musculus.GRCm38.dna.chromosome.9.fa

Mus_musculus.GRCm38.dna.chromosome.18.fa Mus_musculus.GRCm38.dna.chromosome.X.fa

Mus_musculus.GRCm38.dna.chromosome.19.fa Mus_musculus.GRCm38.dna.chromosome.Y.fa

Mus_musculus.GRCm38.dna.chromosome.1.fa

If the genome sequence of all chromosomes is required as a single file, the chromosome-specific (as is typically the case for whole genome sequence analysis) FASTA files can be merged using the cat command, where the asterisk (“*”) denotes “any characters” (Case 11).

username@hostname:~$ cat Mus_musculus.GRCm38.dna.chromosome.*fa > M.musculus.GRCm38.ref.fa

By default, the merged files would be written to the terminal. Given the size of a typical genome, this is not desirable. Instead, the “>” character redirects the output to a file, Mus_musculus.GRCm38.dna.ref.fa merging all files with names starting by “Mus_musculus.GRCm38.dna.chromosome. ” and ending with the “.fa” extension. By default, the cat command concatenates files in alphabetical order, so particular chromosomes will be ordered as: 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1, 2, 3, 4, 5, 6, 7, 8, 9, X, Y. In order to arrange the files in chromosome order, a for loop can be used. An example is provided below (Case 12).

username@hostname:~$ for chr in {1..19} X Y; do cat Mus_musculus.GRCm38.dna.chromosome.$chr.fa done > M.musculus.GRCm38.ref.num.order.fa

The first command says: in the first iteration of the loop, display the content of the Mus_musculus.GRCm38.dna. chromosome.1.fa file, in the next iteration, display the content of the file Mus_musculus.GRCm38.dna. chromosome.2.fa, and so on. As the files are displayed in the terminal, save them in the new file called M.musculus.GRCm38.ref.num.order.fa. The newly created file is in multi-sequence FASTA format.

4.2. Example 2: File Content Searching (The FASTA Format)

The regular expression is a very useful tool for pattern matching in a sequence of characters. The simplest example may be a string of characters such as “bioinfo” which would match, e.g., “bioinfo” or “bioinformatics”, but not “biotechnics”, “Bioinfo” (because of the case sensitivity) or anything else that does not contain “bioinfo” string. The grep command searches files for lines containing a pattern of interest. It may be, for example, used to check the chromosomes order by displaying only those lines of M.musculus.GRCm38.ref.num.order.fa which contain a greater than (“>”) character, which corresponds to lines containing sequence description (only the first five lines are shown here). Please note, that special characters such as a greater-than character (“>”) have to be given in quotes (Case 13). Otherwise, it indicates a redirection of output as in the previous example (Case 12).

username@hostname:~$ grep ’>’ M.musculus.GRCm38. ref.num.order.fa

>1 dna:chromosome chromosome:GRCm38:1:1:195471 971:1 REF

>2 dna:chromosome chromosome:GRCm38:2:1:182113 224:1 REF

>3 dna:chromosome chromosome:GRCm38:3:1:160039 680:1 REF

>4 dna:chromosome chromosome:GRCm38:4:1:156508 116:1 REF

>5 dna:chromosome chromosome:GRCm38:5:1:151834 684:1 REF

4.3. Example 3: File Reformatting By Column Extraction (The VCF Format)

The file used in this example was downloaded from the Ensembl database (ftp://ftp.ensembl.org/pub/release-94/variation/vcf/homo_sapiens/). It contains all known short DNA polymorphisms in Homo sapiens. In the example, we use the first 100,000 lines of this file. Different programmes for data analysis require different inputs, so reformatting the original file is an essential skill. Assume we need some information about SNPs kept in the VCF file (Case 14).

username@hostname:~$ grep -v “^#” 1000GENOMES-phase_3.100000.vcf | awk '$5 !~/,/ && $8 ~/SNV/ {print $1,$2,$2,$4”/”$5}' > 1000GENOMES-phase_3.100000.edited.txt

username@hostname:~$ head -n 5 1000GENOMES-phase_3.100000.edited.txt

1 10505 10505 A/T

1 10506 10506 C/G

1 10511 10511 G/A

1 10539 10539 C/A

1 10542 10542 C/T

The pipe character (|) directs the output of a command to the left of the pipe to become the input of the command to the right of the pipe. In the first command, the file header is removed, by ignoring all lines, which start with a “#”. It is done using the grep -v command. The ^ character denotes the beginning of a line, so the pattern “^#” only matches lines starting with a “#” character. The output of the grep command is not displayed; instead, it forms the input for the next command awk. Next, only lines containing polymorphisms classified as “SNV” (meaning Single Nucleotide Variant) are kept by awk, which checks if column 8 contains the “SNV” pattern ($8 ~/SNV/). Moreover, awk excludes (!~) those lines where column 5 includes the comma character being a separator of alternative alleles ($5 !~/,/). In practice, SNVs with more than one allele are excluded. Columns containing information about chromosome location (column 1), SNV position (column 2), and genotype (columns 4 and 5) are extracted by awk and saved into a new file named 1000GENOMES-phase_3.100000.edited.txt. This file contains two copies (as start and end positions) of column 2; columns are separated by a space character (by default, awk separates fields by a comma). The alleles are separated by a “/” character. The newly created file is in one of the input formats accepted by Variant Effect Predictor (VEP) software. In order to check the file content, five lines (-n flag) of the newly created file are displayed on the terminal by the head command.

4.4. Example 4: Information Extraction (The VEP Output Format)

From the file in the VEP output format, the list of genes (column 6) containing polymorphisms can easily be extracted by awk. Nevertheless, since one polymorphism may be located in multiple transcripts of the same gene, gene names may occur multiple times. To select unique gene names, the sort and uniq commands are used. In the following example (Case 15), gene names are sorted alphabetically, then unique gene names are saved to a file. The uniq command removes all but one copy of adjacent identical lines. To remove all duplicates, the file, therefore, needs to be sorted. This is done using the sort command.

username@hostname:~$ grep -v “^#” 1000GENOMES-phase_3.100000.edited.vep | awk '{print $4}' | grep -v “-” | sort | uniq > 1000GENOMES-phase_3.100000.genes.vep

username@hostname:~$ head -n 5 1000GENOMES-phase_3.100000.genes.vep

ENSG00000008128

ENSG00000008130

ENSG00000067606

ENSG00000078369

ENSG00000078808

The wc –l command counts the number of lines in the input. This can be used to count the number of genes containing polymorphisms without saving them to a file first. As the following example (Case 16) shows, the newly created list of genes is not saved to a file, but rather piped to the command wc -l. This counts the number of lines (genes). It works because, in this particular dataset, each line corresponds to one gene due to the use of sort | uniq.

username@hostname:~$ grep -v “^#” 1000GENOMES-phase_3.100000.edited.vep | awk '{print $4}' | grep -v “-” | sort | uniq | wc -l

4.5. Example 5: Data filtering (the VEP output format)

In the following example (Case 17), only genes containing polymorphisms, causing a frameshift in translation (column 7) are processed. The awk command displays the name of the gene (column 4) with the consequence assigned to “frameshift_variant” (column 7). The sort and uniq commands remove duplicate gene names from the list.

username@hostname:~$ awk '$7~”frameshift_variant” {print $4}' 1000GENOMES-phase_3.100000.edited.vep | sort | uniq

4.6. Example 6: Identifying Variants Shared Between Two VCF Files

Assume we would like to find SNPs common to two separate VCF files. Columns defining SNPs are 1 (CHR), 2 (POS), 4 (REF) and 5 (ALT). All of these have to be taken into account in the comparison (Case 18).

username@hostname:~$ comm -12 <(grep -v '^#' 1000GENOMES-phase_3.100000.vcf | awk '{print $1,$2,$4,$5}' | sort) <(grep -v '^#' investigated_sample.vcf | awk '{print $1,$2,$4,$5}' | sort)

1 1000018 G A

1 1000079 A G

1 1000112 G T

1 1000170 C T

1 1000242 C A

For each file, the header (lines starting with a “#”) is removed by grep –v. Next, columns are extracted by the awk command. The command comm with flags -12 prints only lines present in both files “1000GENOMES-phase_3.100000.vcf” and “investigated_sample.vcf”. Only the five first lines of output are shown.

4.7. Example 7: Counting The Number Of Reads And Their Average Length In A FASTQ File

The command below (Case 19) counts the total number of reads as well as their average length.

username@hostname:~$ awk 'NR%4==2 {reads++; readlen += length($0)} END {print reads; print readlen/reads}' SRR5078057_1.fastq

Since a sequence is stored in the second of four lines corresponding to each read, it may be extracted using modulo operator ‘%’ (remainder of integer division). ‘NR’ is the line number of the current line read. By the condition ‘NR%4==2’, we extract only lines whose remainder after division by 4 (the number of lines describing a single read) is equal to 2 (i.e., the line which contains a sequence line). Then, the number of reads is incremented by reads++. The lengths of reads are summed (readlen += length($0)). The last thing awk does is execute the commands following the ‘END’ condition. Here, awk prints the number of reads (print reads) and the average read length (print readlen/reads) in a given FASTQ file. Note that for compressed files, you must first uncompress them or read them using zcat or bzcat – depending on the type of compression (Case 20).

zcat SRR5078057_1.fastq.gz | awk 'NR%4==2 {reads++; readlen += length($0)} END {print reads; print readlen/reads}'

4.8. Example 8: Conversion Of The FASTQ Format To The FASTA Format

Sed is a stream editor used to transform text. In this example, it is used to convert FASTQ to FASTA file by reducing the number of lines corresponding to each nucleotide sequence. Normally, sed prints the input lines without any modification. With the ‘-n’ option, sed prints only what is defined by the ‘p’ option. In the example below, the 1~4 pattern selects the first line out of 4. In this line, the initial ‘@’ character is replaced by a ‘>’, since ‘s’ before the pattern indicates substitution. After the semicolon, the 2~4p pattern indicates that every fourth line starting from the second line (i.e., the line containing the sequence) should be printed. The effect is written in the new file called SRR5078057_1.fasta (Case 21).

username@hostname:~$ sed -n '1~4s/^@/>/p; 2~4p' SRR5078057_1.fastq > SRR5078057_1.fasta

username@hostname:~$ head -n 6 SRR5078057_1.fasta

>SRR5078057.1 HWI-ST:6:1101:1149:1947/1

NGCGACCTGAACCTCTACAACAAGGAGTCCAAGCTGTCCTACTTCACCGA

>SRR5078057.2 HWI-ST:6:1101:2495:1939/1

NGGCCGTCGGACTGCTCTGTGTATCAGCAGCGCTGCTGGTGCGACAGCGG

>SRR5078057.3 HWI-ST:6:1101:2340:1969/1

TGCGCTTACAAACTAATTAATAAATTAATAGTTAGCTTAAAAAGAGGCTT

4.9. Example 9: Transformation Of Columns Using A “Key File”

A common challenge is that one has to update a column in a text file with information from another source. Awk can be used to transform values in a given column using information from a “key file”, which contains pairs of old and new values. In the following example, the genes.txt file is the “key file” containing the list of Ensembl gene IDs with corresponding genes symbols. The file polymorphisms.txt contains SNPs located in some of these genes, but only Ensembl gene ID is provided. Awk transforms Ensembl gene IDs to gene symbols in polymorphism_in_genes.txt file (Case 22).

username@hostname:~$ cat genes.txt

ENSG00000227232 FAM39F

ENSG00000278267 hsa-mir-6859-1

ENSG00000284332 hsa-mir-1302-2

ENSG00000284332 MIRN1302-2

ENSG00000237613 F379

username@hostname:~$ cat polymorphism_in_genes.txt

1 14404 14404 G/A ENSG00000227232

1 34570 34570 C/G ENSG00000237613

username@hostname:~$ awk 'BEGIN{OFS=”\t”}; NR==FNR{k[$1]=$2;next}; {$5=k[$5]; print}' genes.txt polymorphism_in_genes.txt

The tab character is defined to be an output field separator by the OFS=”\t” option. The syntax NR==FNR{k[$1]=$2;next} concerns the “key file” (genes.txt) and creates the array (k) containing gene symbols from the second column ($2) indexed by Ensembl gene IDs from the first column ($1). Therefore, the number of array elements is the same as the number of lines in the genes.txt file. The syntax {$5=k[$5]; print} assigns the values from array k indexed by column 5 of polymorphism_in_genes.txt to this column and prints the file with new values (gene symbols) instead of Ensembl IDs.

6.1. Example 10: The Sequence Length And Its GC Content Calculations (FASTX)

This command computes the sequence length and its GC content. This is useful in the analysis of sequence complexity (Case 24).

username@hostname:~$ bioawk -c fastx '{print “>” $name $comment; print “Length:”, length($seq); print “GC%:”, gc($seq)}' SRR5078057_1.fasta

The -c option ‘fastx’ file format is defined. Then the ‘>’ character is printed together with the sequence ID (‘$name’) and other information (‘$comment’) stored in the first line of each sequence header. Then, the length and GC content percent for each sequence are respectively calculated using ‘length’ and ‘gc’ functions. Finally, the input file in FASTA format is provided.

6.2. Example 11: Calculating the Average Quality and the Proportion of Nucleotides Satisfying a Quality Threshold (FASTQ)

Quality control is a standard step of high-throughput sequence data analysis. The following command combination (Case 25) computes the mean sequencing quality across all reads in a FASTQ file as well as the proportion of nucleotides with a quality above a certain threshold (in this case, set to 30).

username@hostname:~$ bioawk -c fastx '{threshold+=qualcount($qual,30); seqLen+=length($seq)} END {print meanqual($qual), threshold/seqLen}' SRR5078057_1.fastq