Using bash in bioinformatics

The Unix commands shown in this notebook are quite intimidating at first. There’s a lot of material becaue I try to cover all the stuff you willl need to construct a typicla bioinformatics pipeline.

Try to understand the expressions and play aound to see what varying them will do. You do not need to memorize the commands now - you can always come back to this notebook for reference. You will also get more practice when you do the Bioiinformatics sessions to construct Unix pipelines for processing your RNA-seq data sets.

Hack to handle broken pipes.

[1]:

cleanup () {
    :
}

trap "cleanup" SIGPIPE

Treat unset variables as an error, and immediately exit.

[2]:

set -u

0.1. Files, directories and working with text

Quickly reveiw what the following navigation and file/directory manipulation commands do:

>, >>, <, |, ls, cd, pwd, cp, mv, rm, mkdir, rmdir

This session dives quite deep into Unix shell scripting and invovles quite a bit of text processing. We review some basic commands for working with text in Unix.

This is the file we will warm up with

We create it using a HEREDOC that allows us to construct multi-line files and redirrect it to the file junk.txt

[3]:

cat > junk.txt <<EOF
1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG
EOF

We can see the contents of a file with cat

[4]:

cat junk.txt

1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

Using head and tail

[5]:

head -3 junk.txt

1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA

[6]:

tail -3 junk.txt

2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

[7]:

tail +3 junk.txt

3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

Pipling to chain commnads

[8]:

head -3 junk.txt | tail -2

2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA

Searching within a file with grep

  • grep takes a REGEX (regular expression) as the query and a file (or files) to search

  • Normally grep shows lines that match the REGEX argument

  • The -v flag shows the lines that DO NOT match the REGEX argument

[9]:

grep "[2|4]" junk.txt

2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

[10]:

grep -v "[2|4]" junk.txt

1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA

Use a regular expression with grep

[11]:

grep "CCC.*CCC" junk.txt

1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG

Get rid of numbers with cut

cut extracts columns separated by delimiters - here the delimiter is a space -d ' ' - and we cut the second column -f 2

[12]:

cut -d ' ' -f 2 junk.txt

CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

Get rid of numbers with grep

The -o flag to grep shows only the strings mathcing the regular expression

[13]:

grep -o '[A-Z]*' junk.txt

CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG

Convert DNA to RNA with tr

[14]:

cat junk.txt | tr T U

1. CCCUAACCCCCUAACCCCCUAACCCUCAGUCGGGGAGGCGACAAUAGCUG
2. GUCAUAUGUUCUGUACGUUAUUGGCCAACUGAUCAUACCUGAAUCGAGCC
3. GAACCGGGAUUAUCAAAGACGAACAUGGUCGGGUCCUUGAACCAAACGAA
4. UCUCCGUCCGCUGGCGUGUUUUUCUUUUCUCAAGUGGGCAAGUUACCCGG

Reverse complement with tr

[15]:

grep -o '[ACTG]*' junk.txt |
tr ACGT TGCA |
rev

CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG
GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC
TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC
CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA

Split command over multiple lines

Use \ to continue command to next line

[16]:

grep -o '[ACTG]*' junk.txt | tr \
ACGT TGCA | rev

CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG
GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC
TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC
CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA

Actually the line continuation character is not necesssary if you break at a pipe |

[17]:

grep -o '[ACTG]*' junk.txt |
tr ACGT TGCA |
rev

CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG
GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC
TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC
CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA

0.2 Getting files

We will download Chromosomes 7 and 8 of Crytococcus neoformans in FASTA format from NCBI Nucleotide. We can find their IDs through a web query.

We find that the IDs are

CM000046.1 for chormosome 7
CM000047.1 for chormosome 8

We will use the NCBI Entrez tool efetch to download the FASTA files to a folder seqs.

[18]:

### setup efetch command
cd ~
#/bin/bash
perl -MNet::FTP -e \
    '$ftp = new Net::FTP("ftp.ncbi.nlm.nih.gov", Passive => 1);
     $ftp->login; $ftp->binary;
     $ftp->get("/entrez/entrezdirect/edirect.tar.gz");'
gunzip -c edirect.tar.gz | tar xf -
rm edirect.tar.gz
#builtin exit
export PATH=${PATH}:$HOME/edirect >& /dev/null || setenv PATH "${PATH}:$HOME/edirect"
./edirect/setup.sh


Trying to establish local installations of any missing Perl modules
(as logged in /home/jovyan/edirect/setup-deps.log).
Please be patient, as this step may take a little while.

Entrez Direct has been successfully downloaded and installed.

In order to complete the configuration process, please execute the following:

  echo "export PATH=\${PATH}:/home/jovyan/edirect" >> $HOME/.bashrc

or manually edit the PATH variable assignment in your .bash_profile file.


[19]:

mkdir -p seqs
efetch -db=nuccore -format=fasta -id=CM000046.1 > seqs/CM000046.1.fa
efetch -db=nuccore -format=fasta -id=CM000047.1 > seqs/CM000047.1.fa

[20]:

ls seqs

CM000046.1.fa      cneo040623.scaffs.fa.gz
CM000047.1.fa      Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
CM000047.1.fa.bak

[21]:

head seqs/CM000046.1.fa

>CM000046.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 7, whole genome shotgun sequence
TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTT
CTTTAATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCC
CTAACCCTTTTATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATA
ACCCCCGAAAAGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAGGA
GGTTGGATTATTTTTTCTTTTCTTTAATAAGTTGTTTATTTAAGTAGTTTCTTTTATTCGGGTAACACAC
ACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCCCATCCAAGGTGCCACGGTTGA
GATTGAGCCCCACACACTGTCGGGACAAAGGAGAAACGACCTTCGGGTCAGAGGTTCCAGCGCTCTGGCC
TTCACTGACTACGACCTGAAGGTTTACTCCCTCGGGGACCGAGACGCGAGGAGCACAGCCACCCCCAGCA
CCCCCAACAGCAAACTGGCCGACTTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTCAGGT

[22]:

head seqs/CM000047.1.fa

>CM000047.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 8, whole genome shotgun sequence
ATGGCCTTCCTGACTACGACCTGAGGTTTACTCCCTCGGCGACCGAGACGCGAGGAGCACAGCCACCCCA
AGCACCCCCAACAGCAAGCTGGCCGAATTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTC
AGGTCGTCTCGAAGAACGCTCCGAAAGTCACTGGTGGGGTCTTTAAACCGATCATCCTTTCCACTGGTGG
CCTGATGAGCAGGAGCACAGCAGACGAATGGAAGGAGTGGAGGGAGGCGATGCCGGTGGGGGGGTTCGAG
AAAATGGAGAAACGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAG
GAGGTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTTCTTTAATTCGGGCAACCC
ACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCTCAATAGCTCTGGCTGGTAC
AGATGGGCAGAAACTCACCTTGCTCCCACGGCATACATCCCTCCCAATAAGTATCTGCAGATGGCGAATA
AACCCCATGCTCTGCGCAAATATATTGAAGGCTGACGACTTTGATGTTATGATCCACCCCATAAACGCTC

0.3 Compressing ard archiving files

In Bioinformatics the most common formats are gzip and bgzip.

[23]:

cd ~/seqs
ls

CM000046.1.fa      cneo040623.scaffs.fa.gz
CM000047.1.fa      Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
CM000047.1.fa.bak

Compressing with gzip

[24]:

gzip *fa
ls

CM000046.1.fa.gz   cneo040623.scaffs.fa.gz
CM000047.1.fa.bak  Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
CM000047.1.fa.gz

Uncompressing with gzip

[25]:

gunzip *gz
ls

CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

Compressing with bgzip

bgzip (Blocked GNU Zip Format) is commonly used in Bioinformatics, and allows random access of the content. BAM, BCF and VCF file formats are typically bgzip compressed.

Unfortunately bgzip does not take multiple files, and we have to use the xargs command to pipe in one file at a time.

[26]:

xargs --help | grep '\-n,'

  -n, --max-args=MAX-ARGS      use at most MAX-ARGS arguments per command line

[27]:

ls *fa | xargs -n1 bgzip
ls

CM000046.1.fa.gz   cneo040623.scaffs.fa.gz
CM000047.1.fa.bak  Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa.gz

Note that bgzip files can be uncompressed with gzip

[28]:

gunzip *gz
ls

CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

Viewing compressed files

Some utilities can work directly with gzipped files. One that we will be using extensively is zcat.

[29]:

ls *fa | xargs -n1 bgzip
ls

CM000046.1.fa.gz   cneo040623.scaffs.fa.gz
CM000047.1.fa.bak  Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa.gz

[30]:

zcat CM000046.1.fa.gz | head

>CM000046.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 7, whole genome shotgun sequence
TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTT
CTTTAATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCC
CTAACCCTTTTATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATA
ACCCCCGAAAAGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAGGA
GGTTGGATTATTTTTTCTTTTCTTTAATAAGTTGTTTATTTAAGTAGTTTCTTTTATTCGGGTAACACAC
ACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCCCATCCAAGGTGCCACGGTTGA
GATTGAGCCCCACACACTGTCGGGACAAAGGAGAAACGACCTTCGGGTCAGAGGTTCCAGCGCTCTGGCC
TTCACTGACTACGACCTGAAGGTTTACTCCCTCGGGGACCGAGACGCGAGGAGCACAGCCACCCCCAGCA
CCCCCAACAGCAAACTGGCCGACTTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTCAGGT

[31]:

gunzip *gz
ls

CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

Archival with tar

Often we want to compress an entire direcotry (possibly with subdirectories). In Unix, this is most commonly done with tar. Flags to note

  • c compress

  • x extract

  • z gzip or gunzip

  • f files

  • v verbose

[32]:

tar --help | head -8

Usage: tar [OPTION...] [FILE]...
GNU 'tar' saves many files together into a single tape or disk archive, and can
restore individual files from the archive.

Examples:
  tar -cf archive.tar foo bar  # Create archive.tar from files foo and bar.
  tar -tvf archive.tar         # List all files in archive.tar verbosely.
  tar -xf archive.tar          # Extract all files from archive.tar.

[33]:

cd
ls seqs

CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

[34]:

tar czvf seqs.tar.gz seqs
ls seqs*

seqs/
seqs/CM000047.1.fa.bak
seqs/CM000046.1.fa
seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf
seqs/CM000047.1.fa
seqs/cneo040623.scaffs.fa
seqs.tar.gz

seqs:
CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

View contents of archive

[35]:

tar tzvf seqs.tar.gz

drwxr-xr-x jovyan/users      0 2019-06-26 09:39 seqs/
-rw-r--r-- jovyan/users 1194398 2019-06-26 09:30 seqs/CM000047.1.fa.bak
-rw-r--r-- jovyan/users 1434513 2019-06-26 09:38 seqs/CM000046.1.fa
-rw-r--r-- jovyan/users 27773163 2019-06-26 09:30 seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf
-rw-r--r-- jovyan/users  1194378 2019-06-26 09:38 seqs/CM000047.1.fa
-rw-r--r-- jovyan/users 18907235 2019-06-26 09:39 seqs/cneo040623.scaffs.fa

Remove direcotry and regenerate from archive

[36]:

rm -rf seqs
ls seqs*

seqs.tar.gz

[37]:

tar xzvf seqs.tar.gz
ls seqs*

seqs/
seqs/CM000047.1.fa.bak
seqs/CM000046.1.fa
seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf
seqs/CM000047.1.fa
seqs/cneo040623.scaffs.fa
seqs.tar.gz

seqs:
CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

[38]:

rm seqs.tar.gz

0.4 Checksums

When working with genomic data, we deal with very large files. There is a small risk that these files will be corrupted over time or during data transfer. To ensure that files are not changed, we use a “checksum” function. This is a function that generates an long, essentially random number called a checksum that represents the contents of the file. When the file contents change, so will the checksum. In theory, there is a very small probability that two different files generate the same checksum, but in practice the probability is too small to worry about.

There are several different algorithms for generating the checksums, and at least 3 Unix commands to do so, but they all work very similarly for our purposes.

The strategy is:

  • Generate and store a checksum together with a data file whose integrity you care about

  • When you use or download the data, re-generate the checksum (using the same algorithm e.g. MD5) and compare with the checksum

[39]:

ls ~/seqs

CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa  cneo040623.scaffs.fa

[40]:

md5sum ~/seqs/*fa > MD5_CHECKSUM

We check that the files are unchanged wth the -c flag

[41]:

md5sum -c MD5_CHECKSUM

/home/jovyan/seqs/CM000046.1.fa: OK
/home/jovyan/seqs/CM000047.1.fa: OK
/home/jovyan/seqs/cneo040623.scaffs.fa: OK

Suppsoe one of the files is corrupted

[42]:

echo 'Oops I did it again' >> seqs/CM000047.1.fa

[43]:

md5sum -c MD5_CHECKSUM  \
| cat # this is jusst so we can Run All

/home/jovyan/seqs/CM000046.1.fa: OK
/home/jovyan/seqs/CM000047.1.fa: FAILED
/home/jovyan/seqs/cneo040623.scaffs.fa: OK
md5sum: WARNING: 1 computed checksum did NOT match

Using diff to find changes

Let’s rename the altered file, re-download the orginal and see what was changed.

[44]:

mv seqs/CM000047.1.fa seqs/CM000047.1.fa.bak

[45]:

efetch -db=nuccore -format=fasta -id=CM000047.1 > seqs/CM000047.1.fa

[46]:

diff -u seqs/CM000047.1.fa seqs/CM000047.1.fa.bak \
| cat # this is jusst so we can Run All

--- seqs/CM000047.1.fa  2019-06-26 09:39:10.615210391 +0000
+++ seqs/CM000047.1.fa.bak      2019-06-26 09:39:09.187204139 +0000
@@ -16820,3 +16820,4 @@
 CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG
 ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT
 AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG
+Oops I did it again

Explanaiton from [official documentation]

Header

--- seqs/CM000047.1.fa  2018-06-21 09:58:18.246185758 -0400
+++ seqs/CM000047.1.fa.bak  2018-06-21 09:55:55.967707764 -0400

Explanation

The unified output format starts with a two-line header, which looks like this:

--- from-file from-file-modification-time
+++ to-file to-file-modification-time
The timestamp looks like ‘2002-02-21 23:30:39.942229878 -0800’ to indicate the date, time with fractional seconds, and time zone. The fractional seconds are omitted on hosts that do not support fractional timestamps.

You can change the header’s content with the --label=label option. See Alternate Names.

Hunk

@@ -16820,3 +16820,4 @@
 CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG
 ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT
 AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG
+Oops I did it again

Explanation

``` @@ from-file-line-numbers to-file-line-numbers @@ line-from-either-file line-from-either-file…

If a hunk contains just one line, only its start line number appears. Otherwise its line numbers look like ‘start,count’. An empty hunk is considered to start at the line that follows the hunk.

If a hunk and its context contain two or more lines, its line numbers look like ‘start,count’. Otherwise only its end line number appears. An empty hunk is considered to end at the line that precedes the hunk.

The lines common to both files begin with a space character. The lines that actually differ between the two files have one of the following indicator characters in the left print column:

‘+’ A line was added here to the first file.

‘-’ A line was removed here from the first file. ```

Verify

We can confirm that a line was added at the position given by diff

[47]:

cat seqs/CM000047.1.fa | tail +16820 | head -3

CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG
ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT
AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG

[48]:

cat seqs/CM000047.1.fa.bak | tail +16820 | head -4

CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG
ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT
AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG
Oops I did it again

1. Working with FASTA files

From NCBI FASTA specification

A sequence in FASTA format begins with a single-line description, followed by lines of sequence data. The description line (defline) is distinguished from the sequence data by a greater-than (">") symbol at the beginning. It is recommended that all lines of text be shorter than 80 characters in length. An example sequence in FASTA format is:

        >P01013 GENE X PROTEIN (OVALBUMIN-RELATED)
        QIKDLLVSSSTDLDTTLVLVNAIYFKGMWKTAFNAEDTREMPFHVTKQESKPVQMMCMNNSFNVATLPAE
        KMKILELPFASGDLSMLVLLPDEVSDLERIEKTINFEKLTEWTNPNTMEKRRVKVYLPQMKIEEKYNLTS
        VLMALGMTDLFIPSANLTGISSAESLKISQAVHGAFMELSEDGIEMAGSTGVIEDIKHSPESEQFRADHP
        FLFLIKHNPTNTIVYFGRYWSP

Note that a FASTA file may contain mulitple sequences, each with its own description line.

[49]:

cd seqs
wget ftp://ftp.ensemblgenomes.org/pub/release-43/fungi/gtf/cryptococcus_neoformans/Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
wget http://www-sequence.stanford.edu/group/C.neoformans/files/cneo040623.scaffs.fa.gz
cd -

--2019-06-26 09:39:11--  ftp://ftp.ensemblgenomes.org/pub/release-43/fungi/gtf/cryptococcus_neoformans/Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
           => ‘Cryptococcus_neoformans.ASM9104v1.43.gtf.gz’
Resolving ftp.ensemblgenomes.org (ftp.ensemblgenomes.org)... 193.62.197.94
Connecting to ftp.ensemblgenomes.org (ftp.ensemblgenomes.org)|193.62.197.94|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /pub/release-43/fungi/gtf/cryptococcus_neoformans ... done.
==> SIZE Cryptococcus_neoformans.ASM9104v1.43.gtf.gz ... 1546239
==> PASV ... done.    ==> RETR Cryptococcus_neoformans.ASM9104v1.43.gtf.gz ... done.
Length: 1546239 (1.5M) (unauthoritative)

Cryptococcus_neofor 100%[===================>]   1.47M  2.04MB/s    in 0.7s

2019-06-26 09:39:13 (2.04 MB/s) - ‘Cryptococcus_neoformans.ASM9104v1.43.gtf.gz’ saved [1546239]

--2019-06-26 09:39:13--  http://www-sequence.stanford.edu/group/C.neoformans/files/cneo040623.scaffs.fa.gz
Resolving www-sequence.stanford.edu (www-sequence.stanford.edu)... 171.66.102.5
Connecting to www-sequence.stanford.edu (www-sequence.stanford.edu)|171.66.102.5|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 5856251 (5.6M) [text/plain]
Saving to: ‘cneo040623.scaffs.fa.gz’

cneo040623.scaffs.f 100%[===================>]   5.58M  1.65MB/s    in 3.4s

2019-06-26 09:39:16 (1.65 MB/s) - ‘cneo040623.scaffs.fa.gz’ saved [5856251/5856251]

/home/jovyan

Set up variables for files we will be using

[50]:

#DATA_DIR='/data/hts/2018/Granek_4837_180427A5'
DATA_DIR='./seqs'
FASTA='cneo040623.scaffs.fa.gz'
FASTQ='10_MA_C_S3_L003_R1_001.fastq.gz'
GTF='Cryptococcus_neoformans.ASM9104v1.43.gtf.gz'

This is how you get the value of a shell variable

[51]:

echo $DATA_DIR

./seqs

List files in data direcotry

[52]:

ls $DATA_DIR | head -3

CM000046.1.fa
CM000047.1.fa
CM000047.1.fa.bak

Inspect the contents of a zipped file.

zcat is like cat but works on gzipped files.

[53]:

zcat $DATA_DIR/$FASTA | head

>chr01.b3501.040506 5 contigs (2263073 bp)
CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG
GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC
GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA
TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG
GGATTTATACAACGATAAGAGGCTATCCGCTGACTATAAATTGTGTTAGC
TGATCCAGATCCGCAGTTGGCACAAACCGATGTTTCCTTTTCCTCAGGTC
CTGAATATTTCTCTTCCATCAAGGCTCCCAACCCCGAGGGGAGTATTTCC
ACCCGATCTGATTCCAAGAGGTCCAGCGTCAACCAGGTACGTGAATGTTT
TCTTCTTTAGTGGTGGGGACCCTATGCTTGTCCCACAGAACTGCACTGAA

Show headers only

[54]:

zcat $DATA_DIR/$FASTA   | grep '^>' | head

>chr01.b3501.040506 5 contigs (2263073 bp)
>chr02.b3501.040506 6 contigs (1607892 bp)
>chr03.b3501.040506 6 contigs (2025976 bp)
>chr04.b3501.040623 5 contigs (1762351 bp)
>chr05.b3501.040506 7 contigs (1432964 bp)
>chr06.b3501.040616 4 contigs (1411173 bp)
>chr07.b3501.040506 6 contigs (1356306 bp)
>chr08.b3501.040506 3 contigs (1153453 bp)
>chr09.b3501.040506 5 contigs (1072068 bp)
>chr10.b3501.040506 3 contigs (1017008 bp)

Show headers with line numbers

[55]:

zcat $DATA_DIR/$FASTA| grep -n '^>' | head

1:>chr01.b3501.040506 5 contigs (2263073 bp)
45264:>chr02.b3501.040506 6 contigs (1607892 bp)
77423:>chr03.b3501.040506 6 contigs (2025976 bp)
117944:>chr04.b3501.040623 5 contigs (1762351 bp)
153193:>chr05.b3501.040506 7 contigs (1432964 bp)
181854:>chr06.b3501.040616 4 contigs (1411173 bp)
210079:>chr07.b3501.040506 6 contigs (1356306 bp)
237207:>chr08.b3501.040506 3 contigs (1153453 bp)
260278:>chr09.b3501.040506 5 contigs (1072068 bp)
281721:>chr10.b3501.040506 3 contigs (1017008 bp)

Count number of headers (one for each crhomosome)

[56]:

zcat $DATA_DIR/$FASTA  | grep -c '^>'

14

Extract first part of the sequence from chromosome 7

[57]:

zcat $DATA_DIR/$FASTA | sed -n '210079,210082p'

>chr07.b3501.040506 6 contigs (1356306 bp)
TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAA
GTTGTTTATTTAAGTAGTTTCTTTAATTCGGGTAACACACACGACAACCC
AATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCCCTAACCCTTT

Extract first line of chromosome 7 and save to file

[58]:

zcat $DATA_DIR/$FASTA | sed -n '210080p' > dna.txt

Inspect the file we just saved

[59]:

cat dna.txt

TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAA

Find complement

[60]:

cat dna.txt | tr ACGT TGCA

ACCAGAATACTCCTTCTCCTCAAACCTAATAAAAAAGAAAAGAAATTTTT

Find reverse complement

[61]:

cat dna.txt | tr ACGT TGCA | rev

TTTTTAAAGAAAAGAAAAAATAATCCAAACTCCTCTTCCTCATAAGACCA

Delete A and T

[62]:

cat dna.txt | tr -d 'AT'

GGCGGGGGGGGGCC

Only count non A and T

[63]:

cat dna.txt | tr -d 'AT' | tr -d '\n' | wc -c

14

[64]:

man tr | grep '\-d,' -A2

bash: man: command not found

Count C and G

[65]:

grep '[CG]' -o dna.txt | wc -l

14

[66]:

man grep | grep '\-o,' -A2

bash: man: command not found

Specify data directory as a shell variable

2. Working with FASTQ files

From Wikipedia

A FASTQ file normally uses four lines per sequence.

Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description (like a FASTA title line).
Line 2 is the raw sequence letters.
Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.
Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.

Specify target file as a shell variable

[67]:

DATA_DIR='/data/hts2018_pilot/Granek_4837_180427A5'

[68]:

ls $DATA_DIR

10_MA_C_S3_L001_R1_001.fastq.gz   26_RZ_C_S16_L004_R1_001.fastq.gz
10_MA_C_S3_L002_R1_001.fastq.gz   27_MA_P_S38_L001_R1_001.fastq.gz
10_MA_C_S3_L003_R1_001.fastq.gz   27_MA_P_S38_L002_R1_001.fastq.gz
10_MA_C_S3_L004_R1_001.fastq.gz   27_MA_P_S38_L003_R1_001.fastq.gz
10_RZ_C_S11_L001_R1_001.fastq.gz  27_MA_P_S38_L004_R1_001.fastq.gz
10_RZ_C_S11_L002_R1_001.fastq.gz  27_RZ_P_S46_L001_R1_001.fastq.gz
10_RZ_C_S11_L003_R1_001.fastq.gz  27_RZ_P_S46_L002_R1_001.fastq.gz
10_RZ_C_S11_L004_R1_001.fastq.gz  27_RZ_P_S46_L003_R1_001.fastq.gz
11_MA_J_S20_L001_R1_001.fastq.gz  27_RZ_P_S46_L004_R1_001.fastq.gz
11_MA_J_S20_L002_R1_001.fastq.gz  2_MA_C_S1_L001_R1_001.fastq.gz
11_MA_J_S20_L003_R1_001.fastq.gz  2_MA_C_S1_L002_R1_001.fastq.gz
11_MA_J_S20_L004_R1_001.fastq.gz  2_MA_C_S1_L003_R1_001.fastq.gz
11_RZ_J_S28_L001_R1_001.fastq.gz  2_MA_C_S1_L004_R1_001.fastq.gz
11_RZ_J_S28_L002_R1_001.fastq.gz  2_RZ_C_S9_L001_R1_001.fastq.gz
11_RZ_J_S28_L003_R1_001.fastq.gz  2_RZ_C_S9_L002_R1_001.fastq.gz
11_RZ_J_S28_L004_R1_001.fastq.gz  2_RZ_C_S9_L003_R1_001.fastq.gz
12_MA_P_S36_L001_R1_001.fastq.gz  2_RZ_C_S9_L004_R1_001.fastq.gz
12_MA_P_S36_L002_R1_001.fastq.gz  2_TOT_C_S17_L001_R1_001.fastq.gz
12_MA_P_S36_L003_R1_001.fastq.gz  2_TOT_C_S17_L002_R1_001.fastq.gz
12_MA_P_S36_L004_R1_001.fastq.gz  2_TOT_C_S17_L003_R1_001.fastq.gz
12_RZ_P_S44_L001_R1_001.fastq.gz  2_TOT_C_S17_L004_R1_001.fastq.gz
12_RZ_P_S44_L002_R1_001.fastq.gz  35_MA_P_S39_L001_R1_001.fastq.gz
12_RZ_P_S44_L003_R1_001.fastq.gz  35_MA_P_S39_L002_R1_001.fastq.gz
12_RZ_P_S44_L004_R1_001.fastq.gz  35_MA_P_S39_L003_R1_001.fastq.gz
13_MA_J_S21_L001_R1_001.fastq.gz  35_MA_P_S39_L004_R1_001.fastq.gz
13_MA_J_S21_L002_R1_001.fastq.gz  35_RZ_P_S47_L001_R1_001.fastq.gz
13_MA_J_S21_L003_R1_001.fastq.gz  35_RZ_P_S47_L002_R1_001.fastq.gz
13_MA_J_S21_L004_R1_001.fastq.gz  35_RZ_P_S47_L003_R1_001.fastq.gz
13_RZ_J_S29_L001_R1_001.fastq.gz  35_RZ_P_S47_L004_R1_001.fastq.gz
13_RZ_J_S29_L002_R1_001.fastq.gz  36_MA_J_S24_L001_R1_001.fastq.gz
13_RZ_J_S29_L003_R1_001.fastq.gz  36_MA_J_S24_L002_R1_001.fastq.gz
13_RZ_J_S29_L004_R1_001.fastq.gz  36_MA_J_S24_L003_R1_001.fastq.gz
14_MA_C_S4_L001_R1_001.fastq.gz   36_MA_J_S24_L004_R1_001.fastq.gz
14_MA_C_S4_L002_R1_001.fastq.gz   36_RZ_J_S32_L001_R1_001.fastq.gz
14_MA_C_S4_L003_R1_001.fastq.gz   36_RZ_J_S32_L002_R1_001.fastq.gz
14_MA_C_S4_L004_R1_001.fastq.gz   36_RZ_J_S32_L003_R1_001.fastq.gz
14_RZ_C_S12_L001_R1_001.fastq.gz  36_RZ_J_S32_L004_R1_001.fastq.gz
14_RZ_C_S12_L002_R1_001.fastq.gz  38_MA_P_S40_L001_R1_001.fastq.gz
14_RZ_C_S12_L003_R1_001.fastq.gz  38_MA_P_S40_L002_R1_001.fastq.gz
14_RZ_C_S12_L004_R1_001.fastq.gz  38_MA_P_S40_L003_R1_001.fastq.gz
15_MA_C_S5_L001_R1_001.fastq.gz   38_MA_P_S40_L004_R1_001.fastq.gz
15_MA_C_S5_L002_R1_001.fastq.gz   38_RZ_P_S48_L001_R1_001.fastq.gz
15_MA_C_S5_L003_R1_001.fastq.gz   38_RZ_P_S48_L002_R1_001.fastq.gz
15_MA_C_S5_L004_R1_001.fastq.gz   38_RZ_P_S48_L003_R1_001.fastq.gz
15_RZ_C_S13_L001_R1_001.fastq.gz  38_RZ_P_S48_L004_R1_001.fastq.gz
15_RZ_C_S13_L002_R1_001.fastq.gz  3_MA_J_S19_L001_R1_001.fastq.gz
15_RZ_C_S13_L003_R1_001.fastq.gz  3_MA_J_S19_L002_R1_001.fastq.gz
15_RZ_C_S13_L004_R1_001.fastq.gz  3_MA_J_S19_L003_R1_001.fastq.gz
16_MA_P_S37_L001_R1_001.fastq.gz  3_MA_J_S19_L004_R1_001.fastq.gz
16_MA_P_S37_L002_R1_001.fastq.gz  3_RZ_J_S27_L001_R1_001.fastq.gz
16_MA_P_S37_L003_R1_001.fastq.gz  3_RZ_J_S27_L002_R1_001.fastq.gz
16_MA_P_S37_L004_R1_001.fastq.gz  3_RZ_J_S27_L003_R1_001.fastq.gz
16_RZ_P_S45_L001_R1_001.fastq.gz  3_RZ_J_S27_L004_R1_001.fastq.gz
16_RZ_P_S45_L002_R1_001.fastq.gz  3_TOT_J_S34_L001_R1_001.fastq.gz
16_RZ_P_S45_L003_R1_001.fastq.gz  3_TOT_J_S34_L002_R1_001.fastq.gz
16_RZ_P_S45_L004_R1_001.fastq.gz  3_TOT_J_S34_L003_R1_001.fastq.gz
1_MA_J_S18_L001_R1_001.fastq.gz   3_TOT_J_S34_L004_R1_001.fastq.gz
1_MA_J_S18_L002_R1_001.fastq.gz   40_MA_J_S25_L001_R1_001.fastq.gz
1_MA_J_S18_L003_R1_001.fastq.gz   40_MA_J_S25_L002_R1_001.fastq.gz
1_MA_J_S18_L004_R1_001.fastq.gz   40_MA_J_S25_L003_R1_001.fastq.gz
1_RZ_J_S26_L001_R1_001.fastq.gz   40_MA_J_S25_L004_R1_001.fastq.gz
1_RZ_J_S26_L002_R1_001.fastq.gz   40_RZ_J_S33_L001_R1_001.fastq.gz
1_RZ_J_S26_L003_R1_001.fastq.gz   40_RZ_J_S33_L002_R1_001.fastq.gz
1_RZ_J_S26_L004_R1_001.fastq.gz   40_RZ_J_S33_L003_R1_001.fastq.gz
21_MA_C_S6_L001_R1_001.fastq.gz   40_RZ_J_S33_L004_R1_001.fastq.gz
21_MA_C_S6_L002_R1_001.fastq.gz   45_MA_P_S41_L001_R1_001.fastq.gz
21_MA_C_S6_L003_R1_001.fastq.gz   45_MA_P_S41_L002_R1_001.fastq.gz
21_MA_C_S6_L004_R1_001.fastq.gz   45_MA_P_S41_L003_R1_001.fastq.gz
21_RZ_C_S14_L001_R1_001.fastq.gz  45_MA_P_S41_L004_R1_001.fastq.gz
21_RZ_C_S14_L002_R1_001.fastq.gz  45_RZ_P_S49_L001_R1_001.fastq.gz
21_RZ_C_S14_L003_R1_001.fastq.gz  45_RZ_P_S49_L002_R1_001.fastq.gz
21_RZ_C_S14_L004_R1_001.fastq.gz  45_RZ_P_S49_L003_R1_001.fastq.gz
22_MA_C_S7_L001_R1_001.fastq.gz   45_RZ_P_S49_L004_R1_001.fastq.gz
22_MA_C_S7_L002_R1_001.fastq.gz   47_MA_P_S42_L001_R1_001.fastq.gz
22_MA_C_S7_L003_R1_001.fastq.gz   47_MA_P_S42_L002_R1_001.fastq.gz
22_MA_C_S7_L004_R1_001.fastq.gz   47_MA_P_S42_L003_R1_001.fastq.gz
22_RZ_C_S15_L001_R1_001.fastq.gz  47_MA_P_S42_L004_R1_001.fastq.gz
22_RZ_C_S15_L002_R1_001.fastq.gz  47_RZ_P_S50_L001_R1_001.fastq.gz
22_RZ_C_S15_L003_R1_001.fastq.gz  47_RZ_P_S50_L002_R1_001.fastq.gz
22_RZ_C_S15_L004_R1_001.fastq.gz  47_RZ_P_S50_L003_R1_001.fastq.gz
23_MA_J_S22_L001_R1_001.fastq.gz  47_RZ_P_S50_L004_R1_001.fastq.gz
23_MA_J_S22_L002_R1_001.fastq.gz  4_MA_P_S35_L001_R1_001.fastq.gz
23_MA_J_S22_L003_R1_001.fastq.gz  4_MA_P_S35_L002_R1_001.fastq.gz
23_MA_J_S22_L004_R1_001.fastq.gz  4_MA_P_S35_L003_R1_001.fastq.gz
23_RZ_J_S30_L001_R1_001.fastq.gz  4_MA_P_S35_L004_R1_001.fastq.gz
23_RZ_J_S30_L002_R1_001.fastq.gz  4_RZ_P_S43_L001_R1_001.fastq.gz
23_RZ_J_S30_L003_R1_001.fastq.gz  4_RZ_P_S43_L002_R1_001.fastq.gz
23_RZ_J_S30_L004_R1_001.fastq.gz  4_RZ_P_S43_L003_R1_001.fastq.gz
24_MA_J_S23_L001_R1_001.fastq.gz  4_RZ_P_S43_L004_R1_001.fastq.gz
24_MA_J_S23_L002_R1_001.fastq.gz  4_TOT_P_S51_L001_R1_001.fastq.gz
24_MA_J_S23_L003_R1_001.fastq.gz  4_TOT_P_S51_L002_R1_001.fastq.gz
24_MA_J_S23_L004_R1_001.fastq.gz  4_TOT_P_S51_L003_R1_001.fastq.gz
24_RZ_J_S31_L001_R1_001.fastq.gz  4_TOT_P_S51_L004_R1_001.fastq.gz
24_RZ_J_S31_L002_R1_001.fastq.gz  9_MA_C_S2_L001_R1_001.fastq.gz
24_RZ_J_S31_L003_R1_001.fastq.gz  9_MA_C_S2_L002_R1_001.fastq.gz
24_RZ_J_S31_L004_R1_001.fastq.gz  9_MA_C_S2_L003_R1_001.fastq.gz
26_MA_C_S8_L001_R1_001.fastq.gz   9_MA_C_S2_L004_R1_001.fastq.gz
26_MA_C_S8_L002_R1_001.fastq.gz   9_RZ_C_S10_L001_R1_001.fastq.gz
26_MA_C_S8_L003_R1_001.fastq.gz   9_RZ_C_S10_L002_R1_001.fastq.gz
26_MA_C_S8_L004_R1_001.fastq.gz   9_RZ_C_S10_L003_R1_001.fastq.gz
26_RZ_C_S16_L001_R1_001.fastq.gz  9_RZ_C_S10_L004_R1_001.fastq.gz
26_RZ_C_S16_L002_R1_001.fastq.gz  Granek_4837_180427A5.checksum
26_RZ_C_S16_L003_R1_001.fastq.gz  README.rtf

Show the contents of a zipped file

[69]:

zcat $DATA_DIR/$FASTQ | head -n 10

@NB501800:128:HTYKFBGX5:3:11401:19195:1019 1:N:0:ATTACTCG+AGGATAGG
NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGC
+
#AAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEE/<<
@NB501800:128:HTYKFBGX5:3:11401:24339:1020 1:N:0:ATTACTCG+AGGATAGG
NCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTC
+
#AAAAEAAEEEEEEEEEEEEEEEE/EEE/EEEEEEEEEEEEEEEE<A/EEEE6EEE//EEEAEEEEEEE/EAEEEA
@NB501800:128:HTYKFBGX5:3:11401:24140:1020 1:N:0:ATTACTCG+AGGATAGG
NTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAG

Extracting just the sequence

We need to extract every 4th line starting from the 2nd. We cna use a tool called sed (stream editor) to do so.

[70]:

zcat $DATA_DIR/$FASTQ | sed -n '2~4p;' | head

NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGC
NCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTC
NTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAG
NGGGCCTTGGTCCTGCCCTCAGCCTTGATAGGACCCTTGCCGTCGATGGGGTTACCCAAAGCGTCAACAACTCGGC
NACCTTTCGGAGACCTCGGTGGGTCTTCCTGGGGAGACGGGTAGTACCCCATCGGGCGGTAACACCCTCGTAACCG
NCCGAACTGTTGAAAAGAGCTTCCACGCTCCAAAGGCAGGGACAAGAGTAAAGAGCCACCACAGCCAATCTCCAAA
NGCAAACTGAACGTCCTCAGTCTCCATTCGCTGGAAGTCCGATTCAGGACGGGCAGCGTCATCGATAGAGAATGGA
NGCATTACCAGATGCACTGATAGCAATGATCGTCTCATCTTCCAAGTACTTGCTCACAGCGATACTATCACTCAAA
GGCGATCTTCCTCTTGGGGAAGTGAACGTTCTTCTTGTGCTGCTTCTCGGTGATAGCAGAAGCCTCCTTGGTAGTC
CGGCGATGATGAGAGTGAAGGGATCTTGCCACATGAGAGGGGATGATCGACGGGGACGAGGACGGATGGGGATGAT

Explanation

sed -n '2~4p'

  • -n means to suppress all output

  • 2~4p means print (p) every 4th line starting at line 2

Joining lines

Getting sequence as a single string. We just need to delete newlline (\n) characters.

Note: We can break long commands uisng the line continuation character

[71]:

zcat $DATA_DIR/$FASTQ | \
sed -n '2~4p;' | \
head | \
tr -d '\n'

NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGCNCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTCNTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAGNGGGCCTTGGTCCTGCCCTCAGCCTTGATAGGACCCTTGCCGTCGATGGGGTTACCCAAAGCGTCAACAACTCGGCNACCTTTCGGAGACCTCGGTGGGTCTTCCTGGGGAGACGGGTAGTACCCCATCGGGCGGTAACACCCTCGTAACCGNCCGAACTGTTGAAAAGAGCTTCCACGCTCCAAAGGCAGGGACAAGAGTAAAGAGCCACCACAGCCAATCTCCAAANGCAAACTGAACGTCCTCAGTCTCCATTCGCTGGAAGTCCGATTCAGGACGGGCAGCGTCATCGATAGAGAATGGANGCATTACCAGATGCACTGATAGCAATGATCGTCTCATCTTCCAAGTACTTGCTCACAGCGATACTATCACTCAAAGGCGATCTTCCTCTTGGGGAAGTGAACGTTCTTCTTGTGCTGCTTCTCGGTGATAGCAGAAGCCTCCTTGGTAGTCCGGCGATGATGAGAGTGAAGGGATCTTGCCACATGAGAGGGGATGATCGACGGGGACGAGGACGGATGGGGATGAT

3.Working wtih GTF files

From UCSD Genome Browser

The GTF (General Transfer Format) is an extension of GFF version 2
and used to represent transcription models. GFF (General Feature Format)
consists of one line per feature, each containing 9 columns of data.

Fields

Fields are tab-separated. Also, all but the final field in each
feature line must contain a value; "empty" columns are denoted
with a '.'

    seqname   - name of the chromosome or scaffold; chromosome names
                without a 'chr'
    source    - name of the program that generated this feature, or
                the data source (database or project name)
    feature   - feature type name. Current allowed features are
                {gene, transcript, exon, CDS, Selenocysteine, start_codon,
                stop_codon and UTR}
    start     - start position of the feature, with sequence numbering
                starting at 1.
    end       - end position of the feature, with sequence numbering
                starting at 1.
    score     - a floating point value indiciating the score of a feature
    strand    - defined as + (forward) or - (reverse).
    frame     - one of '0', '1' or '2'. Frame indicates the number of base pairs
                before you encounter a full codon. '0' indicates the feature
                begins with a whole codon. '1' indicates there is an extra
                base (the 3rd base of the prior codon) at the start of this feature.
                '2' indicates there are two extra bases (2nd and 3rd base of the
                prior exon) before the first codon. All values are given with
                relation to the 5' end.
    attribute - a semicolon-separated list of tag-value pairs (separated by a space),
                providing additional information about each feature. A key can be
                repeated multiple times.

Attributes

The following attributes are available. All attributes are semi-colon
separated pairs of keys and values.

- gene_id: The stable identifier for the gene
- gene_version: The stable identifier version for the gene
- gene_name: The official symbol of this gene
- gene_source: The annotation source for this gene
- gene_biotype: The biotype of this gene
- transcript_id: The stable identifier for this transcript
- transcript_version: The stable identifier version for this transcript
- transcript_name: The symbold for this transcript derived from the gene name
- transcript_source: The annotation source for this transcript
- transcript_biotype: The biotype for this transcript
- exon_id: The stable identifier for this exon
- exon_version: The stable identifier version for this exon
- exon_number: Position of this exon in the transcript
- ccds_id: CCDS identifier linked to this transcript
- protein_id: Stable identifier for this transcript's protein
- protein_version: Stable identifier version for this transcript's protein
- tag: A collection of additional key value tags
- transcript_support_level: Ranking to assess how well a transcript is supported (from 1 to 5)

[72]:

DATA_DIR=~/seqs

[73]:

ls ${DATA_DIR}

CM000046.1.fa         cneo040623.scaffs.fa.gz
CM000047.1.fa         Cryptococcus_neoformans.ASM9104v1.43.gtf
CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf.gz
cneo040623.scaffs.fa

Inspect the GTF file

[74]:

zcat $DATA_DIR/$GTF | head

#!genome-build ASM9104v1
#!genome-version ASM9104v1
#!genome-date 2016-06
#!genome-build-accession GCA_000091045.1
#!genebuild-last-updated 2016-06
1       ena     gene    4839    6236    .       +       .       gene_id "CNA00010"; gene_source "ena"; gene_biotype "protein_coding";
1       ena     transcript      4839    6236    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";
1       ena     exon    4839    5178    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-1"; exon_version "1";
1       ena     CDS     4914    5178    .       +       0       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AAW41410";
1       ena     start_codon     4914    4916    .       +       0       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";

Ignore comment lines version 1

This uses tail to take all lines starting from line 6

[75]:

zcat $DATA_DIR/$GTF | tail +6 | head -3

1       ena     gene    4839    6236    .       +       .       gene_id "CNA00010"; gene_source "ena"; gene_biotype "protein_coding";
1       ena     transcript      4839    6236    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";
1       ena     exon    4839    5178    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-1"; exon_version "1";

Ingore comment lines 2

This uses grep to ignore (-v) lines starting with ^#

[76]:

zcat $DATA_DIR/$GTF | grep -v '^#' | head -3

1       ena     gene    4839    6236    .       +       .       gene_id "CNA00010"; gene_source "ena"; gene_biotype "protein_coding";
1       ena     transcript      4839    6236    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";
1       ena     exon    4839    5178    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-1"; exon_version "1";

Find all unique feature types in GTF

[77]:

zcat $DATA_DIR/$GTF | grep -v '^#' | cut -f 3 | sort | uniq

CDS
exon
five_prime_utr
gene
start_codon
stop_codon
three_prime_utr
transcript

Get counts of all unique feature types in GTF

[78]:

zcat $DATA_DIR/$GTF | grep -v '^#' | cut -f 3 | sort | uniq -c

  44079 CDS
  44271 exon
   6738 five_prime_utr
   6813 gene
   6894 start_codon
   2369 stop_codon
   6718 three_prime_utr
   7044 transcript

Find the the attributes of the first 3 exons

[79]:

zcat $DATA_DIR/$GTF | grep -v '^#' | awk '$3 == "exon"' | head -3

1       ena     exon    4839    5178    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-1"; exon_version "1";
1       ena     exon    5241    5377    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "2"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-2"; exon_version "1";
1       ena     exon    5453    5560    .       +       .       gene_id "CNA00010"; transcript_id "AAW41410"; exon_number "3"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AAW41410-3"; exon_version "1";

Store attributes of all genes in a file

It is common to use Unix commands to do simple processing and save the results to a file, then use a scripting language such as R or Python to furhter process the file.

[80]:

zcat $DATA_DIR/$GTF | grep -v '^#' | awk '$3 == "gene"' | cut -f 9 > genes.txt

[81]:

head -6 genes.txt

gene_id "CNA00010"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00020"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00030"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00040"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00050"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00060"; gene_source "ena"; gene_biotype "protein_coding";

4. Looping

Loop throgh first 10 lines

  • Creae a counter variable LINE_NUM and initialize it to 0

  • Redirect genes.txt to a while loop

  • Whlie loop reads each line into the variable LINE

  • The [[ CONDITION ]] tests if CONDITION is true

  • If LINE_NUM is less than -lt 10 print LINE and increment LINE_NUM

  • Note that (( EXPR )) is used to evaluate numeric expressions

[82]:

LINE_NUM=0
while read LINE
do
    if [[ $LINE_NUM -lt 10 ]]; then
        echo $LINE
        (( LINE_NUM++ ))
    fi
done < genes.txt

gene_id "CNA00010"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00020"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00030"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00040"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00050"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00060"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00070"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00075"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00080"; gene_source "ena"; gene_biotype "protein_coding";
gene_id "CNA00090"; gene_source "ena"; gene_biotype "protein_coding";

Loop through first 10 tRNAs

Almost the same as previous program but we add a new condition to check tha the gene_biotype is tRNA using a REGULAr EXPRESSION via the match operator =~

.*gene_biotype[[:space:]]\"tRNA\".* means

  • .*: Match any number of any character

  • gene_biotype: Match the string gene_biotype

  • [[:space:]]: Match a single space cahracter

  • \"tRNA\": Math the string “tRNA” - we escape the double quotes to match linteral double quote characters

  • .*: Match any number of any character

[83]:

LINE_NUM=0
while read LINE
do
    if [[ $LINE_NUM -lt 10 ]] && [[ $LINE =~ .*gene_biotype[[:space:]]\"tRNA\".* ]] ; then
        echo $LINE
        ((LINE_NUM++))
    fi
done < genes.txt

gene_id "CNA02110"; gene_name "CNA02110"; gene_source "ena"; gene_biotype "tRNA";
gene_id "EBG00005237496"; gene_name "tRNA"; gene_source "Rfam"; gene_biotype "tRNA";
gene_id "CNA02190"; gene_name "CNA02190"; gene_source "ena"; gene_biotype "tRNA";
gene_id "CNA02200"; gene_name "CNA02200"; gene_source "ena"; gene_biotype "tRNA";
gene_id "CNA02280"; gene_name "CNA02280"; gene_source "ena"; gene_biotype "tRNA";
gene_id "EBG00005237502"; gene_name "tRNA"; gene_source "Rfam"; gene_biotype "tRNA";
gene_id "CNA02320"; gene_name "CNA02320"; gene_source "ena"; gene_biotype "tRNA";
gene_id "CNA03230"; gene_name "CNA03230"; gene_source "ena"; gene_biotype "tRNA";
gene_id "EBG00005237494"; gene_name "tRNA"; gene_source "Rfam"; gene_biotype "tRNA";
gene_id "CNA03260"; gene_name "CNA03260"; gene_source "ena"; gene_biotype "tRNA";

Loop through files in a directory and print the filename followed by the first line

  • Use a COMMAND EXPRESSION $(EXPR) to list files matching the glob pattern

  • We loop through each file using a for statement

  • For each file, extract the first line with head -1 and save in variable FIRST

  • The 2>/dev/null ignores output from the stderr stream

  • We use another COMMAND EXPRESSION to print the filename withou the directory (using basename)

  • We then print the first line storred in FIRST

[84]:

DATA_DIR='/data/hts2018_pilot/Granek_4837_180427A5'

[85]:

for FILE in $(ls $DATA_DIR/*_MA_J_S20_*gz)
do
    FIRST=$(zcat $FILE 2>/dev/null | head -1)
    echo $(basename $FILE)
    echo $FIRST
    echo
done

11_MA_J_S20_L001_R1_001.fastq.gz
@NB501800:128:HTYKFBGX5:1:11101:26760:1056 1:N:0:CGCTCATT+NGGATAGG

11_MA_J_S20_L002_R1_001.fastq.gz
@NB501800:128:HTYKFBGX5:2:11101:20327:1041 1:N:0:CGCTCATT+AGGATAGN

11_MA_J_S20_L003_R1_001.fastq.gz
@NB501800:128:HTYKFBGX5:3:11401:10896:1024 1:N:0:CGCTCATT+AGGATAGG

11_MA_J_S20_L004_R1_001.fastq.gz
@NB501800:128:HTYKFBGX5:4:11401:8259:1032 1:N:0:CGCTCATT+AGGATAGG