Issue with seq-gc

[peebeenojay]

Hello,

I am trying to run your gggenomes tutorial and I am having some trouble for the GC content as shown in the tutorial (here).
I am using your emales.fna data and have run the protocol (with some hiccups but) successfully until this point.

I believe the issue I am having is related to samtools/bedtools, which are needed for this part of the analysis.

After installing both, I try to run:
seq-gc -Nbw 50 emales.fna > emales-gc.tsv

Actually like this, because I am running on macOS Monterey 12.2.1 Terminal:
perl seq-gc.pl -Nbw 50 emales.fna > emales-gc.tsv

And I get the error messages:

[E::fai_build3_core] Failed to open the file seq-gc
[faidx] Could not build fai index seq-gc.fai
cut: seq-gc.fai: No such file or directory
Error: The requested fasta database file (seq-gc) could not be opened. Exiting!

I have tried running with and without the 'perl' and the '.pl', and this is the only way it recognises there is even a seq-gc to run.

Full disclosure: although I am very comfortable in R, I am not used to using any of the other languages that are present in the tutorial. Nevertheless, I have been able to run everything up until this point, which is what made me think it may be an issue with the code itself.

Thank you for your time in advance, and thank you for gggenomes!
Cátia

[thackl]

Hi Catia,

seq-gc isn't a perl but a bash script. Try

bash seq-gc -Nbw 50 emales.fna > emales-gc.tsv

(if that works on your terminal)

And sorry for the state of the tutorial. It needs an overhaul, but I haven't had time yet.

[peebeenojay]

Hello again,

So I tried that, and still got the exact same error messages.

I ran:

bash seq-gc.sh -Nbw 50 emales.fna > emales-gc.tsv

With the .sh because it said no such file if I ran without it.
And got the messages:

[E::fai_build3_core] Failed to open the file seq-gc
[faidx] Could not build fai index seq-gc.fai
cut: seq-gc.fai: No such file or directory
Error: The requested fasta database file (seq-gc) could not be opened. Exiting!

Thank you again for your time,
Cátia

[thackl]

Hm, so I can't reproduce this on my system and I don't have a mac to try, so I'm pretty much guessing. The error is weird in the sense that the script is trying to open the actual scriptfile (seq-gc) instead of the fasta file (emales.fna) and I have no idea, why that would be. But what I also find strange is that you cannot just run ./seq-gc -h and need to do bash seq-gc.sh. This could be part of the problem. Can you try and run chmod +x seq-gc on the original file. That should make it executable if it isn't already. Then ./seq-gc -h should work.

If all of that fails, you just need the GC values to continue with the tutorial, you can read those from the bundled example data. Something like read_feats(ex("emales/emales-gc.bed")) in R.

[clb21565]

I was able to make it work by changing

"ARGS=getopt --name "seq-gc"
--options "w:s:bNmhV"
-- "$@"" -

which has mac OS specific syntax to

ARGS=$(getopt -o w:s:bNmhV --name "seq-gc" -- "$@")

new script ...

edit: whoops - also need to remove line breaks.

#!/bin/bash



## subs

usage(){

cat <<EOF

Usage:

  seq-gc FASTA > FASTA.gc.tsv



Calculate per sequence / window GC content.

  -w  Window size. Unless set, GC is calculated per sequence.

  -s  Shift between window starts. Default -w.

  -m  Calculate median of windows. Default -w=1000.

  -b  Output bed instead of tsv format for -w mode.

  -N  ignore Ns (and other non ATGC) in per sequence mode

      gives AT+GC = 100% even on scaffolds, 0 for N only.

EOF

exit 0;

}



check_bin(){

    hash $1 || { echo "$1 required in PATH" >&2; exit 1;}

}



## prep

[[ $# -eq 0 ]] && usage;



# Execute getopt
ARGS=$(getopt -o w:s:bNmhV --name "seq-gc" -- "$@")  

#Bad arguments

[ $? -ne 0 ] && exit 1;



# A little magic

eval set -- "$ARGS"



# Now go through all the options

while true; do

    case "$1" in

        -w)

            [ ! -n "$2" ] && (echo "$1: value required" 1>&2 && exit 1);

            WIN=$2;

            shift 2;;

        -s)

            [ ! -n "$2" ] && (echo "$1: value required" 1>&2 && exit 1);

            SHIFT=$2;

            shift 2;;

        -m)

            MEDIAN=1;

            shift;;

        -b)

            BED=1;

            shift;;

        -N)

            NIGNORE=1;

            shift;;

        -h)

	    usage && exit 0;;



        -V)

            grep -m1 'Version' "$0" | sed 's/.*Version\s*//';

            exit 0;;

        --)

            shift

            break;;

        *)

            echo "$1: Unknown option" 1>&2 && exit 1;;

    esac

done





## prep

[[ $# -eq 0 ]] && usage;



check_bin samtools;

check_bin bedtools;



FA=$1;

PRE=`basename $FA`

PRE=.${PRE%.*};

[[ -n $MEDIAN ]] && [[ -z $WIN ]] && WIN=1000;

[[ -z $SHIFT ]] && SHIFT=$WIN



## main

samtools faidx $FA;



if [[ -z $WIN ]]; then

    cut -f 1,2 $FA.fai | sed 's/\t/\t0\t/' > $PRE-gc.bed

else

    cut -f 1,2 $FA.fai > $PRE.len

    bedtools makewindows -g $PRE.len -w $WIN -s $SHIFT > $PRE-gc.bed

fi;



if [[ -z $MEDIAN ]]; then  # no median

    if [[ -z $NIGNORE ]]; then  # dont ignore Ns

        if [[ -z $BED ]]; then  # output bed

            bedtools nuc -fi $FA -bed $PRE-gc.bed |
                cut -f1,5 |
                tail -n +2

        else

            bedtools nuc-fi $FA -bed $PRE-gc.bed |

        gawk -v w=${WIN} 'BEGIN{FS="\t"; OFS="\t"}

{

if (FNR>1) {print $1,$2,$3,"GC.w"w,$5}

}'

        fi;

    else  # ignore Ns

        if [[ -z $BED ]]; then  # output bed

        bedtools nuc -fi $FA -bed $PRE-gc.bed |

            tail -n +2 |

            perl -ane '

              my $gcc = $F[6] + $F[7];

              my $tot = $gcc + $F[5] + $F[8];

              printf "$F[0]\t%0.6f\n", $tot ? $gcc/$tot : 0;

            '

        else

            bedtools nuc -fi $FA -bed $PRE-gc.bed |

                tail -n +2 |

                perl -ane '

                  my $gcc = $F[6] + $F[7];
                  my $tot = $gcc + $F[5] + $F[8];
                  printf "%s\t%s\t%s\tGC.w'$WIN'\t%0.6f\n", @F[0..2], $tot ? $gcc/$tot : 0;

                '

        fi;

    fi;

else

    bedtools nuc -fi $FA -bed $PRE-gc.bed |

    cut -f1,5 |

    tail -n +2 |

    perl -ane '

      if($F[0] eq $id){

          push @w, $F[1]

      }else{

          if(@w){

              @w=sort{$a<=>$b}@w;

              print $id,"\t",$w[@w/2],"\n"

          }

          $id=$F[0];
          @w=($F[1]);

      };

      END{

          if(@w){
              @w=sort{$a<=>$b}@w;
              print $id,"\t",$w[@w/2],"\n"
          }
      }'
fi;


rm -f $PRE-gc.bed $PRE.len;