Hdi2u rendertotsv exercise

Introduction

• From a reference sequence annotation of an organism with 14 chromosomes, coordinates for a gene family KIR have been isolated.
• A mapping of reads to this region, and a bedtools intersect operation was carried out, after which coordinates were converted to lengths.
• File: a list of reference sequence genes and their sizes: kir_ref_sizes.tsv
• 3 files: Read intersects filtered for 3 different coverage levels:

kirreads_covlev1.tsv
kirreads_covlev2.tsv
kirreads_covlev3.tsv

Actually we can refer to these three almost as if they were one file:

kirreads_covlev{1,2,3}.tsv

Key problem

• We want to view this in a spreadsheet, with reference gene lengths and read coverage length in columns, so we can visually calculate percentage coverage.
• Some genes in the higher coverage level reads will have been filtered, i.e. each gene has row: absent rows
• Of course we could manually re-arrange in a spreadsheet, or mentally build up some command lines
• We can lean on grep, it will print out occurences in all files of a certain string (gene name)

Look at our data

head kir_ref_sizes.tsv kirreads_covlev{1,2,3}.tsv

Try:

grep chr01_0290 kir_ref_sizes.tsv kirreads_covlev{1,2,3}.tsv
grep chr02_1410 kir_ref_sizes.tsv kirreads_covlev{1,2,3}.tsv
wc -l kirreads_covlev{1,2,3}.tsv

• Notice how the results come out in separate row, and often not the same number of rows.
• The decreasing coverage sizes show us the filter at work, so we know we are working with the right data at least.

Start with a very basic for-loop

• We can start by picking out the gene names from the reference
• But it more than just a list of names, it has values too which must be cut out by using cut

for i in $(cut -f1 geneskir_ref_sizes.tsv); do echo $i;done

• We prefer a cleaner gene name, and perhaps we not interested in seeing the whole list

for i in $(cut -f1 kir_ref_sizes.tsv); do echo ${i%.*};done | head

${i%.*} is a brace expansion: it strips the variable from the last dot onwards.

• Fine, so we're ready now to use grep instead of echo here.
• We need to feed it the reference file and 3 read coverage files:

for i in $(cut -f1 kir_ref_sizes.tsv); \
   do grep ${i%.*} kir_ref_sizes.tsv kirreads_covlev{1,2,3}.tsv; \
   echo;done | head

• Note above how we inserted an echo which gives us a blank line, this will be useful

• Remove | head and add > kirgenecov.txt to capture into a file.
• Have a look at kirgenecov.txt: how well it looking?
• It's not great on the eyes, but we can indeed verify that all the genes are matched up properly, which is useful.

Crank up the complexity: conditionals

• We only want the gene names for the first column
• We need to identify which column we are on and choose an output depending on which column
• We use a counter called C which we increment for each column (C++) and set to zero on a blank line ((/^$))
• We clean up the gene name by choosing 15 characters starting from character 19: substr($1,19,15)
• printf is formatted print, it is practially available in all languages and gives us high control over how our output looks

awk 'BEGIN{C=0}{if(/^$/) {printf "\n"; C=0} \
   else {if(C==0) {N=substr($1,19,15);printf "%s\t%i\t",N,$2; C++} \
   else {printf "%i\t",$2;C++}}}' \
   kirgenecov.txt > kirgenecov.tsv

• Have a look at this file. It shoulld be ready for opening in a spreasheet.
• Sure, we have to add a header row. That's not hard.
• What about blank cells? Cn do it manuall, but there are other much bigger gene families .. we need the zeros!

Inputting zeros automatically

• Adding zeros so in columns 3 and 4 where no values exist
• We can use awk again, it has a variable NF, number of fields, essentially the number of columns in each row
• We can see that when there is a row missing the final column (NF==4), we need an extra zero
• When there are two rows missing the final column (NF==3), we need two extra zeros
• printf "%s",$0 allows us print the whole line without the final newline

awk '{printf "%s",$0; \
   if(NF==5) printf "\n"; else if(NF==4) printf "0\n"; else if(NF==3) printf "0\t0\n"}' \
   kirgenecov.tsv > kirgenecov0.tsv

• And that's it! All the spreadsheet needs now is a header file, adn we can play around with the number in traditional Excel fashion.
• Final question if there's time: Let's say we wanted comma-separated values, what easy step could we take (in Command-line Navigation section).