# 04 Linux command line

question Questions
objectives Objectives

time Time estimation: 15 minutes

Print the word hello to the screen

echo hello


Print the sentence “hello my friends to the screen (with open quotation and without end quotation) Remember you can use the up arrow to go back to previously typed commands

echo "hello my friends


Now the terminal hangs. We typed an incorrect command and the terminal does not know what to do.

What to type when the terminal hangs ?

Ctrl-C If Ctrl-C fails, try hitting ESC. In most of the cases, this will do the trick.

Open the manual of the echo command ?

man echo


The synopsis of this command is:

echo [-n] [string ...]


Things in between square brackets are optional, so it means that you can use echo without options and arguments.

When the manual page is longer than the terminal, you can scroll down the page one line at a time by pressing the down arrow key, or one page at a time by pressing the spacebar. To exit the man page, press q (for quit).

The manual explains that echo prints its argument by default to the screen and then puts the prompt on a new line. The way it does this is by appending a character called a newline (a special character that literally puts the text on a new line). Because echo is often used in programs to print out a sequence of strings not separated by newlines, there is an option to prevent the newline from being inserted.

By reading the man page, find the command to print hello without a newline, and verify that it works as expected. Again remember to use the up arrow to go to previously typed commands.

echo -n hello


Open the manual of the sleep command, find how to make the terminal “sleep” for 5 seconds, and execute the command.

man sleep


According to the manual sleep has a required argument called number representing the number of seconds to sleep.

sleep 5


Make the terminal sleep for 5000 seconds and rescue the terminal.

sleep 5000


That’s more than an hour so use Ctrl-C to break off the command.

Type the following command in the terminal:

cd


cd stands for change directory and is used for navigating the Linux file system

Which directory are you in ? To view the name of the current working directory, type

pwd


pwd stands for print working directory. You see that using cd without arguments leads you to your home directory, on the BITS laptops this is /home/bits. Which directories are located in your home directory ? To view a list of the files and directories that are located in the current working directory, type

ls


ls stands for list and is used for listing all files and directories in the current directory. On the BITS laptops the home directory /home/bits contains a set of folders like Desktop, Documents, Downloads…

ls D*


D(star) means everything which name starts with a D

A common pattern when using the command line is changing directories using cd and then immediately typing ls to view the contents of the directory.

List the detailed content of your home directory ?

ls -l


the l in -l stands for long output. Among others, the detailed list shows a date and time indicating the last time a file was modified. The number before the date is the size of the file in bytes.

List the content of the /usr/local/bin directory ?

ls /usr/local/bin


/usr/local/bin corresponds to a directory in the file system (/), with bin a subdirectory of local and local a subdirectory of usr.

If you have to reuse a variable often then it can be helpful to create a name for a variable, especially when the variable is long. Suppose you want to work in a directory called Illumina_exp4_20042004_mapping_to_hg19_results. To avoid repeating this long name over and over you can create a variable for it, give it a short name and use that in your commands.

Name the variable folder

Use the following command: folder=Illumina_exp4_20042004_mapping_to_hg19_results

To create a new directory use the mkdir (make directory) command.

Create the folder using the newly created variable

If you want to refer to a named variable in a command you have to preceed the name by a $sign to indicate that what is following is a reference to a variable. So use the following command: mkdir${folder} The curly braces delineate the start and end of the variable name. Check if the folder is created using the ls command.

To remove a directory, use the rm (remove) command. You could use rmdir but this only works on empty folders. To remove a folder with the rm command you need to use the -r option. This stands for recursively which means it will remove the folder and its complete content.

Remove the Illumina_exp4_20042004_mapping_to_hg19_results directory.

Use the variable as an argument of the rm command:

rm -r ${folder}  Check if it’s removed using the ls command. Now navigate to the NGS folder which is located in the /home/bits/ folder. Navigate to this location. Since you want to navigate, you need to use the cd command. Since the NGS folder is located in the folder that you are currently in, you can simply give the name of the folder (NGS) as an argument: cd NGS If you want to move to a folder that’s located in another location of the file system, you have to give the full path to the folder. Go to the /usr/bin folder cd /usr/bin Go back to your home folder cd ### Manipulating files Even without a text editor, there are ways to create a file with text using the redirect operator > Create a file called test1.txt containing the text “Why do bioinformaticians work on the command line?” using echo echo "Why do bioinformaticians work on the command line?" > test1.txt  The redirect operator > takes the text output of echo and redirects its contents to a file called test1.txt Check if it worked by viewing the content of the file on the screen cat test1.txt  The name cat is short for “concatenate”. The command can be used to combine the contents of multiple files, but we use it here to dump the content of a single file to the screen. Cat is as a “quick-and-dirty” way to view the content of a file, less is a neater way. Add the line “Because they don’t want to scare you with huge amounts of data!” to the file and check if it worked To add lines of text to a file, use the append operator »: echo "Because they don't want to scare you with huge amounts of data!" >> test1.txt cat test1.txt  The append operator » appends the text output of echo to the file test1.txt Create an empty file called test2.txt and check if it exists To create an empty file, use the touch command: touch test2.txt ls  List the names of all text files in your current directory ls *.txt  Here *.txt automatically expands to all filenames that match the pattern “any string followed by .txt”. Rename the file test2.txt to test_partII.txt using mv and check if it worked To rename a file use the mv command, short for move: mv test2.txt test_partII.txt ls *.txt  Copy the file test_partII.txt to test2.txt and check if it worked To copy a file use the cp command, short for copy: cp test_partII.txt test2.txt ls *.txt  You don’t have to type out test_partII.txt, instead you can type something like test_-Tab thereby making use of tab completion. Tab completion involves automatically completing a word if there’s only one valid match on the system. For example, if the only file starting with the letters “test_” is test_partII.txt, test_-Tab refers to test_partII.txt Especially with longer names, tab completion can save a huge amount of typing. Remove the file test_partII.txt and check if it worked To remove a file use the rm command, short for remove: rm test_partII.txt ls *.txt  Download the file called exams.txt, containing the results of the spelling and maths exams of all 10-year olds of a school, into your home folder. Use wget to download the file from http://data.bits.vib.be/pub/trainingen/NGSIntro/exams.txt Download the file. wget http://data.bits.vib.be/pub/trainingen/NGSIntro/exams.txt  Show the first 10 lines of the file. head exams.txt  Two complementary commands for inspecting files are head and tail, which allow to view the beginning (head) and end (tail) of a file. The head command shows the first 10 lines of a file. Show the last 10 lines of the file. Similarly, tail shows the last 10 lines of a file. tail exams.txt  Open the manual of head to check out the options of head. Learn how to look at the first n lines of the file. Save the first 30 lines of exams.txt in a file called test.txt head -n 30 exams.txt > test.txt  Look at test.txt using the less command less test.txt  There are many commands to look at the full content of a file. The oldest of these programs is called more, and the more recent and powerful variant is called less. Less lets you navigate through the file in several ways, such as moving one line up or down with the arrow keys, pressing space bar to move a page down… Perhaps the most powerful aspect of less is the forward slash key /, which lets you search through the file from beginning to end. Search for Jasper in test.txt The way to do this in less is to type /Jasper The last three essential less commands are G to move to the end of the file and 1G to move back to the beginning. To quit less, type q (for quit). Look at the last 10 lines of the first 20 lines of exams.txt head -n 20 exams.txt | tail  The command runs head -n 20 exams.txt and then pipes the result through tail using the pipe symbol | The reason the pipe works is that the tail command, in addition to taking a filename as an argument, can take input from “standard in”, which in this case is the output of the command before the pipe. The tail program takes this input and processes it the same way it processes a file. ### Running tools Bioinformatics tools are just commands on the commands line. You use them in exactly the same way as all the commands we have run up to now, by defining options and arguments. A list of options and arguments can be found in the help file. #### Installing and running sl We have seen ls the list command and use it frequently to view the contents of a folder but because of miss-typing sometimes you would result in sl, how about getting a little fun in the terminal and not command not found. This is a general linux command, you can install it from a repository. Install sl For installing you need superuser privileges ! sudo apt-get install sl  Find out in the manual what sl stands for man sl  You can find the solution in the description section of the manual. run the command sl  :o) Try out some of the options !! #### Running blastp In the folder /home/bits/Linux/ you find a file called [http://data.bits.vib.be/pub/trainingen/NGSIntro/sprot.fasta sprot.fasta] containing a set of protein sequences. We will use this file as a database for blast. The query sequence is the following: MLLFAPCGCNNLIVEIGQRCCRFSCKNTPCPMVHNITAKVTNRTKYPKHLKEVYDVLGGSAAWE  Create a fasta file containing the query sequence using echo called seq.fasta echo ">query seq" > seq.fasta cat seq.fasta echo MLLFAPCGCNNLIVEIGQRCCRFSCKNTPCPMVHNITAKVTNRTKYPKHLKEVYDVLGGSAAWE >> seq.fasta cat seq.fasta  Blast can be done via the [https://blast.ncbi.nlm.nih.gov/Blast.cgi blast website], but you can also download the blast tool and run it locally (on your computer) via the command line. For instance if you want to blast against you own database of sequences, you have to do it locally. Blast has been installed on the bits laptops. First you have transform your own database (the sprot.fasta file in our case) into a database that can be searched by blast using the makeblastdb command. Look at the help file of makeblastdb and find the options to define the input fasta file and the database type makeblastdb -help  You have to define the input fasta file using the -in option and the type of sequences using the -dbtype option Create the blast database makeblastdb -in sprot.fasta -dbtype prot  Now you can perform a blastp search using the blastp command. Write the results to a tabular text file with comments called output.txt Look at the help file of blastp and find the options to define input, database, output and output format blastp -help  You need the -query, the -db, the -out and the -outfmt option Perform the blast and open the results with less blastp -query seq.fasta -db sprot.fasta -out output.txt -outfmt 7 less output.txt  #### Running cutadapt In this exercise we’ll do some real NGS analysis on the SRR074262.fastq file that is stored in folder /home/bits/NGS/Intro. Go to this folder and look at the 10 first lines of the file. cd /home/bits/NGS/Intro head SRR0474262.fastq  This data sets contain a high number of adapter sequences. These are reads that consist solely or partly of adapter sequence. You have to remove this adapter contamination using command line tools like [https://code.google.com/p/cutadapt/ cutadapt]. This tool is installed on the bits laptops. It is not a regular bash command (it’s a python program) so it doesn’t have a manual but it does have a help file. Check the help file of cutadapt for the option to define the adapter sequence and trim at the 3’ends of the reads. To open the cutadapt help file type: cutadapt -h  The -a option trims adapter sequences at the 3’ end of the reads. At the top of the help file you see that the standard usage of the command is: cutadapt -a ADAPTER -o output.fastq input.fastq  The sequence of the adapter is GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAA Trim the adapter sequence and store the trimmed sequences in a file called SRR074262trim.fastq Go to the folder where the input file is located and type: cutadapt -a GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAA -o SRR074262trim.fastq SRR074262.fastq  Look at the first lines of the trimmed file Go to the folder where the input file is located and type: head SRR0474262trim.fastq  #### Running Picard The trimmed fastq file is subsequently mapped resulting in a bam file that you can download from http://data.bits.vib.be/pub/trainingen/NGSIntro/1271011_reads.pair.1.list.accepted_hits.bam Download the file via the command line wget http://data.bits.vib.be/pub/trainingen/NGSIntro/1271011_reads.pair.1.list.accepted_hits.bam  Rename the file SRR074262.bam Remember to use tab autocompletion ! mv 1271011_reads.pair.1.list.accepted_hits.bam SRR074262.bam  This is a raw unsorted bam file, if we want to visualize the mapping results in IGV, we need to sort and index the file. We can do the sorting using one of [http://broadinstitute.github.io/picard/ the Picard tools], called SortSam. Picard can be downloaded from https://github.com/broadinstitute/picard/releases/download/2.8.2/picard-2.8.2.jar Download the file Remember to use tab autocompletion ! wget https://github.com/broadinstitute/picard/releases/download/2.8.2/picard-2.8.2.jar ll  For the tools to run properly, you must have Java 1.8 installed. To check your java version run the following command: java -version  Running Java tools from the command line requires a special syntax: you have to start the command with java and then the name of the java tool and its options and arguments. Java jar-files are archives of multiple java files (similar to tar archives of multiple regular files). They require an even more elaborate syntax. You have to start the command with java -jar and then the name of the jar file and its options and arguments. As you can see the picard tools come as a jar-file. Test the installation by opening the help file java -jar picard-2.8.2.jar -h  Bam files are enormous files that are hard to search through. The order of the reads in a bam file is the same as in the original fastq file. However, if you want to visualize the mapping results or if you want to calculate mapping statistics it’s much more efficient to sort the reads according to genomic location. This can be achieved with the SortSam tool. Look in [https://broadinstitute.github.io/picard/command-line-overview.html the picard documentation] for the SortSam tool. Sort the bam file to SRR074262sorted.bam Remember to use tab autocompletion ! java -jar picard-2.8.2.jar SortSam \ I=SRR074262.bam \ O=SRR074262sorted.bam \ SORT_ORDER=coordinate  Bam files contain duplicate reads unless you removed them during the quality control step. MarkDuplicates locates and tags duplicate reads in a bam or sam file. Duplicate reads originate from the same fragment and were typically introduced during library construction using PCR. Duplicate reads can also result from a single cluster on the flow cell, incorrectly detected as multiple clusters by the optical sensor of the sequencing instrument. MarkDuplicates compares sequences of reads and detects duplicates. The tool’s output is a new SAM or BAM file, in which duplicates have been identified in the SAM flags field. If needed, duplicates can be removed using the REMOVE_DUPLICATE and REMOVE_SEQUENCING_DUPLICATES options. (See [https://broadinstitute.github.io/picard/command-line-overview.html#MarkDuplicates the Picard documentation] for more details). Remove duplicates from the sorted bam file Remember to use tab autocompletion ! java -jar picard.jar MarkDuplicates \ I=SRR074262sorted.bam \ O=SRR074262sortednodup.bam \ M=marked_dup_metrics.txt \ REMOVE_DUPLICATES=true  For visualization and easy access you can build an index to the bam file using BuildBamIndex. Look in [https://broadinstitute.github.io/picard/command-line-overview.html the picard documentation] for the BuildBam Index tool. Build the bai file for SRR074262sortednodup.bam Remember to use tab autocompletion ! java -jar picard-2.8.2.jar BuildBamIndex \ I=SRR074262sortednodup.bam  Check if the files were generated. ### File compression Compress the SRR074262.bam file to .gz format Remember to use tab autocompletion ! gzip SRR074262.bam ll  and unzip it again Remember to use tab autocompletion ! gunzip SRR074262.bam.gz ll  ### Writing scripts #### Writing and executing bash scripts We are going to make additions to the bash script you find below: #this program pretends to hack sites !Define a variable str equal to " 0 1 23 45 6 789" clear !Print to screen: "hacking www.bits.vib.be" !Do nothing for 2 seconds !Print to screen: "Server hacking module is loading" !Do nothing for 2 seconds !Print to screen: "Hack module is starting in 2 seconds" !Do nothing for 1 second !Print to screen: "1 second" !Do nothing for 1 second ping -c 3 www.bits.vib.be !Do nothing for 1 second netstat !Do nothing for 1 second for i in {1..1000} do number1=$RANDOM
let "number1 %= ${#str}" number2=$RANDOM
let "number2 %=4"
!Print to screen without newlines and with backslash escapes: "\033[1m${str:number1:1}\033[0m" done !Print to screen: "453572345763425834756376534" !Do nothing for 3 seconds !Print to screen: "www.bits.vib.be succesfully hacked!" !Print to screen: "PASSWORD ACCEPTED: token is 453572345763425834756376534"  Open gedit and paste the code. Replace all lines that start with ! by the appropriate command #this program pretends to hack sites str=" 0 1 23 45 6 789" clear echo "hacking www.bits.vib.be" sleep 2 echo "Server hacking module is loading" sleep 2 echo "Hack module is starting in 2 seconds" sleep 1 echo "1 second" sleep 1 ping -c 3 www.bits.vib.be sleep 2 netstat sleep 1 for i in {1..1000} do number1=$RANDOM
let "number1 %= ${#str}" number2=$RANDOM
let "number2 %=4"
echo -n -e "\033[1m${str:number1:1}\033[0m" done echo "453572345763425834756376534" sleep 3 echo "www.bits.vib.be succesfully hacked!" echo "PASSWORD ACCEPTED: token is 453572345763425834756376534" . > Add a shebang line to the top of the script  #!/usr/bin/env bash #this program pretends to hack sites str=” 0 1 23 45 6 789” clear … Save the script as HackIt.sh > If necessary make executable  chmod 755 HackIt.sh > Run the script  bash HackIt.sh What if you want to "hack" another website ? The easiest way to do allow for this is to enable to give the url as an argument of the bash command so that's what we'll do. Reopen the file in gedit > Replace www.bits.vib.be by$1



#!/usr/bin/env bash #this program pretends to hack sites str=” 0 1 23 45 6 789” clear echo “hacking $1” sleep 2 echo “Server hacking module is loading” sleep 2 echo “Hack module is starting in 2 seconds” sleep 1 echo “1 second” sleep 1 ping -c 3$1 sleep 2 netstat sleep 1 for i in {1..1000} do number1=$RANDOM let “number1 %=${#str}” number2=$RANDOM let “number2 %=4” echo -n -e “\033[1m${str:number1:1}\033[0m” done echo “453572345763425834756376534” sleep 3 echo “$1 succesfully hacked!” echo “PASSWORD ACCEPTED: token is 453572345763425834756376534” > Save and run the script again now giving www.kuleuven.be as an argument  bash HackIt.sh www.kuleuven.be $1 refers to the first argument of the command. If you have two arguments you use $1 and$2 to represent them.

#### Writing and executing Perl scripts

We are going to create and the perl script you find below:


#This program predicts if a sequence is protein, nucleic acid or rubbish $seq =$ARGV[0]; if ($seq =~ /[JO]/) { print “is not a sequence, first illegal character is$&\n”; } elsif ($seq =~ /[EFILPQZ]/) { print “is protein\n”; } else { print “is nucleic acid\n”; }  Open gedit and paste the code. > Add a shebang line to the top of the script  #!/usr/bin/env perl #This program predicts if a sequence is protein, nucleic acid or rubbish$seq = $ARGV[0]; if ($seq =~ /[JO]/) { …

Save the script as SeqIt.pl

> If necessary make executable



chmod 755 SeqIt.pl

> Run the script using your first name in capitals as an argument



perl SeqIt.pl JANICK


#### Writing and executing Python scripts

We are going to make additions to the python script you find below:


#This program counts the number of amino acids in a protein sequence !Define variable mySequence equal to “SFTMHGTPVVNQVKVLTESNRISHHKILAIVGTAESNSEHPLGTAITKYCKQELDTETLGTCIDFQVVPGCGI” !Create a set myUniqueAminoAcids out of mySequence for aaCode in myUniqueAminoAcids: !Print to screen, use format to fill in the values: “Amino acid {} occurs {} times.”


Open gedit and paste the code.

> Replace all lines that start with ! by the appropriate command



#This program counts the number of amino acids in a protein sequence mySequence = “SFTMHGTPVVNQVKVLTESNRISHHKILAIVGTAESNSEHPLGTAITKYCKQELDTETLGTCIDFQVVPGCGI” myUniqueAminoAcids = set(mySequence) for aaCode in myUniqueAminoAcids: print(“Amino acid {} occurs {} times.”.format(aaCode,mySequence.count(aaCode)))

> Add a shebang line to the top of the script



#!/usr/bin/env python #This program counts the number of amino acids in a protein sequence mySequence = “SFTMHGTPVVNQVKVLTESNRISHHKILAIVGTAESNSEHPLGTAITKYCKQELDTETLGTCIDFQVVPGCGI” myUniqueAminoAcids = set(mySequence) …

Save the script as CountIt.py

> If necessary make executable



chmod 755 CountIt.py

> Run the script



python CountIt.py

What if you want to "count" another protein ? The easiest way to do allow for this is to enable to give the sequence as an argument of the python command so that's what we'll do.

Reopen the file in gedit

> Adjust the code to read the first argument of the python command using the sys library



!#/usr/bin/env python #This program counts the number of amino acids in a protein sequence import sys mySequence = sys.argv[1] myUniqueAminoAcids = set(mySequence) for aaCode in myUniqueAminoAcids: print(“Amino acid {} occurs {} times.”.format(aaCode,mySequence.count(aaCode)))

> Save and run the script again now giving QWEERTIPSDFFFGHKKKKLLLLLLLLLLLLLL as an argument



python CountIt.py QWEERTIPSDFFFGHKKKKLLLLLLLLLLLLLL


sys.argv[1] refers to the first argument of the command. If you have two arguments you use sys.argv[1] and sys.argv[2] to represent them.

#### Installing and using Python tools

Installing Python-based tools is not done with apt-get, instead the comand pip is used. If pip is not yet installed, the terminal will show an error message saying that pip is currently not installed. You can install pip using apt-get.

As an example we will install Biopython, a Python library for bioinformatics. See [http://biopython.org/wiki/Download the documentation] for more details.

> Install biopython

You need superuser privileges for this


sudo pip install biopython

We will write a small python script to check if Biopython was successfully installed. In the folder /home/bits/Linux/ you find a file called [http://data.bits.vib.be/pub/trainingen/NGSIntro/sprot.fasta sprot.fasta] containing a set of protein sequences that we will use as input. Move to the folder containing the file.

We will use SeqIO module of Biopython to parse the fasta file with the protein sequences. Check out [http://biopython.org/wiki/SeqIO the tutorial of the module].



!Import the SeqIO module of the Bio library !For every record in the sprot.fasta file do: !Print the id of the seq_record to the screen !Print the length of the sequence to the screen


Open gedit and paste the code.

> Replace all lines that start with ! by the appropriate command



from Bio import SeqIO for seq_record in SeqIO.parse(“sprot.fasta”,”fasta”): print(seq_record.id) print(len(seq_record))

> Add a shebang line to the top of the script



#!/usr/bin/env python from Bio import SeqIO for seq_record in SeqIO.parse(“sprot.fasta”,”fasta”): …

Save the script as ParseIt.py in the folder that contains the input file.

> If necessary make executable



chmod 755 ParseIt.py

> Run the script



python ParseIt.py


### Compressing and decompressing files

Some files or tools come in **.zip** format, how to decompress them ?

In the **/usr/bin/tools** folder you can find the zipped version of the FastQC tool. To unzip it, you have to use the **unzip** command.

The **/usr/bin/** folder belongs to the root user, not to the bits user. Therefore only root is allowed to do manipulations in this folder. Switch to root using the **su** command or type **sudo** in front of your commands. The system will ask for the password: bitstraining on the BITS laptops.

> Decompress the FastQC tool with unzip.
First look at the unzip manual to get an idea about the working of the command.
man unzip
To unzip the file you can use the simple command: unzip name_of_the_zip_file. Remember to use tab autocompletion.
This will generate a folder called FastQC in /usr/bin/tools.
After decompression use **ls** and **cd** to take a look at the content of the newly created **FastQC** folder. You will see the fastqc command in this folder.

> Make sure that you can read, write and execute the fastqc command and that other people can read and execute it.
To see the current permissions of the command:
ls -l
The command that allows you to change the access permissions of files and directories is **chmod** (change mode). chmod has two mandatory arguments:

- A three digit number representing the access permissions you want to set. Each digit refers to a different audience:

- first digit refers to the owner of the file
- second digit refers to the group the owner belongs to
- third digit refers to all others

The numbers themselves represent the permissions:

- 7 full access: read, write and execute
- 3 write and execute
- 2 write only
- 1 execute only
- 0 no access

- The name of the file for which you want to change the access permissions

As you can see **root** is the owner of the file. This is why you need to log on as superuser (= root) to be able to change root's files.

### Sorting files

We want to sort the file exams.txt from highest to lowest score on maths.

> Sort the file based on score on maths. Write results to a file called examssort1.txt
You have to **sort** the lines in the file according to the maths score. So you want to sort the file based on the numbers in the second column: it means that you cannot use the default sort command (this will sort the lines based on the content of the first column) but you have to use an option that allows you to specify the column you wish to sort on.
When you look in the manual you see that you can use the -k option for this:


sort -k2 exams.txt

This will sort the file according to the values in the second column, but it will overwrite the original file. To save the sorted list in a new file, examssort1.txt, use the **redirect operator: >**


sort -k2 exams.txt > examssort1.txt

> Use the head command to look at the sorted file.



You can see that the sorting was not done correctly: it was done alphabetically, treating the numbers in the second column as characters, instead of numbers. This means that we are still missing an option that allows for numerical sorting.

> Sort the file numerically based on score on maths.


sort -k2 -n exams.txt > examssort1.txt head examssort1.txt

This looks a lot better, but we still have to reverse the order since we want the scores from high to low.

> Sort the file numerically from highest to lowest score on maths.
For this we need to add a third option to the **sort** command.
When you look in the manual you see that you can use the -r option for this:



sort -k2 -n -r exams.txt > examssort1.txt head examssort1.txt

> Show the results of the 10 students with the highest scores on the maths exam using a single line of commands.

This means that you have to combine the **head** command and the **sort** command from the previous exercise into one single command. Remember that you can combine commands by writing them in the order they have to be performed, so in our case first **sort** then **head**, separated by the **pipe operator: |**



sort -k2 -n -r exams.txt | head

> Show only the names of the 10 students with the highest scores on the maths exam using a single line of commands.

To leave out gender and scores you have to use the **cut** command. To specify which columns to cut you can use the -f option. Please note that the -f option specifies the column(s) that you want to retain ! As an argument you have to specify the name of the file you want to cut.
In the manual you can see that TAB is the default delimiter for the cut command. So if you have a tab-delimited text file, as in our case, you do not need to specify the delimiter. Only if you use another delimiter you need to specify it.



sort -k2 -n -r exams.txt | head | cut -f3



**The case of chromosomes and natural sorting.**
'sort' will sort chromosomes as text; adding few more parameters allows to get the sort you need.

> Write a list of human chromosomes (values: 22 to 1 X Y MT) to the screen. Use {end..begin} to define a numerical range.

Remember that you can use **echo** to print text to the screen, so to generate text. Try


echo {22..1} X Y MT

and see what happens...
You don't want to numbers next to each other in one row, you want them in a column underneath each other. This means you want to replace the blanks by end-of-lines.

> Replace blanks by end-of-lines. Use the sed command for this.

Look up the command for replacing text in the slides. Blanks are represented by **\ ** (back slash followed by a blank) and end-of-lines are represented by **\n** (back slash followed by n). To replace all blanks by an end-of-line you need to add the **g** option (see [http://sed.sourceforge.net/sed1line.txt sed tutorial] for more info). So


sed “s/\ /\n/g”

should do the replacement. Of course you need to combine the two commands using the output of echo as input in sed. Look in the slides or the cheat sheet how to do this.
However, you do not want to print the text to the screen you want to print the text to a file. Look in the slides or the cheat sheet how to do this and try to combine the three parts of the command.

> Write chromosomes as a column to a file called chroms.txt

The correct solution is:


echo {22..1} X Y MT | sed “s/\ /\n/g” > chroms.txt

The s in the sed argument refers to substitution: you want to substitute blanks by end-of-lines, it is followed by the character you want to replace (a blank or "\ "), then the character you want to replace it with (an end-of-line or "\n"), then you add g to use sed recursively, in other words to do the substitution more than once so each time a blank is encountered.
It prints the chromosome numbers as a column to the file chroms.txt

> Look at the file using the less command.



less chroms.txt

Remember to use q to leave a less page.

> Sort the chromosome file by using a simple sort. Write results to chromssort.txt



sort chroms.txt > chromssort.txt head chromssort.txt

Not good! This is a tricky problem that always comes up when you are working with chromosome numbers e.g. when sorting bam/sam files, annotation files, vcf files...

> Modify the sort command so that the sorting of the chromosomes is done in the correct way.

Most people solve it by specifying that you want sort to do natural sorting using the -V option:


sort -V chroms.txt > chromssort.txt head chromssort.txt

Nice !
Now try with chr in front.

> Create a file with values chr22 to chr1 chrX chrY chrMT into one column called chroms2.txt in one single command



echo chr{22..1} chrX chrY chrMT | sed “s/\ /\n/g” > chroms2.txt head chroms2.txt

> Sort the file into a new file called chromssort2.txt



sort -V chroms2.txt > chromssort2.txthead chroms2.txt


### Getting files from the internet

For NGS analysis you often need to download genomic sequence data from the internet. As an example we are going to download the E.coli genome sequence from the iGenomes website: ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Escherichia_coli_K_12_MG1655/NCBI/2001-10-15/Escherichia_coli_K_12_MG1655_NCBI_2001-10-15.tar.gz

Go to this folder and use the wget command to download the data:


cd /home/bits/NGS?ChIPSeq/ wget ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Escherichia_coli_K_12_MG1655/NCBI/2001-10-15/Escherichia_coli_K_12_MG1655_NCBI_2001-10-15.tar.gz ll

In the same way you can download NGS data from the internet. We are not going to actually do this because

> Decompress the file.



tar -xzvf Escherichia_coli_K_12_MG1655_NCBI_2001-10-15.tar.gz ll

This creates a new folder called Escherichia_coli_K_12_MG1655.
Go into this folder and look at the whole genome fasta sequence

> Look at the fasta sequence.
Use **cd** to navigate the folders and **head** to look at the file


cd Escherichia_coli_K_12_MG1655 ll cd NCBI ll cd 2001-10-15 ll cd Sequence ll cd WholeGenomeFasta ll head genome.fa


### Installing tools

The FastQC tool was installed by unzipping it. Most tools can be installed using the **make** command. There are many ways to install software on Linux:

- via the software manager, an application with a very easy user friendly interface
- via the **apt-get** command
- software packages written in Python are installed via the **pip install** command

These methods handle the installation and removal of software on Linux distribution in a simplified way. They fetch the software from software repositories on the internet. However, these repositories do not always contain the most up-to-date version of software packages, especially not for niche software like bioinformatics tools.

So to be on the safe side, it is recommended that you download the latest version of a tool from its website (using wget) and use **make** to install it. In that way, you have full control over the version of the tool that you are installing.

This is not true for pip. Pip does the difficult steps in the installation for you and accesses an up-to-date package repository, so Python programs can safely be installed using pip.

Download and install all packages in the **tools** folder of the **/usr/bin/** folder. This is a folder owned by root so it is a good idea to switch to superuser again.

#### Installing TopHat

In the Introduction training we use RNA-Seq reads. Mapping RNA-Seq reads is done using the TopHat tool. So we need to install the [http://ccb.jhu.edu/software/tophat/tutorial.shtml TopHat tool]. We are going to do this in the /usr/bin/NGS/ folder so we need to be superuser for this.

- Go to the [http://ccb.jhu.edu/software/tophat/tutorial.shtml TopHat website]

- Go to the terimnal
- Navigate to the /usr/bin/NGS/ folder
- Type **wget **
- Press the Shift and Insert keys simultaneously to paste the url

> Decompress the file
For decompressing a .tar.gz file you need the following command:


tar -xzvf tophat-2.1.1.Linux_x86_64.tar.gz

Remember to use tab autocompletion !

This creates a new folder called tophat-2.1...
Go into the tophat folder and type:


./tophat


If this opens the help of tophat, it means the software has been installed correctly. It does not mean that you can use the software now. Well you can but you will always have to type the commands from inside the tophat folder like we do here or provide the full path to the tophat folder. The dot slash (./) in front of the command means use the tophat **that is located in this folder**. It tells the command line where it can find the script (./ = the current directory = /usr/bin/tools/tophat-2.1.1.Linux_x86_64/).To avoid this we can create a symbolic link for tophat2 (see later).

#### Installing samtools

When you navigate to the **tophat** folder in /usr/bin/NGS/ you see that samtools is automatically installed when TopHat was installed:

If you see the samtools help page when you type


./samtools_0.1.18

it means that samtools is indeed installed

[http://wiki.bits.vib.be/index.php/Introduction_to_ChIP-Seq_analysis Installing tools for the ChIP-Seq training]

It has already been installed on the bits laptops but if you need to install it, use [http://wiki.bits.vib.be/index.php/Installing_cutadapt these instructions].

### Quality control of NGS data

#### Checking the quality of the Introduction training data using FASTQC====

In the /home/bits/NGS/Intro directory you can find a file called SRR074262.fastq (the file containing Arabidopsis RNA-Seq reads), that was used in the exercises on FastQC in Windows. FastQC is a tool that checks the quality of fastq files, containing NGS data.

We will now try to do the same FastQC analysis from command line in Linux. FastQC is a java-based tool that needs java to be able to run.

> Check if the correct version of java is installed
In command line you can check if java is installed on your laptop using the following command:


java -version

You should see something like:


ava version “1.8.0_101” Java(TM) SE Runtime Environment (build 1.8.0_101-b13) Java HotSpot(TM) 64-Bit Server VM (build 25.101-b13, mixed mode)

If you get an error then you don't have java installed. If the version listed on the first line is less than 1.5 then you will have problems running FastQC and you need to update java on your laptop.

> Run FastQC
To run FastQC as a GUI just like in Windows type:


fastqc

This opens the FastQC GUI and you could load a fastq file via the GUI to get its quality report. However, you can also use fastqc as a command via the command line.

> Open the file SRR074262.fastq to obtain the sequence of the contaminating adapter.
Check [http://wiki.bits.vib.be/index.php/Quality_control_of_NGS_data#Exercise_1:_Quality_control_of_the_data_of_the_introduction_training the exercise on FastQC in Windows] for details on the quality report that is generated
The big plus of running FastQC from command line is that command line allows you to combine and run a set of commands as a program by writing a command script.

#### Automating FASTQC analyses
If you have many FASTQ files to check you might prefer running FASTQC from command line so you can loop over your files and process the reports automatically.

> View the help files of the fastqc command
As for most commands the -h option nicely opens the help file:


fastqc -h



To run via command line you can simply specify a list of files to process:


fastqc somefile.fastq someotherfile.fastq

You can specify as many files as you like. If you don't specify any files the program will open the GUI.

However, there are a few options that might be helpful to use. Since FASTQC can process FASTQ, SAM and BAM files, it is always safer to tell him upfront which format to expect.

We will generate FASTQC reports for the two FASTQ files in the /home/bits/NGS/RNASeq/ folder.

> Decompress the files
First you have to decompress the fastq files. In the cheat sheet look up the command for decompressing a .gz file


gunzip chr22_SRR1039509_1.fastq.gz gunzip chr22_SRR1039509_2.fastq.gz

Decompression of the files results in two .fastq files that can be used as inputs generating the FASTQC reports.

> Generate the FASTQC reports for the two fastq files.
As you can see in the help file of fastqc, the -f option allows you to specify the format of the input file(s).


fastqc -f fastq chr22_SRR1039509_1.fastq chr22_SRR1039509_2.fastq

The two .html files contain the FASTQC reports and can be opened in a browser.

> Open the first report in firefox via command line



firefox chr22_SRR1039509_1_fastqc.html

By default, FastQC will create an HTML report with embedded graphs, but also a zip file containing individual graphs and additional data files containing the raw data from which the plots were drawn.

> Remove the .html and the .zip files



rm *.zip rm *.html

If you have many files you might want to use a for-loop instead of typing all file names into the command.

> Write a for-loop to process the two FASTQ files.
First go back to the folder that contains the fastqc command and make sure you are operating as superuser.
Take a close look at the syntax of the for-loop that is described in the slides. We are going to use the syntax for looping over files in a folder. Don' t forget the ***** to loop over all fastq files in the specified folder:


for file in /home/bits/NGS/RNASeq/*.fastq do fastqc -f fastq ${file} done Don't forget the **$** since file is just a variable that refers to the actual files in the folder. Write every line on a different line in the terminal.

When you go to the /home/bits/NGS/RNASeq folder you should see the same html and zip files as in the previous exercise. The two .html files contain the FASTQC reports and can be opened in a browser.
If you want to save your reports in a folder other than the folder which contains your fastq files you can specify an alternative location by using the **-o** option.

> Create a new folder called FASTQCresults



mkdir FASTQCresults

> Create a variable output. Its value is the path to the newly created folder.



output=/home/bits/NGS/RNASeq/FASTQCresults/

> Write a for-loop to analyze the quality of the fastq files and write the report to the new folder

Adjust the code of the for-loop to write the results to the newly created folder


for file in /home/bits/NGS/RNASeq/*.fastq do fastqc -f fastq -o ${output}${file} done

• SRA ID
• GEO
• ENA