More Unix data tools:
regular expressions, awk, sed & more

Optional self-study content

Author
Affiliations

Jelmer Poelstra

CFAES Bioinformatics Core, Ohio State University

Dept. of Plant Pathology, Ohio State University

Published on

October 10, 2026

Page last modified

October 13, 2026

WarningThis page is a draft and is still under construction.

1 Introduction

This page contains optional self-study material that builds on the w4c lecture, where you learned the basics of grep, cut, sort, and uniq. Here, you’ll go further with those commands and learn several new ones:

  • Regular expressions, to search for text patterns rather than fixed text
  • More sort and uniq, plus the related commands column, paste, and comm
  • tr to substitute or delete single characters
  • awk to filter and process tabular data
  • sed to find and replace text, and more
  • join to merge two tables, in a final example that combines many of these tools

Throughout, you’ll work with the files from the Garrigós et al. (2025) dataset, plus the Culex pipiens genome assembly FASTA file.

TipAim for reading fluency, not memorization

awk and sed are small programming languages in their own right, and nobody remembers all of their syntax. Aim to understand what the examples do, so that you can recognize, adapt, and look up such commands later. This also matters for week 8: generative AI is quite good at writing one-liners like the ones on this page, but you’ll need to be able to read them to check that they do what you want.

1.1 Setting up

  1. At https://ondemand.osc.edu, start a VS Code session in /fs/ess/PAS3493/people/$USER.

  2. Create a dir for this page and move into it:

    mkdir -p week04/data-tools
cd week04/data-tools

  3. Copy the metadata file and one FASTQ file, like in w4c:

    cp ../../garrigos/meta/metadata.tsv .
cp ../../garrigos/fastq/S63_R1.fastq.gz .

  4. Create a decompressed copy of the GTF file. This way, you won’t have to start every command with zcat:

    zcat ../../garrigos/ref/annot.gtf.gz > annot.gtf

  5. Download the genome assembly FASTA file from NCBI with wget, create a decompressed copy named genome.fna, and remove the compressed file:

    wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/801/865/GCF_016801865.2_TS_CPP_V2/GCF_016801865.2_TS_CPP_V2_genomic.fna.gz
zcat GCF_016801865.2_TS_CPP_V2_genomic.fna.gz > genome.fna
rm GCF_016801865.2_TS_CPP_V2_genomic.fna.gz

    wget downloads a file from a URL. You’ll learn more about it, and download this same file again, in week 11.

  6. Check that you have these four files:

    ls -lh
    -rw-rw----+ 1 jelmer PAS0471 123M Sep 21 10:12 annot.gtf
-rw-rw----+ 1 jelmer PAS0471 547M Sep 21 10:13 genome.fna
-rw-rw----+ 1 jelmer PAS0471  457 Sep 21 10:11 metadata.tsv
-rw-rw----+ 1 jelmer PAS0471  21M Sep 21 10:11 S63_R1.fastq.gz

The decompressed genome FASTA takes up quite a bit of space, so remove genome.fna once you’re done with this page.

2 Regular expressions with grep

2.1 What are regular expressions?

In w4c, you used grep to search for fixed text, like relictum. But grep actually interprets its search pattern as a regular expression (“regex”): a pattern in which some characters have a special meaning. With regular expressions, you can for example search for:

  • Lines that start with a > (the header lines of a FASTA file)
  • Reads that end in a stretch of 15 or more As
  • Any text between quotes that follows gene_id in a GTF file

Regular expressions are used far beyond grep: sed and awk (later on this page) use them, as do R and Python, and even VS Code’s search box (click the .* button in the search box to turn on regex search). So, time spent learning them pays off in many places.

WarningUse single quotes around regular expressions

Some regex symbols, like $ and \, also have a special meaning to the shell. Inside double quotes, the shell may therefore change the pattern before grep gets to see it. Inside single quotes, the shell leaves everything alone, so on this page, we’ll use single quotes around regular expressions.

2.2 Anchors: ^ and $

By default, a pattern can match anywhere in a line. ^ anchors a pattern to the start of a line, and $ anchors it to the end of a line. To see the difference, search the sample IDs in the metadata file (column 1) for the number 6:

cut -f 1 metadata.tsv | grep '6'
S76
S86
S64
S67
S66
S69
S63
S65
S68
# Only sample IDs that start with 'S6':
cut -f 1 metadata.tsv | grep '^S6'
S64
S67
S66
S69
S63
S65
S68
# Only sample IDs that end with '6':
cut -f 1 metadata.tsv | grep '6$'
S76
S86
S66

Anchors are especially useful with genomics files, in which certain types of lines start with a specific character. For example, to count the sequences in a FASTA file, count its header lines, which start with a >:

grep -c '^>' genome.fna
290

Let’s also look at the first few headers:

grep '^>' genome.fna | head -n 5
>NC_068937.1 Culex pipiens pallens isolate TS chromosome 1, TS_CPP_V2, whole genome shotgun sequence
>NC_068938.1 Culex pipiens pallens isolate TS chromosome 2, TS_CPP_V2, whole genome shotgun sequence
>NC_068939.1 Culex pipiens pallens isolate TS chromosome 3, TS_CPP_V2, whole genome shotgun sequence
>NW_026292818.1 Culex pipiens pallens isolate TS unplaced genomic scaffold, TS_CPP_V2 Cpp_Un0001, whole genome shotgun sequence
>NW_026292819.1 Culex pipiens pallens isolate TS unplaced genomic scaffold, TS_CPP_V2 Cpp_Un0002, whole genome shotgun sequence

So, the genome assembly consists of 290 sequences: the three Culex chromosomes, and many smaller “unplaced scaffolds” (pieces of the genome that could not be assigned to a chromosome).

Combining both anchors, ^$ matches empty lines, so grep -v '^$' removes all empty lines from its input.

2.3 Any character (.), repetition (*), and grep -o

Two more key regex symbols are:

  • . — matches any single character.
  • * — matches the preceding character any number of times, including zero.

Together, .* therefore matches any number of any characters, like the * shell wildcard does.

To see exactly what a regular expression matches, use grep’s -o option, which prints only the matching part of each line instead of the full line. Let’s try to extract the gene ID (the quoted text after gene_id) from the first line of the GTF file’s table part:

grep -v '^#' annot.gtf | head -n 1 | grep -o 'gene_id ".*"'
gene_id "LOC120427725"; transcript_id ""; db_xref "GeneID:120427725"; description "homeotic protein deformed"; gbkey "Gene"; gene "LOC120427725"; gene_biotype "protein_coding"

That’s a lot more than we wanted! The problem is that .* is “greedy”: it matches as much as possible, so the match runs all the way to the last " on the line. Character classes, discussed next, offer a solution.

To match one of these special characters literally, put a backslash before it (“escape” it). For example, NC_068937\.1 only matches an actual period between 7 and 1, while NC_068937.1 would also match, e.g., NC_06893701.

2.4 Character classes: [...] and [^...]

A character class matches one character from a set, which you list between square brackets:

  • [ACGT] matches a single A, C, G, or T.
  • Ranges are allowed: [0-9] matches any digit, and [a-z] any lowercase letter.
  • A ^ as the first character inside the brackets negates the class: [^ACGT] matches any single character that is not A, C, G, or T.

Note the two different meanings of ^: outside brackets, it anchors to the start of a line, while as the first character inside brackets, it negates the character class.

Negated character classes are the solution to our greedy matching problem. Instead of “any characters, followed by a quote” (.*"), we ask for “any characters except a quote, followed by a quote” ([^"]*"). This way, the match can’t extend beyond the next quote:

grep -v '^#' annot.gtf | head -n 1 | grep -o 'gene_id "[^"]*"'
gene_id "LOC120427725"

The same trick works to extract just the sequence IDs from the FASTA headers — everything after the > up until the first space:

grep -o '^>[^ ]*' genome.fna | head -n 4
>NC_068937.1
>NC_068938.1
>NC_068939.1
>NW_026292818.1

Extracting and counting matches

Combining grep -o with the sort | uniq -c frequency-table approach from w4c is a very handy way to summarize information that is “hidden” inside a column. For example, the GTF file’s attribute column (column 9) has a gene_biotype for every gene. How many genes of each type are there?

grep -o 'gene_biotype "[^"]*"' annot.gtf | sort | uniq -c | sort -nr
  16297 gene_biotype "protein_coding"
   1620 gene_biotype "lncRNA"
    842 gene_biotype "pseudogene"
    686 gene_biotype "tRNA"
    155 gene_biotype "rRNA"
     58 gene_biotype "snRNA"
      9 gene_biotype "snoRNA"
      5 gene_biotype "transcribed_pseudogene"
      1 gene_biotype "misc_RNA"

Most genes are protein-coding, but the annotation also includes, e.g., about 1,600 long non-coding RNA (lncRNA) genes.

2.5 Extended regular expressions: grep -E

Several other useful regex symbols are only available as “extended regular expressions”, which you can turn on with grep -E:

Symbol Matches Example Example matches
+ The preceding character one or more times [0-9]+ 7, 2046
? The preceding character zero or one time genes? gene, genes
{n} The preceding character exactly n times A{3} AAA
{n,} The preceding character at least n times A{3,} AAA, AAAAAA
{n,m} The preceding character n to m times A{3,5} AAA, AAAAA
| Either the pattern on the left or the one on the right cecropin|defensin cecropin, defensin
( ) Groups a pattern, e.g. to repeat it (CA)+ CA, CACACA

For example, how many reads in the FASTQ file end in at least 15 As, which may be the poly-A tails of mRNAs?

zcat S63_R1.fastq.gz | grep -cE 'A{15,}$'
199

And which GTF features belong to genes for the antimicrobial peptides cecropin or defensin, which are part of the mosquito immune response?

grep -E 'cecropin|defensin' annot.gtf | cut -f 3 | sort | uniq -c
     16 CDS
     17 exon
      7 gene
      8 start_codon
      8 stop_codon
      8 transcript

Without -E, grep uses “basic” regular expressions, in which +, ?, { }, |, and ( ) are just literal characters. In the GNU version of grep (which is what OSC uses), you can still use their special meaning by putting a backslash before them:

Extended (grep -E) Basic (grep)
A+ A\+
A{15,} A\{15,\}
cecropin|defensin cecropin|defensin

That gets hard to read quickly, so it’s easiest to simply use -E whenever you need these symbols. The same goes for sed, which also has an -E option. awk always uses extended regular expressions.

When searching online, you’ll also come across “Perl-compatible” regular expressions (grep -P), which offer additional shortcuts. For example, \t matches a Tab with grep -P, but not with grep or grep -E:

# Print lines with 'gene' as the feature type (column 3), surrounded by Tabs:
grep -P '\tgene\t' annot.gtf | head -n 2

When you search a FASTQ file for a sequence, grep only prints the matching sequence line. To see the full 4-line FASTQ record of each matching read, also print 1 line before (-B 1) and 2 lines after (-A 2) each match. The --no-group-separator option prevents grep from printing -- lines between non-adjacent matches, so that the output is valid FASTQ:

zcat S63_R1.fastq.gz | grep --no-group-separator -B 1 -A 2 -E 'A{15,}$' | head -n 8
@ERR10802863.17755922 17755922 length=73
AGAATGCATACACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@ERR10802863.21692370 21692370 length=73
AGATCAACTACTGAGAAATGTGTTGAGCGATAAATACAGTTCATTAAAATCAACCAGAAAAAAAAAAAAAAAA
+
AAAAAEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

2.6 Regular expressions vs. shell wildcards

Regular expressions and the shell wildcards from w4b (and the More Unix shell page) look similar, but they are different things:

  • Shell wildcards match file names, and are expanded by the shell before a command runs.
  • Regular expressions match text in a command’s input, and are interpreted by the command itself (e.g. grep).

Confusingly, some symbols have a different meaning in each:

Shell wildcard Regex equivalent Matches
? . Any single character
* .* Any number of any characters
[...] and [^...] Same One character that is / is not in the set

Also recall from w4c that grep patterns don’t need to match the entire line, whereas a glob pattern has to match the entire file name.

Exercise 1: Regular expressions

  1. Thinking of the shell wildcard S6*, a classmate expects the command below to count the 7 samples whose ID starts with S6. What count do you think it will report instead, and why? Run the command to check your answer.

    cut -f 1 metadata.tsv | grep -c 'S6*'

      Click to see the solution

    It reports 22: every single sample ID.

    In a regular expression, * means “the preceding character any number of times, including zero”. So, S6* matches an S followed by zero or more 6s — and every sample ID contains an S.

    The regex equivalent of the wildcard is S6.*, but since grep patterns don’t need to match the full line, plain S6 (or, more precisely, ^S6) is all you need here.

  2. Count the number of FASTA headers that start with >NC, and the number that start with >NW, and print the >NC headers. What do the NC sequences represent, and what do the NW ones represent?

      Click to see the solution

     
    grep -c '^>NC' genome.fna
    4
    grep -c '^>NW' genome.fna
    286
    grep '^>NC' genome.fna
    >NC_068937.1 Culex pipiens pallens isolate TS chromosome 1, TS_CPP_V2, whole genome shotgun sequence
>NC_068938.1 Culex pipiens pallens isolate TS chromosome 2, TS_CPP_V2, whole genome shotgun sequence
>NC_068939.1 Culex pipiens pallens isolate TS chromosome 3, TS_CPP_V2, whole genome shotgun sequence
>NC_015079.1 Culex pipiens pipiens mitochondrion, complete genome

    The 4 NC sequences are the 3 chromosomes plus the mitochondrial genome, while the 286 NW sequences are the unplaced scaffolds (in NCBI’s RefSeq database, NC_ IDs are used for complete molecules like chromosomes, and NW_ IDs for scaffolds).

  3. A labmate counts the reads in a FASTQ file by counting the header lines, which start with @:

    zcat S63_R1.fastq.gz | grep -c '^@'

    For this file, that gives the correct answer of 500,000 reads. But why could this approach overcount the reads in other FASTQ files? (Hint: think about the quality scores, and see the table in the w4a FASTQ quality score slides.)

      Click to see the solution

    The @ character is also used as a quality score: it represents a Phred score of 31 (one more than ?, which represents 30). So, a quality score line can start with an @, and would then be counted as a read.

    It happens to work for our file because its quality scores are “binned” into just a few values, none of which is 31. More robust approaches are to divide the number of lines by 4, or to search for a longer, more specific header pattern like '^@ERR'.

  4. How many genes have a description that contains “heat shock”, and are any of these descriptions shared by multiple genes? (Hint: description attributes are only present on gene lines, so you can adapt the gene_biotype command above.)

      Click to see the solution

    Allow any characters except a quote both before and after “heat shock”, so that the match stays within the description’s quotes:

    grep -o 'description "[^"]*heat shock[^"]*"' annot.gtf | sort | uniq -c | sort -nr
          3 description "heat shock protein 70 B2"
      2 description "heat shock 70 kDa protein cognate 4-like"
      1 description "heat shock protein beta-1"
      1 description "heat shock protein 83-like"
      1 description "heat shock protein 83"
      1 description "heat shock protein 75 kDa, mitochondrial"
      1 description "heat shock protein 70 A2-like"
      1 description "heat shock protein 60A-like"
      1 description "heat shock factor-binding protein 1"
      1 description "heat shock factor protein"
      1 description "heat shock 70 kDa protein cognate 5"
      1 description "heat shock 70 kDa protein cognate 4"
      1 description "heat shock 70 kDa protein cognate 2"
      1 description "heat shock 70 kDa protein 4"
      1 description "activator of 90 kDa heat shock protein ATPase homolog 1"
      1 description "10 kDa heat shock protein, mitochondrial"

    There are 16 distinct descriptions, two of which are shared by multiple genes. To count the genes, count the matching lines:

    grep -c 'description "[^"]*heat shock' annot.gtf
    19

3 More sort and uniq, and friends

3.1 Sorting by column, numerically

In w4c, you saw that sort -k sorts by specific columns, and that sort -n sorts numerically. You can combine the two by adding an n to a specific column (key).

To try this, create a small file with the sequence ID, feature type, and start and end positions (columns 1 and 3-5) of the first 6 genes in the GTF file. (The -w option makes grep only match gene as a whole word.)

cut -f 1,3-5 annot.gtf | grep -w 'gene' | head -n 6 > genes6.tsv
cat genes6.tsv
NC_068937.1     gene    2046    110808
NC_068937.1     gene    48457   50554
NC_068937.1     gene    76295   82374
NC_068937.1     gene    82671   86331
NC_068937.1     gene    160196  177583
NC_068937.1     gene    299209  303777

These genes are already sorted by their start position (column 3 in this file). But see what happens when you sort them by that column without telling sort that it contains numbers:

sort -k3,3 genes6.tsv
NC_068937.1     gene    160196  177583
NC_068937.1     gene    2046    110808
NC_068937.1     gene    299209  303777
NC_068937.1     gene    48457   50554
NC_068937.1     gene    76295   82374
NC_068937.1     gene    82671   86331

Sorting alphabetically, sort compares character by character, so 160196 comes before 2046 because 1 comes before 2. With -k3,3n, you get the correct numeric order back:

sort -k3,3n genes6.tsv
NC_068937.1     gene    2046    110808
NC_068937.1     gene    48457   50554
NC_068937.1     gene    76295   82374
NC_068937.1     gene    82671   86331
NC_068937.1     gene    160196  177583
NC_068937.1     gene    299209  303777

A common use of this is to sort an annotation file by position, first by sequence ID (alphabetically) and then by start position (numerically). For the GTF file, whose start positions are in column 4, that would be:

# [Example only, no need to run this]
sort -k1,1 -k4,4n annot.gtf

Many genomics tools require their input files to be sorted this way.

3.2 Other useful sort options

  • -c checks whether a file is already sorted. It prints nothing if the file is sorted, and reports the first out-of-order line if not:

    sort -c -k3,3n genes6.tsv
    sort -c -k3,3 genes6.tsv
    sort: genes6.tsv:5: disorder: NC_068937.1     gene    160196  177583

  • -u only keeps unique lines, so sort -u does the same as sort | uniq.

  • -V (“version sort”) sorts numbers within text by their numeric value:

    # printf prints text like echo does, and '\n' starts a new line:
printf "sample1\nsample10\nsample2\n" | sort -V
    sample1
sample2
sample10

    (Without -V, sample10 would come before sample2.)

  • -g sorts numbers in scientific notation (like 1e-5) correctly, which -n doesn’t. This is useful for, e.g., p-values.

  • -t sets the column delimiter, which by default is any whitespace. For example, sort -t , -k2,2 sorts a CSV file by its second column.

3.3 Finding duplicates with uniq -d

With the -d option, uniq only prints lines that occur more than once (one copy of each). For example, how many gene descriptions are shared by multiple genes?

grep -o 'description "[^"]*"' annot.gtf | sort | uniq -d | wc -l
1424

To see which descriptions are shared by the most genes, use the frequency table approach instead:

grep -o 'description "[^"]*"' annot.gtf | sort | uniq -c | sort -nr | head -n 5
    156 description "5S ribosomal RNA"
    141 description "uncharacterized protein K02A2.6-like"
     44 description "pickpocket protein 28-like"
     41 description "transfer RNA valine (anticodon AAC)"
     39 description "transfer RNA tyrosine (anticodon GUA)"

(The 156 5S ribosomal RNA genes are not an error: rRNA genes typically occur in many copies in a genome.)

3.4 column -t for viewing tabular files

When the values in a column differ in length, tabular files can look messy in the terminal:

grep -v '^#' annot.gtf | cut -f 1-8 | head -n 3
NC_068937.1     Gnomon  gene    2046    110808  .       +       .
NC_068937.1     Gnomon  transcript      2046    110808  .       +       .
NC_068937.1     Gnomon  exon    2046    2531    .       +       .

column -t lines up the columns nicely:

grep -v '^#' annot.gtf | cut -f 1-8 | head -n 3 | column -t
NC_068937.1  Gnomon  gene        2046  110808  .  +  .
NC_068937.1  Gnomon  transcript  2046  110808  .  +  .
NC_068937.1  Gnomon  exon        2046  2531    .  +  .

column is only meant for viewing: it separates columns with spaces, so don’t save its output as a TSV file. For a CSV file, specify the delimiter with -s: column -s , -t my-data.csv.

3.5 paste to combine lines

paste combines lines side by side, separated by Tabs. Given two files, it pastes line 1 of the first file next to line 1 of the second, and so on.

A neat trick is to give paste a dash (-) for each line it should combine from its standard input. With four dashes, paste combines every four lines into one — which turns each 4-line FASTQ record into a single line:

zcat S63_R1.fastq.gz | paste - - - - | head -n 2
@ERR10802863.8435456 8435456 length=74  CAACGAATACATCATGTTTGCGAAACTACTCCTCCTCGCCTTGGTGGGGATCAGTACTGCGTACCAGTATGAGT      +       AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@ERR10802863.27637245 27637245 length=74        GCCACACTTTTGAAGAACAGCGTCATTGTTCTTAATTTTGTCGGCAACGCCTGCACGAGCCTTCCACGTAAGTT      +       AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE<EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

Now that each read is one line, with the sequence in column 2, you can use all the column-based tools you know. For example, which read sequences occur most often?

zcat S63_R1.fastq.gz | paste - - - - | cut -f 2 | sort | uniq -c | sort -nr | head -n 5
  32589 NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    255 ATAAATTAAAGGGCCGCAGTATTTTGACTGTGCGAAGGTAGCATAATCACTAGTCTTTTAATTGGAGGCTTGT
    223 AACAGCGTAATTTTTTTTTAGAGTTCATATCGACAAAAAAGATTGCGACCTCGATGTTGGATTAAGAGTTATTT
    204 AAATTAAAGGGCCGCAGTATTTTGACTGTGCGAAGGTAGCATAATCACTAGTCTTTTAATTGGAGGCTTGTATG
    138 GCAAAAGTACGTCGACGTGCACTGCCACGACAGTGACGACCAGAGCGTGAAGCTGTACAAGAACCAGATCGCC

By far the most common “sequence” is that of the all-N reads you found in w4c. The FastQC program, which you’ll use in week 6, reports such “overrepresented sequences” too.

Exercise 2: sort, uniq, and paste

  1. The experiment has 2 time points and 3 treatments. Create a frequency table of the number of samples for each combination of time point and treatment. Is the design balanced, i.e. does every combination have the same number of samples?

      Click for a hint

    Select both the time point and the treatment column with cut.

      Click to see the solution

    Select columns 2 and 3 with cut -f 2,3, and keep the rest of the w4c frequency table pipeline the same:

    tail -n +2 metadata.tsv | cut -f 2,3 | sort | uniq -c
          4 10dpi   cathemerium
      4 10dpi   control
      4 10dpi   relictum
      3 24hpi   cathemerium
      3 24hpi   control
      4 24hpi   relictum

    Because uniq compares entire lines, it counts each distinct combination of time point and treatment. The design is not fully balanced: at 24hpi, there are only 3 samples for both the cathemerium and control treatments.

  2. A labmate wants to sort the GTF file by start position (column 4). They remove the header lines and then sort as follows:

    grep -v '^#' annot.gtf | sort -k4

    Without running it, can you spot two problems with this command?

      Click to see the solution
    1. Without n, the start positions are sorted alphabetically, so that, e.g., 160196 would come before 2046.
    2. The sequence ID (column 1) is ignored, so features from different chromosomes and scaffolds get mixed up. To sort by position within each sequence, use sort -k1,1 -k4,4n.

    (A minor third point: -k4 means “from column 4 to the end of the line”, rather than just column 4 as -k4,4 would.)

    You can see the second problem right away in the output, whose first lines come from three different scaffolds:

    grep -v '^#' annot.gtf | sort -k4 | head -n 3 | cut -f 1-5
    NW_026292908.1  Gnomon  CDS     10      11481
NW_026292919.1  Gnomon  CDS     10      1167
NW_026292909.1  Gnomon  CDS     10      1587

  3. paste - - - - is a safe way to put each FASTQ read on a single line. Would paste - - similarly be a safe way to put each sequence in genome.fna on a single line, next to its header? Why or why not?

      Click to see the solution

    No: in FASTQ files, every record consists of exactly 4 lines, but in FASTA files, a sequence can be spread across any number of lines (recall exercise 2 in the w4a slides). In genome.fna, each line holds 80 bases, so paste - - would mostly paste one 80-base line next to the next one.

4 tr to translate characters

4.1 Substituting characters

The tr (translate) command substitutes characters one-for-one: the first character in its first argument is replaced by the first character in its second argument, and so on. For example, to replace every a by a b:

echo "aaaaccc" | tr 'a' 'b'
bbbbccc

Oddly enough, tr does not accept a file name as an argument: it only reads from standard input. So, to use it on a file, pipe the output of cat into tr. For example, to turn the Tab-separated metadata.tsv into a comma-separated file, replace each Tab (\t) by a comma:

cat metadata.tsv | tr '\t' ',' > metadata.csv
head -n 3 metadata.csv
sample_id,time,treatment
S82,10dpi,cathemerium
S75,10dpi,cathemerium

tr also understands ranges of characters. This is handy because many genome FASTA files, including ours, are “soft-masked”: repetitive regions are written in lowercase letters. For example, lines 12 and 13 of the genome FASTA file contain both:

head -n 13 genome.fna | tail -n 2
CCTAATTCATGGAGGTGAATAGAGCATACGTAAATACAAAACTCATGACATTAGCCTGTAAGGATTGTGTaattaatgca
aaaatattgaTAGAATGAAAGATGCAAGTCccaaaaattttaagtaaatgaATAGTAATCATAAAGATAActgatgatga

To convert everything to uppercase, translate the range a-z into A-Z:

head -n 13 genome.fna | tail -n 2 | tr 'a-z' 'A-Z'
CCTAATTCATGGAGGTGAATAGAGCATACGTAAATACAAAACTCATGACATTAGCCTGTAAGGATTGTGTAATTAATGCA
AAAATATTGATAGAATGAAAGATGCAAGTCCCAAAAATTTTAAGTAAATGAATAGTAATCATAAAGATAACTGATGATGA

4.2 Reverse complementing a sequence

Combined with rev, which reverses each line, tr can create the reverse complement of a DNA sequence: rev takes care of the “reverse” part, and tr of the “complement” part by swapping A with T and C with G:

echo "ATGCCGTA" | rev | tr 'ACGT' 'TGCA'
TACGGCAT

4.3 Deleting and “squeezing” characters

  • The -d option deletes characters. For example, by deleting all newlines (\n) from the sequence lines of the genome FASTA, you can count the total number of bases in the genome assembly (this takes a few seconds):

    grep -v '^>' genome.fna | tr -d '\n' | wc -c
    566354144

    So, the Culex pipiens genome assembly is about 566 million bases (566 Mb) long.

  • The -s option “squeezes” runs of the same character into one. For example, to turn a variable number of spaces into a single space:

    echo "a     b   c" | tr -s ' '
    a b c

Exercise 3: tr

  1. A labmate wants to rename the control treatment to mock in the metadata file, and runs the command below. What do you think the output will look like? Run the command to check.

    cat metadata.tsv | tr 'control' 'mock'

      Click to see the solution

    The output is mangled, also outside of the control lines:

    sampke_id       kime    kkeakmeck
S82     10dpi   makhemekium
S75     10dpi   makhemekium
[...]
S78     10dpi   mkckkkk
S84     10dpi   mkckkkk
[...]

    tr translates individual characters, not words: every c became an m, every n a c, and so on, anywhere in the file. Because the second argument is shorter than the first, the remaining characters (t, r, o, and l) all became its last character, k.

    To replace a word, use sed, which is covered below.

  2. Use tr to calculate which percentage of the genome assembly is soft-masked (i.e., written in lowercase).

      Click for a hint

    Combine -d with the -c option, which takes the complement of the (first) set of characters: e.g., tr -cd 'xyz' deletes everything except x, y, and z. Then compare with the total of 566,354,144 bases.

      Click to see the solution

    Delete everything except lowercase bases, and count the remaining characters:

    grep -v '^>' genome.fna | tr -cd 'acgtn' | wc -c
    315245493

    315,245,493 / 566,354,144 = 0.56, so 56% of the assembly is soft-masked. Large fractions of many genomes, including those of mosquitoes, consist of repetitive DNA such as transposable elements.

    (The grep -v '^>' step is needed because the header lines contain lowercase letters too.)

  3. A labmate wants the reverse complement of the entire genome, and runs:

    rev genome.fna | tr 'ACGT' 'TGCA' > genome_revcomp.fna

    Without running it, can you name at least two problems with this command?

      Click to see the solution
    1. rev reverses each line separately, but leaves the order of the lines unchanged. Because the sequences are spread across many lines, they are therefore not properly reversed.
    2. The header lines are also reversed and “complemented”: e.g., the first header would end in [...] sneipip xeluG 1.739860_GN>.
    3. Lowercase (soft-masked) bases are not complemented, since tr is case-sensitive.

    This is a task for a dedicated bioinformatics tool.

5 awk for tabular data processing

awk is a small programming language in its own right, which is especially good at quickly processing tabular data. For example, it can select rows based on the values in specific columns, reorder columns, or do calculations. To understand awk, you need to know about two core concepts: records and fields, and pattern-action pairs.

5.1 Records and fields

awk processes its input one record at a time, and by default, each line is one record. Within a record, each column is a field. The entire record is available as $0, and the fields as $1 (first column), $2 (second column), and so on.

An awk command is placed between single quotes, and actions go between curly braces { }. For example, to print columns 3 and 1 of the metadata file:

awk '{ print $3, $1 }' metadata.tsv | head -n 3
treatment sample_id
cathemerium S82
cathemerium S75

This already goes beyond cut, which can’t change the order of columns. Note how the comma , between the columns made awk separate them with a space. Without the comma, the columns are pasted together, and to separate them with a Tab, you can add a "\t":

awk '{ print $3 $1 }' metadata.tsv | head -n 3
treatmentsample_id
cathemeriumS82
cathemeriumS75
awk '{ print $3 "\t" $1 }' metadata.tsv | head -n 3
treatment       sample_id
cathemerium     S82
cathemerium     S75
WarningSet the column delimiter with -F '\t' for TSV files

By default, awk splits fields on any whitespace: Tabs as well as spaces. That works for the metadata file, which contains no spaces, but not for the GTF file, whose ninth column contains spaces. See what happens when you try to print column 9 of a GTF line:

grep -v '^#' annot.gtf | head -n 1 | awk '{ print $9 }'
gene_id

The built-in variable NF (number of fields) shows that awk split this line into 24 fields, rather than the 9 columns of a GTF file:

grep -v '^#' annot.gtf | head -n 1 | awk '{ print NF }'
24

With -F '\t', you tell awk to only split fields on Tabs:

grep -v '^#' annot.gtf | head -n 1 | awk -F '\t' '{ print NF }'
9

To be safe, always use -F '\t' with Tab-separated files.

5.2 Pattern-action pairs

The general syntax of an awk command is awk 'pattern { action }' file, where:

  • The pattern is a condition that is tested for each record.
  • The action is performed for each record for which the pattern is true.

Either of the two can be omitted:

  • Omit the pattern, like in the examples above, and the action is performed for every record.
  • Omit the action, and every record that matches the pattern is printed in full.

Patterns

Patterns are often comparisons, such as $3 == "control" (column 3 equals “control”; note the double == and the quotes around the text). Patterns can also be combined, for example with && (“and”). To print all infected (i.e., non-control) samples from the 24hpi time point:

awk '$2 == "24hpi" && $3 != "control"' metadata.tsv
S64     24hpi   cathemerium
S67     24hpi   cathemerium
S70     24hpi   cathemerium
S71     24hpi   relictum
S74     24hpi   relictum
S65     24hpi   relictum
S68     24hpi   relictum

Here are the most important comparison and logical operators:

Operator Meaning
a == b a is equal to b
a != b a is not equal to b
a < b / a > b a is less than / greater than b
a <= b / a >= b a is less than or equal to / greater than or equal to b
a ~ /b/ a matches regular expression b
a !~ /b/ a does not match regular expression b
a && b a and b
a || b a or b
!a not a

As the ~ operator shows, awk can also use regular expressions, which go between forward slashes. For example, $1 ~ /^NW_/ would select features on the unplaced scaffolds.

You can also calculate in patterns. For example, the length of a GTF feature is its end position (column 5) minus its start position (column 4), plus 11. So, to print all genes longer than 1 Mb (1,000,000 bases):

awk -F '\t' '$3 == "gene" && $5 - $4 + 1 > 1000000' annot.gtf
NC_068939.1     Gnomon  gene    10362825        11542514        .       +       .       gene_id "LOC120432606"; transcript_id ""; db_xref "GeneID:120432606"; description "RNA-binding protein Musashi homolog Rbp6"; gbkey "Gene"; gene "LOC120432606"; gene_biotype "protein_coding";

Only one gene, on chromosome 3, is that long.

Actions

Patterns and actions can of course be combined. For example, to print the sequence ID, start, end, and length of each gene:

awk -F '\t' '$3 == "gene" { print $1, $4, $5, $5 - $4 + 1 }' annot.gtf | head -n 3
NC_068937.1 2046 110808 108763
NC_068937.1 48457 50554 2098
NC_068937.1 76295 82374 6080

5.3 BEGIN and END

Two special patterns are BEGIN and END, whose actions are performed once: before the first record and after the last record, respectively.

BEGIN is often used to set options. For example, setting the variable OFS (output field separator) to a Tab makes the commas in print produce Tabs instead of spaces:

awk -F '\t' 'BEGIN { OFS = "\t" } $3 == "gene" { print $1, $4, $5, $5 - $4 + 1 }' annot.gtf | head -n 3
NC_068937.1     2046    110808  108763
NC_068937.1     48457   50554   2098
NC_068937.1     76295   82374   6080

END is useful for printing a summary after processing all records. For example, to compute the mean gene length, add up the lengths of all genes in a variable sum, count the genes in a variable n, and divide the two at the end:

awk -F '\t' '$3 == "gene" { sum += $5 - $4 + 1; n++ } END { print sum / n }' annot.gtf
18961.9

In this command:

  • sum += x is shorthand for sum = sum + x, and n++ for n = n + 1.
  • Two statements within one action are separated by a semicolon ;.
  • You don’t need to create (initialize) sum and n beforehand: awk variables start out as 0.

5.4 Processing FASTQ files with awk

The built-in variable NR holds the current record (line) number. This is very handy for FASTQ files, where the sequence is always on the 2nd line of each 4-line record: the lines for which NR % 4 == 2 (the remainder after dividing the line number by 4 is 2). Combined with the length() function, this gives you the length of every read — and then, a frequency table of read lengths:

zcat S63_R1.fastq.gz | awk 'NR % 4 == 2 { print length($0) }' | sort | uniq -c | sort -nr | head -n 5
 303297 74
 123453 73
  32678 35
  28617 72
   8483 71

Most reads are 72–74 bases long, but 32,678 reads are only 35 bases long: the all-N reads you saw in w4c and in the paste section above.

5.5 awk reference tables

Some useful built-in variables:

Variable Meaning
$0 The entire record (by default, a line)
$1, $2, … The first, second, … field (by default, a column)
NR Number of the current record (in END: total number of records)
NF Number of fields in the current record
FS Input field separator (set with -F; default: any whitespace)
OFS Output field separator (default: a space)

Some useful functions and statements:

Function / statement Meaning
length(s) Number of characters in s (length() alone: in the whole record)
tolower(s) / toupper(s) s converted to lowercase / uppercase
substr(s, start, n) The part of s that starts at position start and is n characters long
split(s, arr, sep) Splits s into pieces on sep, and stores the pieces in array arr
sub(regex, repl, s) Replaces the first match of regex in s by repl
gsub(regex, repl, s) Replaces every match of regex in s by repl
next Skips to the next record, without considering any remaining pattern-action pairs
exit Stops processing records (any END action is still performed)

Exercise 4: awk

  1. A classmate wants to print the control samples, but makes a small typo:

    awk '$3 = "control"' metadata.tsv

    What do you think this will print? Run the command to check.

      Click to see the solution

    It prints every line, with the third column replaced by control (and with spaces instead of Tabs between the columns):

    sample_id time control
S82 10dpi control
S75 10dpi control
[...]

    A single = doesn’t compare but assigns a value: it sets column 3 to control for every line. Furthermore, the result of that assignment, the text control, counts as “true”, so every line is printed. When awk changes a field, it rebuilds the line using OFS, which is a space by default.

    With ==, the command works as intended:

    awk '$3 == "control"' metadata.tsv

  2. Compute the mean length of exons. Why do you think it is so much shorter than the mean gene length of nearly 19,000 bases?

      Click to see the solution

    Change the pattern from "gene" to "exon":

    awk -F '\t' '$3 == "exon" { sum += $5 - $4 + 1; n++ } END { print sum / n }' annot.gtf
    542.133

    Exons are only about 540 bases long on average, because most genes consist of multiple exons separated by introns, which are often much longer than the exons. The length of a gene spans all of its exons and introns.

  3. Remove all reads shorter than 50 bases from S63_R1.fastq.gz, and save the remaining reads in a new FASTQ file S63_R1_filtered.fastq. How many reads are left?

      Click for a hint

    Put each read on one line with paste, filter with awk, and then turn each line back into a 4-line FASTQ record with tr. (The length() function can also take a specific field, like length($2).)

      Click to see the solution

    After paste - - - -, the sequence is in column 2, so keep only lines with length($2) >= 50. Then, replace the Tabs that paste inserted by newlines:

    zcat S63_R1.fastq.gz |
    paste - - - - |
    awk -F '\t' 'length($2) >= 50' |
    tr '\t' '\n' > S63_R1_filtered.fastq

    (A line that ends with a pipe | continues on the next line, which makes long pipelines easier to read.)

    Count the reads by dividing the number of lines by 4:

    wc -l S63_R1_filtered.fastq
    1869040 S63_R1_filtered.fastq

    1,869,040 / 4 = 467,260 reads are left, so 32,740 reads were shorter than 50 bases.

    (In practice, you would use a dedicated read-trimming program for this, which can filter reads by length and quality in one go.)

6 sed to find and replace

6.1 Substitution with s/pattern/replacement/

sed (stream editor) is most often used to find and replace text, using the syntax s/pattern/replacement/, where s stands for “substitute”. Unlike tr, sed replaces whole patterns rather than single characters. So, here is how you can rename the control treatment to mock (recall Exercise 3a):

sed 's/control/mock/' metadata.tsv | head -n 8 | tail -n 4
S83     10dpi   cathemerium
S78     10dpi   mock
S84     10dpi   mock
S77     10dpi   mock

The pattern is a regular expression, and by default, sed only replaces the first match on each line. To replace all matches, add the g (“global”) modifier at the end:

echo "AAACCCAAA" | sed 's/A/T/'
TAACCCAAA
echo "AAACCCAAA" | sed 's/A/T/g'
TTTCCCTTT
WarningSaving the output of sed

Like the other commands on this page, sed does not change the input file, but prints its output to the screen. To save the output, redirect it to a new file (recall from w4c that you should never redirect to the input file itself):

sed 's/control/mock/' metadata.tsv > metadata_mock.tsv

sed also has an -i option to edit a file in place, e.g. sed -i 's/control/mock/' metadata.tsv. This is convenient, but be careful: there is no undo, so make sure your command works before adding -i.

6.2 Selecting lines with addresses

By default, sed operates on every line. You can restrict it to certain lines by putting an “address” in front of the command: either a line number or range of line numbers, or a regular expression between slashes. For example, many tools use only the first “word” of a FASTA header as the sequence ID, and long headers can be a nuisance. To remove everything from the first space onwards, but only on header lines:

sed '/^>/ s/ .*//' genome.fna | grep '^>' | head -n 3
>NC_068937.1
>NC_068938.1
>NC_068939.1

Besides s, two other handy sed commands are d (delete) and p (print):

  • sed '/^#/d' annot.gtf deletes all lines starting with #, i.e. it does the same as grep -v '^#' annot.gtf.
  • sed '1d' metadata.tsv deletes the first line, i.e. it does the same as tail -n +2 metadata.tsv.
  • sed -n '12,13p' genome.fna prints lines 12 to 13. The -n option turns off sed’s default printing of every line, so only the lines selected with p are printed.

6.3 Using a different delimiter

The / in s/pattern/replacement/ can be replaced by (almost) any other character. This is useful when the pattern or replacement contains slashes, like file paths. For example, to remove the dirs from a list of file paths (the .*/ greedily matches everything up to and including the last /):

ls ../../garrigos/fastq/*fastq.gz | sed 's#.*/##' | head -n 3
S63_R1.fastq.gz
S63_R2.fastq.gz
S64_R1.fastq.gz

To apply multiple substitutions in one sed command, put -e before each one. For example, to also remove the read direction and extension, and get a list of sample IDs:

ls ../../garrigos/fastq/*fastq.gz | sed -e 's#.*/##' -e 's/_R[12]\.fastq\.gz$//' | sort -u | head -n 3
S63
S64
S65

6.4 Capture groups and backreferences

A powerful regex feature is that you can capture parts of a match by putting them in parentheses ( ), and then recall them in the replacement with \1 (the first group), \2 (the second group), and so on. (Use sed -E, so that the parentheses don’t need backslashes.)

For example, genomic regions are often written as <seqID>:<start>-<end>. To turn such a region into three Tab-separated columns:

echo "NC_068937.1:2046-110808" | sed -E 's/^([^:]+):([0-9]+)-([0-9]+)$/\1\t\2\t\3/'
NC_068937.1     2046    110808

Here, ([^:]+) captures everything up to the colon, and each ([0-9]+) captures a number.

Capture groups are also a way to extract information from a line. For example, to extract the gene ID from the GTF’s gene lines, match the entire line, and replace it by just the captured ID:

awk -F '\t' '$3 == "gene"' annot.gtf | sed -E 's/.*gene_id "([^"]+)".*/\1/' | head -n 3
LOC120427725
LOC120427726
LOC128092783

Let’s try the same with the transcript IDs, on all lines of the GTF file:

grep -v '^#' annot.gtf | sed -E 's/.*transcript_id "([^"]+)".*/\1/' | head -n 3
NC_068937.1     Gnomon  gene    2046    110808  .       +       .       gene_id "LOC120427725"; transcript_id ""; db_xref "GeneID:120427725"; description "homeotic protein deformed"; gbkey "Gene"; gene "LOC120427725"; gene_biotype "protein_coding";
XM_052707445.1
XM_052707445.1

What went wrong in the first line? Gene lines have an empty transcript ID (transcript_id ""), which doesn’t match [^"]+ (one or more non-quote characters). And when the pattern doesn’t match, sed simply prints the line unchanged.

The solution is to combine the -n option with the p modifier, which prints a line only if a substitution was made:

grep -v '^#' annot.gtf | sed -E -n 's/.*transcript_id "([^"]+)".*/\1/p' | head -n 3
XM_052707445.1
XM_052707445.1
XM_052707445.1

Exercise 5: sed

  1. Above, sed 's#.*/##' removed the dirs from a file path. What do you think the command below will print instead? Run it to check.

    echo "../../garrigos/fastq/S63_R1.fastq.gz" | sed 's#/.*##'

      Click to see the solution

     
    ..

    The pattern /.* matches the first / and everything after it, so only the first .. is left.

  2. Create a version of the metadata file in which the time column holds the number of hours post-infection: 240 instead of 10dpi, and 24 instead of 24hpi. Also rename the column header from time to time_h — make sure that only the header is affected.

      Click for a hint

    Use -e for multiple substitutions, and an address to restrict the header change to line 1.

      Click to see the solution

     
    sed -e '1s/time/time_h/' -e 's/10dpi/240/' -e 's/24hpi/24/' metadata.tsv > metadata_hours.tsv
head -n 3 metadata_hours.tsv
    sample_id       time_h  treatment
S82     240     cathemerium
S75     240     cathemerium

    In this case, s/time/time_h/ without the 1 address would have worked too, because no other line contains time. But with the address, you can be sure that only the header line is changed.

  3. Do the reverse of the sed capture group example above: for each gene in the GTF file, print its region in the format <seqID>:<start>-<end>, e.g. NC_068937.1:2046-110808. Can you think of a way to do this with awk instead?

      Click to see the solution

    With sed, first select the gene lines and the relevant columns, and then capture and rearrange the three columns:

    awk -F '\t' '$3 == "gene"' annot.gtf | cut -f 1,4,5 |
    sed -E 's/^([^\t]+)\t([0-9]+)\t([0-9]+)$/\1:\2-\3/' | head -n 3
    NC_068937.1:2046-110808
NC_068937.1:48457-50554
NC_068937.1:76295-82374

    With awk, this is a lot simpler: print the columns separated by the desired text (between quotes):

    awk -F '\t' '$3 == "gene" { print $1 ":" $4 "-" $5 }' annot.gtf | head -n 3

    As you can see, there is often more than one way to do something with these tools.

7 Putting it together: genes per chromosome

As a final example, let’s combine several of the tools on this page to answer a biological question: how gene-dense is each chromosome? To compute the number of genes per million bases (Mb) for each sequence, you’ll need information from two files: sequence lengths from the genome FASTA, and gene counts from the GTF. You’ll then merge these with a new command, join.

7.1 Step 1: Sequence lengths from the FASTA file

Computing the length of each sequence in a FASTA file takes a slightly more complex awk command, because each sequence is spread across many lines:

awk '/^>/ { if (id) print id "\t" len; id = substr($1, 2); len = 0; next }
     { len += length($0) }
     END { print id "\t" len }' genome.fna | sort -k1,1 > seq_lengths.tsv

Here is what this does:

  • For header lines (pattern /^>/):
    • If this isn’t the first header, print the ID and length of the previous sequence.
    • Store the new sequence ID: the first field without the > (substr($1, 2) takes everything from the 2nd character onwards).
    • Reset the length to 0, and skip to the next line (next).
  • For sequence lines (no pattern, but header lines never get here because of next): add the number of bases on the line to the length.
  • At the end, print the ID and length of the last sequence, which hasn’t been printed yet.
  • Finally, sort by sequence ID, which join will need.
head -n 4 seq_lengths.tsv
NC_015079.1     14856
NC_068937.1     139153999
NC_068938.1     224097512
NC_068939.1     186194774

7.2 Step 2: Gene counts from the GTF file

To count the genes on each sequence, make a frequency table of column 1 for the gene lines. Then, swap the two columns of uniq -c’s output with awk, so that the sequence ID comes first, as in the other file:

awk -F '\t' '$3 == "gene"' annot.gtf | cut -f 1 | sort | uniq -c |
    awk '{ print $2 "\t" $1 }' > gene_counts.tsv

head -n 4 gene_counts.tsv
NC_015079.1     13
NC_068937.1     5173
NC_068938.1     7532
NC_068939.1     6248

7.3 Step 3: Merging the two files with join

join merges two files based on a shared column — by default, the first column of each file. For every sequence ID that occurs in both files, it prints the ID followed by the remaining columns of both files:

join seq_lengths.tsv gene_counts.tsv | head -n 4
NC_015079.1 14856 13
NC_068937.1 139153999 5173
NC_068938.1 224097512 7532
NC_068939.1 186194774 6248

join requires both files to be sorted on the join column. If they aren’t, it will miss matches (usually, but not always, with a warning).

The joined output has fewer lines than seq_lengths.tsv. Which sequences are missing, and why? join only prints lines whose ID occurs in both files (this is called an “inner join”). Sequences without any annotated genes don’t occur in gene_counts.tsv, so they’re left out. To also keep the lines that only occur in the first file, use join -a 1.

7.4 Step 4: Computing gene density

Finally, use awk to divide the number of genes (column 3) by the length in Mb (column 2 divided by 1,000,000), and only keep the sequences with NC IDs:

join seq_lengths.tsv gene_counts.tsv |
    awk '{ print $1, $2, $3, $3 / ($2 / 1000000) }' |
    grep '^NC'
NC_015079.1 14856 13 875.067
NC_068937.1 139153999 5173 37.1746
NC_068938.1 224097512 7532 33.6104
NC_068939.1 186194774 6248 33.5563

So, the three chromosomes have a similar gene density of about 34–37 genes per Mb. The mitochondrial genome (NC_015079.1) is far more gene-dense: its 13 genes are packed into fewer than 15,000 bases.

For more complex merges of tables, R is a better tool: its dplyr package, which you’ll learn about in week 13, has a set of flexible join functions.

8 Wrap-up: which tool when?

Task Tool(s)
Find lines that contain a (regex) pattern grep
Extract the matching part of lines grep -o, or sed -E with capture groups
Select columns cut; awk to also reorder them
Select rows based on the values in a column awk
Calculate with values in columns awk
Replace text sed; tr for single characters
Count or summarize values sort and uniq -c
Merge two tables join

These tools shine for quick, one-off questions about (large) files, and for small processing steps in shell scripts. But when a pipeline gets long or hard to read, it’s often better to switch to:

  • R or Python, for more complex data wrangling, statistics, and plotting (you’ll start learning R in week 12).
  • Dedicated bioinformatics tools, for format-specific tasks like filtering or reverse complementing reads and sequences. Such tools handle the details of file formats (like multi-line FASTA sequences, or soft-masked bases) for you.

Finally, now that you’re done with this page, remove the large genome FASTA file:

rm genome.fna


Back to top

References

Garrigós, Marta, Guillem Ylla, Josué Martínez-de la Puente, Jordi Figuerola, and María José Ruiz-López. 2025. “Two Avian Plasmodium Species Trigger Different Transcriptional Responses on Their Vector Culex pipiens.” Molecular Ecology 34 (15): e17240. https://doi.org/10.1111/mec.17240.

Footnotes

  1. The plus 1 is needed because GTF coordinates are “1-based and inclusive”: a feature with start and end position 2046 is 1 base long. Watch out: in other formats, like BED, this is not the case.↩︎