This is a test version of Biostars. For the public version, visit https://www.biostars.org.
grep unexpected behaviour

I have a list of sequence names:

c32026_g2_i1
c43297_g1_i9
c45863_g2_i2
c43297_g1_i10
c35765_g2_i1
c44444_g3_i1
...

I want to take each sequence from the list and output the annotation of that sequence extracting it from the annotations file

This is what I do and it doesn't work well, for example in the list there is:

c38478_g1_i5
c38478_g1_i4
c38478_g1_i17
c38478_g1_i9
c38478_g1_i18
c38478_g1_i1

and it outputs:

c38478_g1_i5    1018    Q3B724    ...
c38478_g1_i4    1000    Q3B724    ...
c38478_g1_i17    887    Q3B724    ...
c38478_g1_i9    1007    Q3B724    ...
c38478_g1_i18    738    Q3B724    ...
c38478_g1_i1    496    -    -    -    - ...
c38478_g1_i10    950    -    -    -   ...
c38478_g1_i11    249    Q3B724    ...
c38478_g1_i12    706    -    -    -    ...
c38478_g1_i13    654    -    -    -    ...
c38478_g1_i14    809    -    -    -  ...
c38478_g1_i15    217    -    -    -  ...
c38478_g1_i16    788    Q3B724    ...
c38478_g1_i17    887    Q3B724    ...
c38478_g1_i18    738    Q3B724    ...
c38478_g1_i19    548    -    -    -    ...

for f in $(cat list.OE.txt); do grep $f Trinity_uniref_2015_02_filt_ann_out.txt; done > OE.annocript.txt
grep
• 3,582 views • link
updated by Ram Ram • written by Illinu Illinu

What are the expected input and output? If you are just greping the list from a file, and your list are store in a file, let's say, list.txt, then you can always do

grep -wf list.txt Trinity_uniref_2015_02_filt_ann_out.txt > OE.annocript.txt
• 0 views • link
updated by Ram Ram • written by Sam Sam

You need "-wFf". When list.txt is huge, "-F" will be much faster.

• 0 views • link
written by lh3 lh3

Yes, that worked :) Thank you

• 0 views • link
written by Illinu Illinu

3 answers

This is expected behavior since you didn't use the -w flag. man grep

For example,

grep c38478_g1_i1

Would return such lines:

c38478_g1_i11
c38478_g1_i1111111
c38478_g1_i122343454
as_long_as_c38478_g1_i1_is_on_the_line
• 0 views • link
updated by Ram Ram • written by 5heikki 5heikki

Great, I see! thank you

• 0 views • link
written by Illinu Illinu

Avoid grep in a loop. This is O(MxN) time complexity.

A more general solution: use hash to store the ids, then check if the data file has that id in a specific column.

awk 'BEGIN{FS=OFS="\t";while(getline<"ids.lst") ids[$1] =1}{if(ids[$1]==1) print}' full_anno.txt
• 0 views • link
updated by Ram Ram • written by tszn1984 tszn1984

A more general description for such case:

  • Data file: a matrix file, with ID in the Nth column
  • ID list file: a list of IDs
>selectIn datafile  IDfile N

selectIn source file

###################
#!/bin/sh
#Last-modified: 13 Apr 2013 03:02:05 PM
USAGE=" Usage: $0 Data.txt id.lst [col=1]"
case $# in
    0) echo $USAGE
       exit;;
    1) echo $USAGE
       exit;;
    2) col=1;;
    *) col=$3;;
esac
data=$1
awk -vLST=$2 -vCOL=$col 'BEGIN{FS=OFS="\t";while(getline<LST) count[$1]=1}{if(count[$(COL)]==1) print}' $data
###################
• 1 views • link
updated by Ram Ram • written by tszn1984 tszn1984

Just to update, this form:

for f in $(cat list.OE.txt); do grep -w $f Trinity_uniref_2015_02_filt_ann_out.txt; done > OE.annocript.txt

works A LOT faster than this one:

grep -wf list.OE.txt Trinity_uniref_2015_02_filt_ann_out.txt > OE.annocript.txt
• 0 views • link
written by Illinu Illinu

If you have large files it will be extremely slow anyway since you're searching the entire file again and again. For such tasks, there's man join

• 0 views • link
written by 5heikki 5heikki

Log in to answer this question.

More posts like this

No similar posts found.