Awk: conditions

Just a small working example to remember:

awk '{ if ($3 == "20" || $3 == "17") print $0}' Homo_sapiens.GRCh37.68.simple_names.refGene | less


Exome coverage and per chromosome coverage

This post is actually inspired by the discussion on the Qualimap google groups.

Assume we need to calculate mean coverage inside of the exons or any other arbitrary regions defined by an annotation file. In this post I will present 3 approaches to perform this task. The impatient can check my favourite (and hopefully the easiest) method.

After googling a little bit on the given problem it's easy to find that one can do coverage calculation using Samtools mpileup or coverageBed. The third option that I will show here will be the Qualimap tool.

So, here we go.

Using coverageBed

Given file test.bam and regions defined in file test.bed one can calculate coverage using the following command:

coverageBed -abam test.bam -b test.bed -d | awk '{c+=$8;len+=1}END{print "mappedBases=" c ,"RegionsSize=" len, "meanCoverage=" c/len}'

What's going on here? coverageBed outputs coverage at each position of each region intersecting with the BAM file; output is piped to a simple awk script which will accumulate the total length of regions and the number of mapped bases (column 8); finally awk will output mean coverage.

Unfortunately, there is a minor problem in this script. Since some of the regions could intersect, we count intersecting bases twice. It does not matter if you are interested in absolute coverage of your regions, but this could bias the mean coverage.

Using samtools mpileup

This option is quite popular. Here is the command:

samtools mpileup -l test.bed -A -Q 0 test.bam | awk 'BEGIN{regionLen=317355}{c+=$4}END{print "mappedBases=" c,"meanCoverage=" c/regionLen}'

Mpileup outputs coverage at each position of the reference intersected with region which has non-zero coverage. To compute the mean coverage one has to calculate the "efficient" length of the regions: the length including non covered bases and counting the intersections only once. In this script the efficient length is the magic number 317355.

More things to keep in mind: since mpileup is originally designed to discover SNPs it may discard a lot of reads based on the minimum coverage, read quality and other parameters. So one has to turn off the turnings by using appropriate command line arguments: -A, -Q0, -d etc. Refer to mpileup documentation for more details.

Using Qualimap:

Qualimap calculates coverage similar to mpileup, which means it calculates coverage at each position of the reference intersected with a region. Moreover it calculates the "efficient length" of the regions (taking into account non covered bases and intersection of regions) when computing mean coverage. Some other advantages of using Qualimap: it can work both with BED and GFF files and output coverage of "outside" regions. Finally, Qualimap also reports per chromosome coverage. Ok, enough of this advertising :) Here's your command:

qualimap bamqc -bam test.bam -gff test.gff

So, what's your preferable way of calculating coverage?

See ya later, folks :)

Source of power

Came across a cool website with a lot of info on Unix tools such as bash, awk, sed:
http://www.grymoire.com/Unix/

Another crazy one-liner: list of genes in GTF file

Well, here is another UNIX one-liner (thanks to my supervisor), which is really crazy:  

grep codon hg19.genes.gtf | awk 'BEGIN{FS="\t";OFS="\t"}{split($9,a,"\"");print a[4]}' | sort | uniq | grep -f - hg19.genes.gtf > hg19.genes.withCDS.gtf 

Ok, what's going on here?

Well, I want to sort the GTF file in the following way: only leave records of genes, which have the CDS regions defined. The lines defining CDS have the type "start_codon" or "stop_codon". One can write a script in Python which will parse the GTF file, collect all lines for single gene record and output them only if CDS is present. However, it is also possible to do this using awk and grep.

Have fun deciphering this one :)

A couple of hints:

• split() is an awk function to split strings
• the "-f" param for grep takes a dash, which is a pointer for standard input

BED to GFF in one line

Here you go:

awk 'BEGIN{OFS="\t"}{ print $1,"qualimap",$4,$2,$3,$5,$6,".","bed" }' CpGIslandsByTakai.human.bed > file.gff

NOTE:
Here I use on 6 fields of BED file, but you can use this script as a basis and allow more fields to be used.

Sum up a column in a tab delimited file

Suppose we have a tab-delimited file:

kokonech@ultor:~/playgrnd/densityAnalysis/density_repeats.64.LTRs.bed$ head 24h-i-input.bam.bed
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 1 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 2 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 3 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 4 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 5 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 6 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 7 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 8 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 9 0
1 83886031 83886750 ERVL-E-int 980 -1 N 187 LTRs 10 0


We want to sum up a column, let's say the 11th.
A piece of cake using Awk:

awk '{a+=$11}END{printf "%i\n",a}' 24h-i-input.bam.bed


counting the sizes of cpg-islands

Small task: we need to calculate the histogram of CpG island sizes.

The CpG islands are stored in a BED file:

kokonech@ultor:~$ head CpGIslandsByTakai.wihtNames.bed
track name="CpG islands by Takai et al." description="CpG islands by Takai et al." color=0,200,100 CpG_0 0 +
chr1 10231 11413 CpG_1 0 +
chr1 26760 27059 CpG_2 0 +
chr1 28597 29942 CpG_3 0 +
chr1 51388 52085 CpG_4 0 +
chr1 88289 88610 CpG_5 0 +
chr1 121148 121635 CpG_6 0 +
chr1 134035 134389 CpG_7 0 +
chr1 134896 135424 CpG_8 0 +
chr1 175186 175545 CpG_9 0 +


First we get the sizes using simple awk command:

awk {'print $3 - $2 + 1'} CpGIslandsByTakai.wihtNames.bed > CpGislands.sizes


Then I use a small python script to draw histogram.
The listing of the script can be found below.

 import numpy as np  
 import matplotlib.pyplot as plt  
 import matplotlib.mlab as mlab  
 import sys  
   
 x = []  
   
 file = open(sys.argv[1])  
   
 for line in file:  
   val = int(line)  
   if (val > 3000):   
     val = 3000  
   x.append ( val )  
   
 fig = plt.figure()  
 ax = fig.add_subplot(111)  
   
 # the histogram of the data  
 n, bins, patches = ax.hist(x, 50, normed=1, facecolor='green', alpha=0.75)  
   
 ax.set_xlabel('CpG island sizes')  
 ax.set_ylabel('Count')  
 ax.grid(True)  
   
 plt.show()  

Small tricks for tab-delimeted files

Working on my current project includes a lot of manipulation with tab-delimited files (SAM,GFF, etc) Here are some nice tricks to remember.

Output only second field of a delimited file:
cat samout | cut -f 2

The same with awk:
awk '{print $2}'

Not related to tab-delimited files but still very useful :)

Compare two files based on their content:
cat file1 file2 | sort | uniq -u