1

I am trying to run the following bash script:

#!/bin/bash

file=$1
filename=${file%%.*}
line1=$(sed -n 1~2p ${file})
seqs=$(grep -v '^>' ${file})
pos=$(echo "${line1}" | awk -F"[__]" 'NF>2{print $2}')

( 
    awk -v str="${seqs}" -v str2="${pos}" -v str3="${line1}" -v name=${filename} -v sep="[$IFS]" '
        BEGIN {
            n = split(str, a, sep)
            m = split(str2, b, sep)
            k = split(str3, c, sep)
            for (i=1;i<=n;i++) {o=10;d[$i]=b[i]-o;s[$i]=d[i]>0?d[i]:1; print c[i] "\n" substr(a[i],d[$i],2*o+(d[$i]<0?d[$i]:1)) > name"_flanks.fasta"}
        }
    '
)

however I am getting:

$ ./test.sh myfile.fasta
./test.sh: line 10: /usr/bin/xargs: Argument list too long

Slap on the wrist for not using version control, but this was working in a previous version of my code. what seems to be the problem?

edit: have noticed that if I pipe "head ${file} |" into the sed and grep commands then this runs fine, but doing "cat ${file} |" re-produces the original error. could it really be a file size limitation? do I have to split the computation up into smaller file chunks?

output of "$seqs$ is around 6,000 of these elements

MEDEAVLDRGASFLKHVCDEEEVEGHHTIYIGVHVPKSYRRRRRHKRKTGHKEKKEKERISENYSDKSDIENADESSSSILKPLISPAAERIRFILGEEDDSPAPPQLFTELDELLAVDGQEMEWKETARWIKFEEKVEQGGERWSKPHVATLSLHSLFELRTCMEKGSIMLDREASSLPQLVEMIVDHQIETGLLKPELKDKVTYTLLRKHRHQTKKSNLRSLADIGKTVSSASRMFTNPDNGSPAMTHRNLTSSSLNDISDKPEKDQLKNKFMKKLPRDAEASNVLVGEVDFLDTPFIAFVRLQQAVMLGALTEVPVPTRFLFILLGPKGKAKSYHEIGRAIATLMSDEVFHDIAYKAKDRHDLIAGIDEFLDEVIVLPPGEWDPAIRIEPPKSLPSSDKRKNMYSGGENVQMNGDTPHDGGHGGGGHGDCEELQRTGRFCGGLIKDIKRKAPFFASDFYDALNIQALSAILFIYLATVTNAITFGGLLGDATDNMQGVLESFLGTAVSGAIFCLFAGQPLTILSSTGPVLVFERLLFNFSKDNNFDYLEFRLWIGLWSAFLCLILVATDASFLVQYFTRFTEEGFSSLISFIFIYDAFKKMIKLADYYPINSNFKVGYNTLFSCTCVPPDPANISISNDTTLAPEYLPTMSSTDMYHNTTFDWAFLSKKECSKYGGNLVGNNCNFVPDITLMSFILFLGTYTSSMALKKFKTSPYFPTTARKLISDFAIILSILIFCVIDALVGVDTPKLIVPSEFKPTSPNRGWFVPPFGENPWWVCLAAAIPALLVTILIFMDQQITAVIVNRKEHKLKKGAGYHLDLFWVAILMVICSLMALPWYVAATVISIAHIDSLKMETETSAPGEQPKFLGVREQRVTGTLVFILTGLSVFMAPILKFIPMPVLYGVFLYMGVASLNGVQFMDRLKLLLMPLKHQPDFIYLRHVPLRRVHLFTFLQVLCLALLWILKSTVAAIIFPVMILALVAVRKGMDYLFSQHDLSFLDDVIPEKDKKKKEDEKKKKKKKGSLDSDNDDSDCPYSEKVPSIKIPMDIMEQQPFLSDSKPSDRERSPTFLERHTSC

The file contains many repeats of data like:

>Q9UM01_334_L_R
MVDSTEYEVASQPEVETSPLGDGASPGPEQVKLKKEISLLNGVCLIVGNMIGSGIFVSPKGVLIYSASFGLSLVIWAVGGLFSVFGALCYAELGTTIKKSGASYAYILEAFGGFLAFIRLWTSLLIIEPTSQAIIAITFANYMVQPLFPSCFAPYAASRLLAAACICLLTFINCAYVKWGTLVQDIFTYAKVLALIAVIVAGIVRLGQGASTHFENSFEGSSFAVGDIALALYSALFSYSGWDTLNYVTEEIKNPERNLPLSIGISMPIVTIIYILTNVAYYTVLDMRDILASDAVAVTFADQIFGIFNWIIPLSVALSCFGGLNASIVAASRLFFVGSREGHLPDAICMIHVERFTPVPSLLFNGIMALIYLCVEDIFQLINYYSFSYWFFVGLSIVGQLYLRWKEPDRPRPLKLSVFFPIVFCLCTIFLVAVPLYSDTINSLIGIAIALSGLPFYFLIIRVPEHKRPLYLRRIVGSATRYLQVLCMSVAAEMDLEDGGEMPKQRDPKSN

I want to read in the header (beginning with ">"), strip out the position number (334) and then with line 2 being the "sequence" I want to:

go to position pos[i] in seqs[i] and choose a substring of seqs[i] that is up to 10 positions either side of pos[i]. for example if pos[i] = 15 I would return:

EYEVASQPEVETSPLGDGAS

I can do this when not using the entire file, however it appears that reading everything directly into awk would make the program more efficient than loading everything via shell variables.

  • The contents of seqs could be too large to pass on the command line (but could probably be passed to awk using a pipe). What is the output of grep -v '^>' ${file}|wc -c ? – Mark Plotnick Mar 10 '15 at 12:04
  • 1
    you can pass the filename to awk, and read if from within awk to value awk's vars, instead of passing the values via shell variables. less size limitations that way. – Olivier Dulac Mar 10 '15 at 12:26
  • @Mark Plotnick - have included output of seqs – brucezepplin Mar 10 '15 at 12:49
  • @Olivier Dulac - but would it still be possible to process using the sed and grep commands after reading the variables directly into awk? – brucezepplin Mar 10 '15 at 12:58
  • 2
    @brucezepplin the idea would be to do it all from within awk (which can neatly emulate sed and grep). the only caveat would be if you need the external (= shell level) variables later, in which case you should continue as you are trying to do. where is the xargs, by the way? – Olivier Dulac Mar 10 '15 at 13:00
1

Why you do not do it just with awk only as offered by @Olivier Dulac:

awk '/^>/{split($0,N,"_");n=N[2];print;next}{print substr($0,n-10,20)}' file > file_flanks.fasta

same:

awk -F'_' '/^>/{n=$2;print;next}{print substr($0,n-10,20)}' file > file_flanks.fasta

Or without array:

awk '/^>/{print;sub("[^_]*_","");n=$0+0;next}{print substr($0,n-10,20)}' file > file_flanks.fasta
  • have been trying for the last hour! and you have just done it for me! – brucezepplin Mar 10 '15 at 14:57

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.