Compare some specific columns of lines within a file using bash script

bobia9193 0
Viewed 37

I want to compare the 2nd and 4th columns of lines in a file. In detail, line 1 with 2,3,4...N, then line 2 with 3,4,5...N, and so on.

I have written a script, it worked but running so long, over 30 minutes.

Let the number of lines is 1733 with header, my code is:

for line1 in {2..1733}; do \
    for line2 in {$((line1+1))..1733}; do \
        i7_diff=$(cmp -bl \
            <(sed -n "${line1}p" Lab_primers.tsv | cut -f 2) \
            <(sed -n "${line2}p" Lab_primers.tsv | cut -f 2) | wc -l);
        i5_diff=$(cmp -bl \
            <(sed -n "${line1}p" Lab_primers.tsv | cut -f 4) \
            <(sed -n "${line2}p" Lab_primers.tsv | cut -f 4) | wc -l);
        if [ $i7_diff -lt 3 ]; then
            if [ $i5_diff -lt 3 ]; then
                echo $(sed -n "${line1}p" Lab_primers.tsv)"\n" >> primer_collision.txt
                echo $(sed -n "${line2}p" Lab_primers.tsv)"\n\n" >> primer_collision.txt
            fi;
        fi;
    done
done

I used nested for loops then using sed to print exactly the $line, next using cut to extract the desired column. Finally, the cmp and wc command to count the number of differences of two columns of a pair lines.

If meeting the condition (both 2nd and 4th columns of pair of lines have the number of differences less than 3), the code will print a pair lines to output file.

Here is an excerpt of the input (it has 1733 lines):

I7_Index_ID  index     I5_Index_ID  index2    primer
D703         CGCTCATT  D507         ACGTCCTG  27
D704         GAGATTCC  D507         ACGTCCTG  28
D701         ATTACTCG  D508         GTCAGTAC  29
S6779        CGCTAATC  S6579        ACGTCATA  559
D708         TAATGCGC  D503         AGGATAGG  44
D705         ATTCAGAA  D504         TCAGAGCC  45
D706         GAATTCGT  D504         TCAGAGCC  46
i796         ATATGCGC  i585         AGGATAGC  R100
D714         TGCTTGCT  D510         AACCTCTC  102
D715         GGTGATGA  D510         AACCTCTC  103
D716         AACCTACG  D510         AACCTCTC  104
i787         TGCTTCCA  i593         ATCGTCTC  R35

Then the expected output is:

D703 CGCTCATT D507 ACGTCCTG 27
S6779 CGCTAATC S6579 ACGTCATA 559


D708 TAATGCGC D503 AGGATAGG 44
i796 ATATGCGC i585 AGGATAGC R100


D714 TGCTTGCT D510 AACCTCTC 102
i787 TGCTTCCA i593 ATCGTCTC R35

My question is what the better code to deal with it, how to reduce the running time?

Thank you for your help!

shell bash
1 Answers

You could start to sort by fields 2 and 4. Then, no need for double loop: if a pair exist, they should be adjacent.

sort -k 2,2 -k 4,4 myfile.txt

Then, we need to print only bunch of consecutive lines that share the same 2 and 4 fields.

first=yes
sort -k 2,2 -k 4,4 test.txt | while read l 
do 
   fields=(${l})
   new2=${fields[1]}
   new4=${fields[3]} # Fields 2 and 4, bash-way
   if [[ "$new2" = "$old2" ]] && [[ "$new4" = "$old4" ]]
   then 
      if [[ $first ]]
      then
         # first time we print something for this series: we need 
         # to also print the previous line (the first of the series)
         echo; echo "$oldl"
         # But if the next line is identical (series of 3, no need to repeat this line)
         first=
      fi
      echo "$l"
   else 
      # This line is not identical to the previous. So nothing to print
      # If the next one is identical to this one, then, this one will
      # be the first of its series
      first=yes
   fi
   old2=$new2
   old4=$new4
   oldl="${l}"
done

One frustrating thing: uniq -D almost does all the job we did in this script. Except that it is unable to filter on specific lines.

But we could also rewrites lines so that uniq can work. Not very fluent in awk (if I were, I am pretty sure awk could do the uniq work for me), but well

sort -k 2,2 -k 4,4 test.txt  | awk '{print $0" "$2" "$4}' | uniq -D -f 5 | awk '{printf "%-12s %-9s %-12s %-9s %s\n",$1,$2,$3,$4,$5}'

does the job.

sort sort the lines by fields 2 and 4. awk add to the end of each line a copy of fields 2 and 4. Which then make uniq usable, since uniq is able to ignore N first fields. So here, we use uniq ignoring the 5 1st fields, that is working only on the copies of fields 2 and 4. With -D uniq display only duplicate lines. Then, the last awk remove the copies of field we don't need anymore.

chrslg0
Related