1

The input file that needs to be edited is as below (can have more rows):

bundle_id   target_id   length  eff_length  tot_counts  uniq_counts est_counts  eff_counts  ambig_distr_alpha   ambig_distr_beta    fpkm    fpkm_conf_low   fpkm_conf_high  solvable    tpm
1   intron_FBgn0035847:4_FBgn0035847:3  61  0   0   0   0   0   0.00E+00    0.00E+00    0.00E+00    0.00E+00    0.00E+00    F   0.00E+00
2   intron_FBgn0032515:2_FBgn0032515:4  72  0   0   0   0   0   0.00E+00    0.00E+00    0.00E+00    0.00E+00    0.00E+00    F   0.00E+00
3   intron_FBgn0266486:5_FBgn0266486:4  58  0   0   0   0   0   0.00E+00    0.00E+00    0.00E+00    0.00E+00    0.00E+00    F   0.00E+00
4   intron_FBgn0031359:10_FBgn0031359:7 4978    1430.739479 91  0   30.333333   105.539363  1.00E+00    1.00E+00    6.30E+00    1.77E+00    1.08E+01    F   1.42E+01
4   intron_FBgn0031359:10_FBgn0031359:8 4978    1430.739479 91  0   30.333333   105.539363  1.00E+00    1.00E+00    6.30E+00    1.77E+00    1.08E+01    F   1.42E+01
4   intron_FBgn0031359:10_FBgn0031359:9 4978    1430.739479 91  0   30.333333   105.539363  1.00E+00    1.00E+00    6.30E+00    1.77E+00    1.08E+01    F   1.42E+01
536 intron_CR31143:1_CR31143:2  40  0   0   0   0   0   0.00E+00    0.00E+00    0.00E+00    0.00E+00    0.00E+00    F   0.00E+00

For each ID in the 2nd column intron_XXXXXXXX:X_XXXXXXXX:X, I want to extract the string between intron_ and the 1st : (in between the string usually but not always starts with FBgn).

Then I have a list as following (one column for FBgn and the other column for corresponding name I want the FBgn to be converted into):

## FlyBase Gene Mapping Table   
## Generated: Fri Dec 20 12:37:29 2013  
## Using datasource: dbi:Pg:dbname=fb_2014_01_reporting;host=flysql9;port=5432...   
FBgn0035847 mthl7
FBgn0032515 loqs
FBgn0266486 CG45085
FBgn0031359 CG18317

Then I want to search the extracted string in the list's 1st column.

If the extracted string has corresponding value in the 2nd column, I want to replace the whole ID intron_FBgnXXXXXX:X_FBgnXXXXXX:X with the corresponding name in the 2nd column.

If the extracted string does not exist in the 1st column, I want to replace the whole ID intron_XXXXXXXX:X_XXXXXXXX:X with the extracted string.

I have a script as the following:

ref="gene_map_table_fb_2014_01_short.tsv"
target="HC25_LNv_ZT02_intron_results.txt"
output="temptemp.txt"

declare -A map
while read line
do
if [[ ! -z "$line" ]] && [[ ! "$line" =~ ^#.* ]]
then
key=$(echo "$line" | cut -f 1)
value=$(echo "$line" | cut -f 2)
map[$key]=$value
fi
done < $ref

while read line
do
  key=$(echo "$line" | sed -n 's/.*_\([^\.]*\)\:.*/\1/p' | head -1)
if [ ! -z "$key" ]
then
  echo "$line" | sed 's/intron_[^[:space:]]*/'${map[$key]}'/g' >> $output
else
  echo "$line" | sed 's/intron_[^[:space:]]*/'$key'/g' >> $output
fi
done < $target

Everything seems to work fine except that the output file lack the lines whose ID does not start with FBgn.

  • It would benefit you, if you provide a simpler but concrete example that covers all the possible scenarios. For instance, you could remove most of the unwanted information in your input file so that the question will look lot more appealing. – Ramesh Oct 9 '14 at 21:16
0

You can do it:

cat gene_map_table_fb_2014_01_short.tsv  |sed '1d' |awk {'print $2'} |awk 'BEGIN{FS=":"} {print $2}' |sed s/._//g

At first cat your file, Then delete first line (header of columns with d1) , then print all of column, then seprate 4_FBgn0035847 with awk 'BEGIN{FS=":"} {print $2}' Then eliminate number_ with sed s/._//g

Output is:

FBgn0035847
FBgn0032515
FBgn0266486
1FBgn0031359
1FBgn0031359
1FBgn0031359
CR31143

However if your end line is extra and if you want to remove it , you can do it:

cat gene_map_table_fb_2014_01_short.tsv  |sed '1d' |awk {'print $2'} |awk 'BEGIN{FS=":"} {print $2}' |sed s/._//g |sed '$d'

So, output is :

FBgn0035847
FBgn0032515
FBgn0266486
1FBgn0031359
1FBgn0031359
1FBgn0031359
0

Using awk

This creates tab-separated output:

$ awk -v OFS="\t" 'NR==FNR{a[$1]=$2;next} FNR==1{print;next} {sub(/intron_/, "", $2); sub(/:.*/,"",$2);if ($2 in a) $2=a[$2];print}' gene_map_table_fb_2014_01_short.tsv HC25_LNv_ZT02_intron_results.txt
bundle_id   target_id   length  eff_length  tot_counts  uniq_counts est_counts  eff_counts  ambig_distr_alpha   ambig_distr_beta    fpkm    fpkm_conf_low   fpkm_conf_high  solvable    tpm
1       mthl7   61      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
2       loqs    72      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
3       CG45085 58      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
536     CR31143 40      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00

Explanation:

  • -v OFS="\t"

    This makes the output field separator into a tab.

  • NR==FNR{a[$1]=$2;next}

    This creates an associative array, a, based on the first file on the command line, with the the first column as key and the second column as value. The next command instructs awk to skip the rest of the commands and jump to the next line.

    The mapping file contains some comment lines. We could have easily added an extra if statement to keep them from getting added into the array a. However, since they do no harm, we skipped that complication.

  • FNR==1{print;next}

    This prints the header line unchanged.

  • {sub(/intron_/, "", $2); sub(/:.*/,"",$2)

    This removes the fluff from the second field leaving behind only the string that we want.

  • `if ($2 in a) $2=a[$2]

    If the string from the second field is present as a key in array a, then we substitute its corresponding value.

  • print

    The revised line is printed.

Using bash

In the script, replace

if [ ! -z "$key" ]

With:

if [[ "$key" && "${map[$key]}" ]]

What the script appears to need to know at this point is whether key was present in map or not. The revised test makes sure not merely that key is non-empty but that it is in map.

With that one change, I get the output:

$ cat temptemp.txt 
bundle_id   target_id   length  eff_length  tot_counts  uniq_counts est_counts  eff_counts  ambig_distr_alpha   ambig_distr_beta    fpkm    fpkm_conf_low   fpkm_conf_high  solvable    tpm
1       mthl7   61      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
2       loqs    72      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
3       CG45085 58      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
4       CG18317 4978    1430.739479     91      0       30.333333       105.539363      1.00E+00        1.00E+00        6.30E+00        1.77E+00        1.08E+01        F       1.42E+01
536     CR31143 40      0       0       0       0       0       0.00E+00        0.00E+00        0.00E+00        0.00E+00        0.00E+00        F       0.00E+00

As an aside on text, [ ! -z "$key" ] returns true if key is non-empty. This is equivalent to [ -n "$key" ]. Since this is such a common test, it can be further shortened to [ "$key" ]. This could be used to simplify several lines in the bash script.

  • sorry, mapping table is my map list? and where should I put input file? – Karli Oct 9 '14 at 21:40
  • Sorry for the confusion. I just updated the answer with the names from your shell script. I also added a line which, I think, fixes your bash script. – John1024 Oct 9 '14 at 21:43
  • Thanks!the awk worked. I will test will full size file and maptable now. for the bash script I changed to if [[ "$key" && "${map[$key]}" ]], but result looks the same. I will check if it is due to something else or not. – Karli Oct 9 '14 at 21:55
  • I combined this part with the other part I asked earlier into one script and tested the full size file after making small changes. Everything is working now. awk is much faster than read the table into the map. I just need to working on naming the output files a little more. Thank you so much for all the help and explanations. Those are important for me so I can learn how to use later. – Karli Oct 9 '14 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.