Tutorial: InsideDNA: 5 Unix commands every bioinformatician should know
gravatar for shilparaopradeep
6 months ago by
shilparaopradeep110 wrote:

Bioinformaticians often have to manage large text files containing reads, sequences of genomes, alignments, genetic variants, and so on. Sometimes the files you receive do not have quite (or indeed at all) the same format required by tools used in the downstream analysis. In such cases you have to rearrange the data in your files – change the delimiters, order of columns or, quite likely, discard some unnecessary values. Manual processing of large files will take plenty of time, so bioinformaticians need some programs for formatting the files. Some of the most common formatting tasks (for example, editing the header of a sam file) can be solved easily by using special tools. However, formatting tasks are so diverse, that the entire spectrum of such tasks will hardly ever be addressed by ready-to-use tools. So, bioinformaticians often have to write small programs in order to rearrange data in a specific file in a certain way. Luckily, Terminal (or Console) provides opportunities to perform rather complex manipulations of text files through the use of very simple commands. In this tutorial we will discuss some useful recipes for bioinformatics file processing. 1. Merging multiple files into one To begin with, log in to the InsideDNA, navigate to the Files tab, and create a folder called Formatting. Upload the example files that we will use in this tutorial (click on the button Download data at the top of the page), or use your own files.

We will work with the virtual Terminal of the InsideDNA platform. To activate it, click on the Terminal button on the left.

To start work, navigate to the Formatting folder by typing the command: cd Formatting To see the folder content, use the ‘ls’ command: ls We need to unzip the file with raw data. Use unzip command for this purpose: unzip source_formatting.zip Run ls again. Now you can see the list of files that you uploaded to the Formatting folder.

Firstly, we will study a very simple but common task – merging multiple files into a single one. Genomes are often stored in the form of separate contigs or chromosomes, but sometimes it is necessary to have them in the form of a single file (for example, to use a genome as a reference). This can be done using a single command: cat *.fa > genome.fa This particular variant of command will merge all files in the current folder with the extension .fa in the into a single file named genome.fa (asterisk is equivalent to “anything”). You find the new file if you type the command ‘ls’ once again. Note that all the original .fa files are still in place.

  1. Extracting a subset of data into a new file Sometimes files contains information that is unnecessary for the next stage in processing. You may want to discard the unnecessary information in order to make files smaller (which could be essential when working with large amounts of data) or to prevent misinterpretation of extra information by downstream software. For such types of processing it is useful to study a language called “awk”. For example, to select certain columns from an input file and write them into new one, you can use a command such as this one: awk '{ print $1,$2}' /data/userXXX/Formatting/filtered_snps.vcf > /data/userXXX/Formatting/snp_positions.txt In this case awk will take the first two columns of the input file filtered_snps.vcf (don’t forget to replace XXX with your own userID, which can be found at the top of the Terminal folder) and save them into the new file snp_positions.txt. Note that you can change the order in which columns will be written into the output file: for example, if you modify from'{ print $1,$2}' to '{ print $2,$1}', the order of columns in the output file will be reversed. You can check that the new file now appears in the Formatting folder, by using the ‘ls’ command:

Awk takes a small program, marked by ‘{}’, and runs it in all lines of an input file. By default, it takes an arbitrary number of spaces or tabs as a delimiter. To show awk that your input file has another type of delimiter, for example, a comma, you can use the option –F: awk -F, '{ print $3 }' /data/userXXX/Formatting/filename.csv > /data/userXXX/Formatting/new_filename.csv What if you need to extract not columns, but rows from the file? For example, a very common task is to select variants with a special filter value from a VCF file. In this case you need to take all the lines where the required value occurs in a specified field: awk '$7 == "PASS" { print }' /data/userXXX/Formatting/filtered_snps.vcf > /data/userXXX/Formatting/pass_snps.vcf

This command will take lines that have the word PASS in the 7th column (which is the one with filter flags) and write them to the file pass_snps.vcf. With awk, you can combine requirements into more complex commands using the && (double ampersand) sign and other logical signs. For instance, this command: awk '$1 == chr3 && $7 == "PASS" { print }' /data/userXXX/Formatting/filtered_snps.vcf > /data/userXXX/Formatting/pass_snps_chr3.vcf

will take just the lines featuring SNPs on the 3rd chromosome, which have also passed all the filters (so they have corresponding values in the 1st and 7th columns of the VCF file), and write these to the file named pass_snps_chr3. 3. Changing the specified patterns Sometimes you need to replace a specific word or value in an input file with another one. For example, a very common requirement of tools, when working with alignments, is to have all the chromosome names in the format chr1, chr2, and so on. For such purposes, the “sed” language is very useful. For instance, this short command: sed 's/chromosome/chr/g' < input.txt > output.txt will replace each instance of the word “chromosome” in the input file with the letters “chr” and will write the result in the file output.txt (option g when using the sed command specifies that each instance of “chromosome” in the file should be replaced, rather than just the first occurrence in each line, which is the default for sed). Sed is commonly used to fix the headers of sam and bam files. There is a special tool called samtools reheader for this purpose, which takes bam files and a sam file with a correct header as two inputs and amend the header in a bam file. The reason for this is that bam files are human unreadable and therefore you cannot manually modify the header in a bam file. Below we are listing a sequence of commands which allows you to modify bam header. The first command we are going to use is samtools view –H. It will create a file called header.sam using the header of the input candida.bam file: idna_submit.py -t reheader -c 1 -r 0.6 -e "idna_samtools_view -H /data/userXXX/Formatting/candida.bam > /data/userXXX/Formatting/header.sam"

Note that in this case we used container idna_submit.py to send our task to the virtual machine. We gave this machine one core and 0.6 Gb of RAM, which are very small amounts of resources. Other tasks may require many more cores and much more RAM to be completed in reasonable time. You can monitor the progress of tasks and send them to virtual machines, using the Tasks tab:

When the task is finished, you will find the new file header.sam in the Formatting folder. To see what was written in this file, use the less command: cd Formatting less header.sam

As you can see, the names of chromosomes in this header have a very complex structure, for example “Ca21chr1_C_albicans_SC5314”. To exit “less” mode, press the Q key. To change chromosome names to the standard “chr1” format, we will use a sed command like the following: sed -e 's/SN:Ca21chr([0-9R])_C_albicans_SC5314/SN:chr\1/' -e 's/SN:Ca19-mtDNA/SN:chrM/' < /data/userXXX/Formatting/header.sam > /data/userXXX/Formatting/header2.sam

This command replaces all the complex chromosome names with the letters “chr”, which are followed with the number of chromosomes, or the letter R. Sed remembers this number or letter because of the specific pattern ([0-9R]), we used in the command. The program will search for this character after the sequence SN:Ca21chr and before the sequence _C_albicans_SC5314, as we specified in the first part of command; it will then replace all instances of this pattern with the symbols SN:chr and the “remembered” character. Notice that we need to put \ before each _ sign in the described pattern, in order to prevent sed from interpreting them as delimiters. We need to tackle the case of mitochondrial chromosomes separately, since these have a different structure in our source header file. The second part of the command that comes after the –e option is related to this task. After the sed run, you can use the “less” command to check that new file ‘header2.sam’ has the desired structure for chromosomes names: less header2.sam

Finally, you need to replace the header of the original bam file with the newly generated correct header. We will use samtools reheader command for this purpose: idna_submit.py -t rehead2 -c 1 -r 0.6 -e "samtools_reheader /data/userXXX/Formatting/header2.sam /data/userXXX/Formatting/candida.bam"

Note that you don’t need to apply the steps that use samtools view and samtools reheader if you need to fix a file in sam or in other human-readable formats – you can edit them directly by applying the sed command. Now you have mastered the basic techniques of formatting that are essential for processing large files containing bioinformatics data. Since fixing the format of files comprises a substantial part of a bioinformatician’s work, hopefully such techniques will help you to save plenty of time. Stay tuned!

ADD COMMENTlink modified 6 months ago by Devon Ryan60k • written 6 months ago by shilparaopradeep110

Dear Sir or Madam, I think this tutorial needs a great deal of work. There are many issues I would like to address but have decided not to spend too much of my own time. The one comment I will make is regarding the formatting of this text, which makes it very difficult for me to glean any useful information from it.

ADD REPLYlink written 4 months ago by mcc60
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1418 users visited in the last hour