1

There are two input files; both are formatted as TSV.

File1: treated.bam.tsv 
File2: untreated.bam.tsv

Both files have the same fields, as listed below. (I numbered then for presentation purposes) – each file has 23 fields.

1chrom           9mismatches_pp       17C_pp
2pos             10deletions          18T
3ref             11deletions_pp       19T_pp
4reads_all       12insertions         20G
5reads_pp        13insertions_pp      21G_pp
6matches         14A                  22N
7matches_pp      15A_pp               23N_pp
8mismatches      16C

If values in the 1st and 2nd columns (chrom, pos) are the same in both files, I want to extract some of the fields in the records, then make a new output file like this below. The output file has 15 fields, combining data from the two input files, as below.

From file1:
  1chrom
  2pos
  3ref
  4reads_all
  8mismatches
  10deletions
  12insertions
  pct_file1 (the values from file1: (8mismatches+10deletions+12insertions)/4reads_all)
From file2:
  3ref
  4reads_all
  8mismatches
  10deletions
  12insertions
  pct_file2 (the values from file2: (8mismatches+10deletions+12insertions)/4reads_all)
-New values from extractions.
  pct_sub  (the values from pct_file1 - pct_file2: ((8mismatches+10deletions+12insertions)/4reads_all) - ((8mismatches+10deletions+12insertions)/4reads_all))

In the output file, the first eight columns come from File1, treated.bam.tsv (the 8th column is the value calculated with 8mismatches, 10deletions, 12insertions and 4reads_all from File1).

The rest of them come from File2, untreated.bam.tsv, and 13th column is also the value calculated with 8mismatches, 10deletions, 12insertions and 4reads_all from File2.

The last field, pct_sub, is calculated from the subtraction value with fields from File1 ((8mismatches+10deletions+12insertions)/4reads_all) and from File2 ((8mismatches+10deletions+12insertions)/4reads_all).

And how can I add new column names in output files, like pct_file1, pct_file2, pct_sub?

This is what I made for the above output file. (Input and output files have all same format: TSV.)

awk 'FNR==NR{array[$1,$2]=$0;next} { if ( $1 $2 in array ) print $1, $2, array[$3], array[$4], array[$8], array[$10], array[$12], (array[$8]+array[$10]+array[$12])/array[$4], $3, $4, $8, $10, $12, ($8+$10+$12)/$4, ((array[$8]+array[$10]+array[$12])/array[$4])-(($8+$10+$12)/$4) > "awkoutput.bam.tsv" }' treated.bam.tsv untreated.bam.tsv

(Actually, $1, $2 are not a problem from File1 or File2)

FILE1 (treated)

chrom pos ref reads_all reads_pp matches matches_pp mismatches mismatches_pp deletions deletions_pp insertions insertions_pp A A_pp C C_pp T T_pp G G_pp N N_pp

chrY 59363551 G 8 0 7 0 0 0 1 0 5 0 0 0 0 0 0 0 7 0 0 0
chrY 59363552 G 7 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0
chrY 59363553 T 7 0 7 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0
chrY 59363554 G 7 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0
chrY 59363555 T 7 0 7 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0

FILE2 (untreated)

chrom pos ref reads_all reads_pp matches matches_pp mismatches mismatches_pp deletions deletions_pp insertions insertions_pp A A_pp C C_pp T T_pp G G_pp N N_pp
chrY 59363551 G 2 0 2 0 0 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0
chrY 59363552 G 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
chrY 59363553 T 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
chrY 59363554 G 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
chrY 59363555 T 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0

OUTPUT

chrom pos ref reads_all mismatches deletions insertions pct_file1 ref reads_all mismatches deletions insertions pct_file2 pct_sub
chrY 59363551 G 8 0 1 5 0.75 G 2 0 0 1 0.5 0.25
chrY 59363552 G 7 0 0 0 0 G 1 0 0 0 0 0
chrY 59363553 T 7 0 0 0 0 T 1 0 0 0 0 0
chrY 59363554 G 7 0 0 0 0 G 1 0 0 0 0 0
chrY 59363555 T 7 0 0 0 0 T 1 0 0 0 0 0
0

2 Answers 2

1

Bioinformatics seems to be fun. If you are open to non-awk solutions, an easy job for miller:

mlr --itsv join -u -j chrom,pos --lp tr_ --rp untr_ -f treated.bam.tsv untreated.bam.tsv | # join data from treated and untreated files by fields chrom and pos
mlr put '$tr_pct=($tr_mismatches+$tr_deletions+$tr_insertions)/$tr_reads_all' | # calculate pct for treated data
mlr put '$untr_pct=($untr_mismatches+$untr_deletions+$untr_insertions)/$untr_reads_all' | # calculate pct for untreated data
mlr cut -o -f chrom,pos,tr_ref,tr_reads_all,tr_mismatches,tr_deletions,tr_insertions,tr_pct,untr_ref,untr_reads_all,untr_mismatches,untr_deletions,untr_insertions,untr_pct | # remove superfluous fields
mlr --otsv put '$pct_sub=$tr_pct-$untr_pct' # append pct subtraction field

chrom   pos tr_ref  tr_reads_all    tr_mismatches   tr_deletions    tr_insertions   tr_pct  untr_ref    untr_reads_all  untr_mismatches untr_deletions  untr_insertions untr_pct    pct_sub
chrY    59363551    G   8   0   1   5   0.750000    G   2   0   0   1   0.500000    0.250000
chrY    59363552    G   7   0   0   0   0   G   1   0   0   0   0   0
chrY    59363553    T   7   0   0   0   0   T   1   0   0   0   0   0
chrY    59363554    G   7   0   0   0   0   G   1   0   0   0   0   0
chrY    59363555    T   7   0   0   0   0   T   1   0   0   0   0   0

It looks scarier than it is. Really.

0
  1. if ( $1 $2 in array ) doesn’t work; you have to do if (($1,$2) in array).
  2. You can’t use array[$3] and array[$4] like that.  Your array looks like
    array[chrY,59363551]=“chrY 59363551 G 8 0 7 0 0 0 1 0 5 0 0 0 0 0 0 0 7 0 0 0”
    array[chrY,59363552]=“chrY 59363552 G 7 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0”
                 ︙
    and when you say array[$3] and array[$4], you’re referring to array[G] and array[2], etc., which do not exist.
  3. The ability to specify > "filename" in awk code is a great feature when you want to write to multiple files.  It’s not so useful when you have only one output file — why don’t you just redirect the output from the awk command?
  4. Long lines are bad.  Break long commands into short lines.  Reduce duplication by reusing variables.
  5. Don’t use an array called array.  That’s like having a variable called variable, a file called file, a person called Person, etc.  Use descriptive names.

OK, that said,

awk 'FNR==NR {file1data[$1,$2]=$0; next}
        {       if (($1,$2) in file1data) {
                        # Save desired values from file2.
                        file2arg03=$3
                        file2arg04=$4
                        file2arg08=$8
                        file2arg10=$10
                        file2arg12=$12
                        pct_file2=($8+$10+$12)/$4
                        # Get data from file1.
                        $0=file1data[$1,$2]
                        pct_file1=($8+$10+$12)/$4
                        print $1, $2, $3, $4, $8, $10, $12, pct_file1, \
                                file2arg03, file2arg04, file2arg08, file2arg10, file2arg12, \
                                pct_file2, pct_file1-pct_file2
                } else printf "(%s,%s) in file2 but not file1.%s", $1, $2, ORS
        }' treated.bam.tsv untreated.bam.tsv > awkoutput.bam.tsv

Like your version, this saves the file1 data in an array and then does all the work/output while it’s reading file2.  Upon receiving a line from file2, it saves the desired fields from that line into named variables (we could also use another array, five elements long) and then it restores the data from the corresponding line in file1.  By assigning the entire line to $0, it causes $1, $2, $3, $4, etc., to be restored to their original values.

Are you really having a problem writing a header line in the output?  Try:

        {       if (FNR == 1) {
                        print "chrom pos ref reads_all mismatches deletions insertions pct_file1 …"
                } else if (($1,$2) in file1data ) {
                        file2arg03=$3
                              ︙

OK, here’s a version that’s closer in spirit to your attempt, and that handles the header line:

awk 'FNR==NR {file1line[$1,$2]=$0; next}
        {       if (FNR == 1) {
                        print "chrom pos ref reads_all mismatches deletions insertions pct_file1 ref reads_all mismatches deletions insertions pct_file2 pct_sub …"
                } else if (($1,$2) in file1line ) {
                        # Get data from file1.
                        split(file1line[$1,$2], file1arg)
                        pct_file1=(file1arg[8]+file1arg[10]+file1arg[12])/file1arg[4]
                        pct_file2=($8+$10+$12)/$4
                        print $1, $2, file1arg[3], file1arg[4], file1arg[8], \
                                file1arg[10], file1arg[12], pct_file1, \
                                $3, $4, $8, $10, $12, pct_file2, pct_file1-pct_file2
                } else printf "(%s,%s) in file2 but not file1.%s", $1, $2, ORS
        }' treated.bam.tsv untreated.bam.tsv > awkoutput.bam.tsv

This retrieves the line from file1 (from file1line) and passes it to split to break it down into its constituent 23 values, which get stored in the array file1arg.  Then it is able to use file1arg[3], file1arg[4], …, the way you were using array[$3], array[$4], …

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .