use getline to capture system command output in awk

I just learnt this today: awk also has its pipe (|) and getline, just like unix. If I want to call a system command in awk and capture its output (Note: system() won't work as it only return the exit status), I can use pipe the output to getline. 

For example,

$cat > test.txt

aa bb cc

11 22 33

44 55 cc

$awk 'BEGIN{cmd="grep cc test.txt | sort -k1,1 | head -n1 | cut -f2 -d\" \""; cmd | getline a; print a}'

55

or 

$awk 'BEGIN{cmd="grep cc test.txt | sort -k1,1 | head -n1 | cut -f2 -d\" \""; system(cmd);}'

55

Note that if only use cmd, it won't print out anything, because cmd itself won't switch to console (unlike system). 

$awk 'BEGIN{cmd="grep cc test.txt | sort -k1,1 | head -n1 | cut -f2 -d\" \""; cmd;}'

If you want to output the multiple lines and process them in awk, you can do

$awk 'BEGIN{cmd="grep cc test.txt | sort -k1,1 | cut -f2 -d\" \""; while( (cmd | getline a) >0) print a;}'

55

bb

Reference: http://stackoverflow.com/questions/1960895/awk-assigning-system-commands-output-to-variable

median filter in AWK

Here is what median filter does from wikipedia:

To demonstrate, using a window size of three with one entry immediately preceding and following each entry, a median filter will be applied to the following simple 1D signal:
x = [2 80 6 3]
So, the median filtered output signal y will be:
y[1] = Median[2 2 80] = 2
y[2] = Median[2 80 6] = Median[2 6 80] = 6
y[3] = Median[80 6 3] = Median[3 6 80] = 6
y[4] = Median[6 3 3] = Median[3 3 6] = 3
i.e. y = [2 6 6 3].
Note that, in the example above, because there is no entry preceding the first value, the first value is repeated, as with the last value, to obtain enough entries to fill the window. This is one way of handling missing window entries at the boundaries of the signal.

Here is my awk code to implement this:

#!/bin/awk -f
# awk script to filter noise by sliding a window and taking mean or median per window
# Authos: Xianjun Dong
# Date: 2015-06-23
# Usage: _filter.awk values.txt
# bigWigSummary input.bigwig chr10 101130293 101131543 104 | _filter.awk -vW=5 -vtype=median
BEGIN{
  if(W=="") W=5; 
  if(type=="") type="median";
  half=(W-1)/2;
} 
{
  for(i=1;i<=NF;i++) {
    array[half+1]=$i
    for(k=1;k<=half;k++){
      array[half+1-k]=(i<=k)?$1:$(i-k);
      array[half+1+k]=((i+k)>NF)?$NF:$(i+k);
    }
    if(type=="median") {asort(array); x=array[half+1];}
    if(type=="mean") {x=0; for(j=1;j<=W;j++) x+=array[j]; x=x/W;}
    printf("%s\t", x);
  }
  print "";
}

Using one line command as input for LSF bsub

In a simple case, you can use bsub command arguments to submit your command job to LSF cluster.

If you have a complicated script with many commands,  you can save into a lsf script (including the shell pathname in the first line) and then submit that script to LSF cluster, e.g.  bsub yourscript arguments

In your script, you wrote something like this:

#!/bin/bash
myFirstArgument = $1

Here I found I can also use pipe to connect multiple commands into one line and simply quote them as one command and works in bsub. Here is an example:

bsub "echo -ne 'ab\tcss' | awk '{print \$2}'"

So far, I found I have to add "\" (backslash) to escape the special character, such as $ in awk. Wondering there might be a way in bsub options to set this.

Trimmed mean and median in AWK

This can be easily done in R, but sometime you want to get it in scripting language like awk or perl in order to process the big data line by line.

Here is the code snip:

# for median
# Note: thanks to anonymous reply below, put c, j in arguments to define that they are local variables.

function median(v, c, j) 

{ 

    c=asort(v,j); 

    if (c % 2) return j[(c+1)/2]; 

    else return (j[c/2+1]+j[c/2])/2; 

}

# for trimmed mean (where p is the percentage of data to be trimmed)
function trimmedMean(v, p, c, j) 

{ 

    c=asort(v,j); 

    a=int(c*p);

    for(i=a+1;i<=(c-a);i++) s+=j[i];

    return s/(c-2*a); 

}

To use it, for example if we want to generate a merged bigwig track for a list of samples, we can take median value of all samples at each genomic position, here is it:

unionBedGraphs -i `ls *normalized.bedGraph` | awk 'function median(v) {c=asort(v,j); if (c % 2) return j[(c+1)/2]; else return (j[c/2+1]+j[c/2])/2.0; } {OFS="\t"; n=1; for(i=4;i<=NF;i++) S[n++]=$i; print $1,$2,$3, median(S)}' > AllSamples.median.normalized.bedGraph
bedGraphToBigWig AllSamples.median.normalized.bedGraph ChromInfo.txt AllSamples.median.normalized.bw

faster way to get number of reads species

Let's call the unique reads sequences as 'species' in this topic.  So, an easy way to get the number of reads species for a RNAseq (for example) sam file can be

cut -f10 alignment.sam | sort -u | wc -l

But it's not the faster one if the sam file is large. We can improve it by not sorting it (which is unnecessary for this task). Here is it:

cut -f10 alignment.sam | awk '{r[$1]++;}END{for(i in r)j++; print "number of species:", j;}'
or

awk '{r[$10]++;}END{for(i in r)j++; print "number of species:", j;}' alignment.sam

The point is, you don't really need to "sort" the key, but instead to "group" the lines by the key. For example, when using groupBy in bedtools, it requires the input is pre-sorted by the key which is to be grouped. To do that in a quick way (esp. important for large file), we can use array of array for awk or hash of array in perl. Here is the hash of array solution in perl for "sort -k10" of a large file:

perl -e 'while (<>) {$l=$_; @a=split("\t", $l); push(@{$HoA{$a[17]}}, $l);}{foreach $i (sort keys %HoA) {print join("", @{$HoA{$i}});}}'

btw, I found this awk array of array does not work:

awk '{r[$10][length(r[$10])+1]=$0;}END{for(i in r) for (j in r[i]) print r[i][j];}' alignment.sam

get intron, UTR, CDS from bed12 format

Input: gene annotation in BED12 format (If your input is GTF/GFF, you want to convert into bed12 ahead, see my another post for howto)
Output: intron/UTR/CDS in bed12 format (which can be further converted into bed6 via "bedtools bed12tobed6" command)

# Introns (if any) of gene annotation in BED12 format
cat annotation.bed | awk '{OFS="\t";split($11,a,","); split($12,b,","); A=""; B=""; for(i=1;i<length(a)-1;i++) {A=A""(b[i+1]-b[i]-a[i])",";B=B""(b[i]+a[i]-(b[1]+a[1]))",";} if($10>1) print $1,$2+a[1], $3-a[length(a)-1], $4,$5,$6,$2+a[1], $3-a[length(a)-1],$9,$10-1,A,B;}'

# 5´ UTR (if any) of gene annotation in BED12 format
cat annotation.bed | awk '{OFS="\t";split($11,a,","); split($12,b,","); A=""; B=""; if($7==$8) next; if($6=="+" && $2<$7) {for(i=1;i<length(a);i++) if(($2+b[i]+a[i])<=$7) {A=A""a[i]",";B=B""b[i]",";} else {A=A""($7-$2-b[i])",";B=B""b[i]","; break; } print $1,$2,$7,$4,$5,$6,$2,$7,$9,i,A,B;} if($6=="-" && $8<$3) {for(i=length(a)-1;i>0;i--) if(($2+b[i])>=$8) {A=a[i]","A;B=($2+b[i]-$8)","B;} else {A=($2+b[i]+a[i]-$8)","A;B=0","B; break; } print $1,$8,$3,$4,$5,$6,$8,$3,$9,length(a)-i,A,B;}}'
## Update: the code crossed above did not consider the special case that start codon is spliced (e.g. in intron within in the start codon). 
awk '{OFS="\t";split($11,blockSizes,","); split($12,blockStarts,","); blockCount=$10;A=""; B=""; if($7==$8) next;N=0;if($6=="+" && $2<$7) {start=$2;end=$7; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i]+blockSizes[i])<=$7) {A=A""blockSizes[i]",";B=B""blockStarts[i]","; end=($2+blockStarts[i]+blockSizes[i]); N++;} else { if(($2+blockStarts[i])<$7) {A=A""($7-$2-blockStarts[i])",";B=B""blockStarts[i]","; N++; end=$7;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;} if($6=="-" && $8<$3) {start=$8;end=$3; for(i=blockCount;i>0;i--) if(($2+blockStarts[i])>=$8) {A=blockSizes[i]","A;B=($2+blockStarts[i]-$8)","B;start=($2+blockStarts[i]); N++;} else {if(($2+blockStarts[i]+blockSizes[i])>$8) {A=($2+blockStarts[i]+blockSizes[i]-$8)","A;B=0","B; N++; start=$8;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;}}'

## Update (2017-Jan-28, see comment below)
awk '{OFS="\t";split($11,blockSizes,","); split($12,blockStarts,","); blockCount=$10;A=""; B=""; if($7==$8) next;N=0;if($6=="+" && $2<$7) {start=$2;end=$7; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i]+blockSizes[i])<=$7) {A=A""blockSizes[i]",";B=B""blockStarts[i]","; end=($2+blockStarts[i]+blockSizes[i]); N++;} else { if(($2+blockStarts[i])<$7) {A=A""($7-$2-blockStarts[i])",";B=B""blockStarts[i]","; N++; end=$7;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;} if($6=="-" && $8<$3) {start=$8;end=$3; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i])>=$8) {if(start==0) {A=blockSizes[i];B=0; start=$2+blockStarts[i];} else {A=A","blockSizes[i];B=B","($2+blockStarts[i]-start);} N++;} else { if(($2+blockStarts[i]+blockSizes[i])>$8) { A=($2+blockStarts[i]+blockSizes[i]-$8);B=0; N++; start=$8;} if(($2+blockStarts[i]+blockSizes[i])==$8) start=0;} print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;}}'

# CDS (if any) of gene annotation in BED12 format
grep "protein_coding\.protein_coding" annotation.bed | awk '{OFS="\t";split($11,a,","); split($12,b,","); A=""; B=""; if($7==$8) next; j=0; for(i=1;i<length(a);i++) if(($2+b[i]+a[i])>$7 && ($2+b[i])<$8) {j++; start=$2+b[i]-$7; size=a[i]; if(($2+b[i])<=$7) {start=0;size=size-($7-($2+b[i]));} if(($2+a[i]+b[i])>=$8) {size=size-($2+a[i]+b[i]-$8);} A=A""size",";B=B""start",";} print $1,$7,$8,$4,$5,$6,$7,$8,$9,j,A,B;}'

### Note1: The thickStart and thickEnd in the bed12 format don't always indicate CDS. "When there is no thick part, thickStart and thickEnd are usually set to the chromStart position." So, we cannot use bed12 to infer CDS (or coding exons), esp. for lincRNA. The correct way is to grep all CDS lines from the GTF file, or directly using UCSC Table Browser to download coding exons.

### Note2: The "Polymorphic pseudogene" in GENCODE can also have a protein-coding transcript, even though the gene itself is classified as a pseudogene (YES, polymorphic pseudogenes are also pseudogenes, they "are coding gene that are pseudogenic due to the presence of a polymorphic premature stop codon in the reference genome" (http://www.genomebiology.com/2012/13/9/R51). For more details, see https://gencodegenes.wordpress.com/toolbox/.

# 3´ UTR (if any) of gene annotation in BED12 format
cat annotation.bed | awk '{OFS="\t";split($11,blockSizes,","); split($12,blockStarts,","); blockCount=$10;A=""; B=""; if($7==$8) next;N=0;if($6=="-" && $2<$7) {start=$2;end=$7; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i]+blockSizes[i])<=$7) {A=A""blockSizes[i]",";B=B""blockStarts[i]","; end=($2+blockStarts[i]+blockSizes[i]); N++;} else { if(($2+blockStarts[i])<$7) {A=A""($7-$2-blockStarts[i])",";B=B""blockStarts[i]","; N++; end=$7;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;} if($6=="+" && $8<$3) {start=$8;end=$3; for(i=blockCount;i>0;i--) if(($2+blockStarts[i])>=$8) {A=blockSizes[i]","A;B=($2+blockStarts[i]-$8)","B;start=($2+blockStarts[i]); N++;} else {if(($2+blockStarts[i]+blockSizes[i])>$8) {A=($2+blockStarts[i]+blockSizes[i]-$8)","A;B=0","B; N++; start=$8;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;}}'

cat annotation.bed | awk '{OFS="\t";split($11,blockSizes,","); split($12,blockStarts,","); blockCount=$10;A=""; B=""; if($7==$8) next;N=0;if($6=="-" && $2<$7) {start=$2;end=$7; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i]+blockSizes[i])<=$7) {A=A""blockSizes[i]",";B=B""blockStarts[i]","; end=($2+blockStarts[i]+blockSizes[i]); N++;} else { if(($2+blockStarts[i])<$7) {A=A""($7-$2-blockStarts[i])",";B=B""blockStarts[i]","; N++; end=$7;} break; } print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;} if($6=="+" && $8<$3) {start=$8;end=$3; for(i=1;i<=blockCount;i++) if(($2+blockStarts[i])>$8) { if(start==0) {A=blockSizes[i];B=0; start=$2+blockStarts[i];} else {A=A","blockSizes[i];B=B","($2+blockStarts[i]-start);} N++; } else { if(($2+blockStarts[i]+blockSizes[i])>$8) { A=($2+blockStarts[i]+blockSizes[i]-$8);B=0; N++; start=$8;} if(($2+blockStarts[i]+blockSizes[i])==$8) start=0;} print $1,start,end,$4,$5,$6,start,end,$9,N,A,B;}}'

Update: Thanks to Heather's comment below, when the thickEnd (end of CDS) is at the end of an exon (or thickStart at the first nt of an exon), it will generate an invalid BED12 format. It's fixed now.  (2017-Jan-28)

BTW, the above codes are integrated into a neat script in Github: https://github.com/sterding/RNAseq/blob/master/bin/bed12toAnnotation.awk

Also, I fixed some bugs in gtf2bed (originally by Erik) and host here:
https://github.com/sterding/RNAseq/blob/master/bin/gtf2bed

FS does not work if putting in awk body

Just noticed this by double-checking the weird result I got:

For example , you have a text file with some empty columns in lines, like:


JUNCTION.100681___1___6___1___6___JUNCTION.100681___ ___17___ ___12

(which is tab-separated, and the 7th and 9th columns are empty)

awk '{FS="\t";OFS="\t"; print $1,$2+$7,$3+$8,$4+$9,$5+$10}' will output:

JUNCTION.100681___18___18___1___6

It takes $7 as 17 and $8 as 12, which is wrong. The right output should be:

JUNCTION.100681___1___23___1___18

It seems FS="\t" does not work properly in the body of AWK.

The right solution is:


awk 'BEGIN{FS="\t"}{OFS="\t";print $1,$2+$7,$3+$8,$4+$9,$5+$10}'

or

awk -F"\t" '{OFS="\t";print $1,$2+$7,$3+$8,$4+$9,$5+$10}'

or



awk -v FS="\t" '{OFS="\t";print $1,$2+$7,$3+$8,$4+$9,$5+$10}'

I am not sure if this is a bug.
btw, I am using GNU Awk 3.1.5.

awk script to correct XS:A tag of Tophat output for strand-specific paired-end reads

It's been noticed that the current Tophat (v2.0.3) can assign wrong XS:A tag for strand-specific paired-end library, at least for some reads. Here is such an example:

HWI-ST560:74:D14ELACXX:4:1201:8827:168386 pr1 chr1 10024104 50 50M = 10020756 -3398 TGGTTCTTGAAACTGCTGGTTCAGCATCTGTGTACTAACATCAATCCCGG IJJJJJJJIIIGJJJJJJJJJJJJJJJJJJJJIJJJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:50 YT:Z:UU NH:i:1 XS:A:+
HWI-ST560:74:D14ELACXX:4:1201:8827:168386 pR2 chr1 10020756 50 36M1628N14M = 10024104 3398 GATGATTTGAAATATGAGACTTCTAAGGCATAATATTGTTTGCAGTGCAC CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIIIIJJJJIJIJHJGEC AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:50 NM:i:0 XS:A:- NH:i:1

This is a dUTP protocol where "R2/r1" FLAG for the read pair indicate the reads are a transcript on the + strand, which means both reads should be assigned as XS:A:+. However /2 is assigned to XS:A:-.

I've written an awk script to solve the problem:

#!/bin/awk -f

BEGIN{

    if(save_discrepancy_to_file!="") system("[ -e " save_discrepancy_to_file " ] && rm " save_discrepancy_to_file);

}

{

    if($1 ~ /^@/) print;

    else

    {

        for(i=1;i<=NF;i++) if($i!~/^XS/) printf("%s\t",$i); else XS0=$i;

        XS1=XS0;

        if($2~/^0x/ || $2~/^[0-9]+$/){   # FLAG in HEX or Decimal format

            if(libtype=="fr-firststrand") XS1=((and($2, 0x10) && and($2, 0x40)) || (and($2,0x80) && !and($2,0x10)))?"XS:A:+":"XS:A:-";

            if(libtype=="fr-secondstrand") XS1=((and($2, 0x10) && and($2, 0x80)) || (and($2,0x40) && !and($2,0x10)))?"XS:A:+":"XS:A:-";

        }

        else if($2~/^[:alpha:]/){   # FLAG in string

            if(libtype=="fr-firststrand") XS1=($2~/r.*1/ || ($2~/2/ && $2!~/r/))?"XS:A:+":"XS:A:-";

            if(libtype=="fr-secondstrand") XS1=($2~/r.*2/|| ($2~/1/ && $2!~/r/))?"XS:A:+":"XS:A:-";

        }

        print XS1;

        if(save_discrepancy_to_file!="" && XS1!=XS0) print >> save_discrepancy_to_file;

    }

}

awk redirecting output in AWK: override or not?

print items > output-file

This redirection prints the items into the output file named output-file. The file name output-file can be any expression. Its value is changed to a string and then used as a file name (see Expressions).When this type of redirection is used, the output-file is erased before the first output is written to it. Subsequent writes to the same output-file do not erase output-file, but append to it. (This is different from how you use redirections in shell scripts.) If output-file does not exist, it is created. For example, here is how an awk program can write a list of BBS names to one file named name-list, and a list of phone numbers to another file named phone-list:

          $ awk '{ print $2 > "phone-list"
          >        print $1 > "name-list" }' BBS-list
          $ cat phone-list
          -| 555-5553
          -| 555-3412
          ...
          $ cat name-list
          -| aardvark
          -| alpo-net
          ...

Each output file contains one name or number per line.

print items >> output-file

This redirection prints the items into the pre-existing output file named output-file. The difference between this and the single-‘>’ redirection is that the old contents (if any) of output-file are not erased. Instead, the awk output is appended to the file. If output-file does not exist, then it is created.

Here is the example code:

cat > test.txt

a

b

c

d

awk 'BEGIN{print "e1" >> "test.txt"; print "e2" >> "test.txt"; print "e3" >> "test.txt";}'

$ cat test.txt 

a

b

c

d

e1

e2

e3

awk 'BEGIN{print "e1" > "test.txt"; print "e2" > "test.txt"; print "e3" > "test.txt";}'

$ cat test.txt 

e1

e2

e3

But I just don't understand why the test.txt is gone in code below:

awk 'BEGIN{print "e1" > "test.txt"; system("rm test.txt");print "e2" > "test.txt"; print "e3" > "test.txt";}'

simple way to get reads length distribution of FASTQ files

• Using perl: 

cat input.fq | perl -ne '$s=<>;<>;<>;chomp($s);print length($s)."\n";' > input.readslength.txt

• Using awk:

cat input.fq | awk '{if(NR%4==2) print length($1)}' > input.readslength.txt

• if zipped file, using:

 zcat input.fq.gz | ...

• get length statistics:

sort input.readslength.txt | uniq -c

• or using fancy code from Jim Kent's utility:

textHistogram


So, one line code for all input fastq files would be:

find *.fq.gz -not -name \*raw\* -printf "zcat %p | awk '{if(NR%%4==2) print length(\$1)}' | textHistogram -maxBinCount=59 stdin \n" | sh


Note that you have to use double % to escape the % character for printf formatting control, just like in C. (Thanks for zencuke's answer here)

You will get something like this:

RNAseq.20E_library.result_primary.clean.fa

large values truncated: need 35 bins or larger binSize than 1

Maximum value 53.000000

 18 ********* 27730

 19 ******************* 58997

 20 ************************************************************ 186919

 21 ************************ 74536

 22 ************** 45171

 23 ************* 39107

 24 *************** 45560

 25 *********************** 70452

 26 ************************************************* 154030

 27 ************************************************************ 187704

 28 ************************** 81198

 29 ***** 17016

 30 ** 5341

 31 * 2439

 32  0

 33  1

 34  0

 35  0

 36  0

 37  0

 38  0

 39  1

<minVal or >= 40  173

RNAseq.Day_13_library.result_primary.clean.fa

large values truncated: need 35 bins or larger binSize than 1

Maximum value 53.000000

 18 ***** 22570

 19 ********* 40335

 20 ***************************** 127999

 21 ********************* 95179

 22 ****************** 79808

 23 ********* 39596

 24 ******* 29423

 25 ************* 55438

 26 ************************************* 164868

 27 ************************************************************ 265722

 28 *************************** 120353

 29 ********* 38625

 30 *** 14684

 31 * 5214

 32  0

 33  0

 34  1

<minVal or >= 35  140

split a file at specific line

I was facing a task: extracting the FASTA sequence from a GFF (or GFF3) file. Well, I just noticed there is a Bio::Perl script for that (http://www.bioperl.org/wiki/Getting_Fasta_sequences_from_a_GFF). But why not using a simple bash/awk for that? I like one-line coding :)

Here it is:

awk '/^>/,/############/' scaffold.gff

the grammer is awk '/from/, /to/' filename, to get lines from the one containing "from" to the line containing "to".  From my tests, I did not get very clear when there are multiple "from" and/or "to". So, be careful!

Another option is:

awk '/^>/ {p=1}; p==1 {print}' scaffold.gff

So, the final code for catenate all extracted fasta sequences is:

find -name scaffold\*.gff -exec sh -c "awk '/^>/ {p=1}; p==1 {print}' {} >> Manduca_gff_files_version_1.scaffold.fa" \; &

If you care the order of output, e.g. scaffold0001, scaffold0002, scaffold0003, ... then you have to sort the result of find (because find does not sort) first:

for i in `find -name scaffold\*.gff | sort`; do awk '/^>/ {p=1}; p==1 {print}' $i >> Manduca_gff_files_version_1.scaffold.fa; done

References:

1. http://www.unix.com/shell-programming-scripting/6959-split-file-specified-string.html
2. http://www.bioperl.org/wiki/Getting_Fasta_sequences_from_a_GFF

print all columns except last one

in awk:

print all columns except the first one column:
awk '{$1=""; print $0}' file

print all columns except the last one column:
awk '{$NF=""; print $0}' file

or acting as a geek, you can use the "rev" command to reverse the lines, then cut from the second field, then rev again:

rev file | cut -f2- | rev

Motivated from:http://lowfatlinux.com/linux-columns-cut.html#ixzz1nFIa9DZm

The unix command 'column' seems working for this purpose, but I did not figure it out yet. Here is column document:
http://www.eskimo.com/~scs/src/column.man.html

Ternary ifelse ( ?: ) in different languages

• AWK

$ awk 'ORS=NR%3?",":"\n"' student-marks

• Perl /PHP

$result = ($a > $b) ? $x : $y;

In Per6, use double ? and ! instead.

$result = ($a > $b) ?? $x !! $y;

• R

ifelse(a>0,a,0)

Ternary operator (if?true:false)

• bash/linux

ternary operator ? : is just short form of if/else

case "$b" in
 5) a=$c ;;
 *) a=$d ;;
esac

Or

 [[ $b = 5 ]] && a="$c" || a="$d"

Reference:

http://en.wikipedia.org/wiki/Ternary_operation

http://en.wikipedia.org/wiki/%3F:

awk.gsub, awk.system

Tips in Computer: "awk '{system('grep ' $4 ' file > result.txt'); getline result < 'result.txt'; close('result.txt'); print result;}' infile

e.g File
/scratch/dongx/jobid_2622600/18_2_62LPFAAXX_092810_MyersLab_4242.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam
/scratch/dongx/jobid_2622600/18_2_62LPFAAXX_092810_MyersLab_4242.fastq.tophat.mm_1.unique.hg19.out/junctions.bed
/scratch/dongx/jobid_2622604/20_2_62MA2AAXX_101210_MyersLab_4430_4.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam
/scratch/dongx/jobid_2622604/20_2_62MA2AAXX_101210_MyersLab_4430_4.fastq.tophat.mm_1.unique.hg19.out/junctions.bed
/scratch/dongx/jobid_2622606/20_5_62MA2AAXX_101210_PhiX.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam

You want to copy the .bam and .bed file to another folder...

$ cat > tophat.log
$ grep "s\.b" tophat.log | awk '{a=$1; gsub(".*[0-9]/","", a); sub("/",".", a); system("cp " $1 " " a);}'

- Sent using Google Toolbar"

PASS parameter to AWK - The UNIX and Linux Forums

PASS parameter to AWK - The UNIX and Linux Forums: "You can also assign variables on the command line, e.g....

awk -v a=1 -v b=2 -v c=3 -f test.awk

The -v option means that the variables are assigned before the BEGIN section.

- Sent using Google Toolbar"