Run program with different chromosome number every iteration
1
0
Entering edit mode
4.4 years ago
akoizumi • 0

I'm running the FUSION.assoc_test.R program from the fusion_twas package and I'm trying to write a simple piece of code that will allow me to loop through the program and substitute with a different chromosome number with every iteration.

So far I have:

for i in {1..22..1}
do
   Rscript ./FUSION.assoc_test.R \
    --sumstats PGC2.SCZ.sumstats \
    --weights ./WEIGHTS/GTEx.Whole_Blood.pos \
    --weights_dir ./WEIGHTS/ \
    --ref_ld_chr ./LDREF/1000G.EUR. \
    --chr $i \
    --out "PGZ_"$i".dat"
done

exit 0

When I try this I get an error:

Error in file(file, "rt") : cannot open the connection
Calls: read_plink -> read.table -> file
In addition: Warning message:
In file(file, "rt") :
  cannot open file './LDREF/1000G.EUR..bim': No such file or directory
Execution halted

which I assume arises because there wasn't a chromosome defined.

Any help would be appreciated.

Thanks!

twas software error • 1.6k views
ADD COMMENT
1
Entering edit mode

and without the loop ? does it work ?

 Rscript ./FUSION.assoc_test.R \
    --sumstats PGC2.SCZ.sumstats \
    --weights ./WEIGHTS/GTEx.Whole_Blood.pos \
    --weights_dir ./WEIGHTS/ \
    --ref_ld_chr ./LDREF/1000G.EUR. \
    --chr 1 \
    --out "PGZ_1.dat"
ADD REPLY
0
Entering edit mode

Thanks for the tip on readability- It works fine without the loop.

ADD REPLY
1
Entering edit mode

and what is the shell ? it's not bash because {1..2..1} is not a bash syntax.

using 'echo', what is the output of ?

for i in {1..22..1}
do
   echo Rscript ./FUSION.assoc_test.R \
    --sumstats PGC2.SCZ.sumstats \
    --weights ./WEIGHTS/GTEx.Whole_Blood.pos \
    --weights_dir ./WEIGHTS/ \
    --ref_ld_chr ./LDREF/1000G.EUR. \
    --chr $i \
    --out "PGZ_"$i".dat"
done
ADD REPLY
0
Entering edit mode

I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

101010 Button

ADD REPLY
0
Entering edit mode

I assume you checked if the file './LDREF/1000G.EUR..bim' exists?
The .. looks a bit strange, so I guess --ref_ld_chr ./LDREF/1000G.EUR. should be changed to --ref_ld_chr ./LDREF/1000G.EUR

ADD REPLY
0
Entering edit mode
4.4 years ago

try to use the bash syntax:

for i in $(seq 1 22 )
do
   Rscript ./FUSION.assoc_test.R \
    --sumstats PGC2.SCZ.sumstats \
    --weights ./WEIGHTS/GTEx.Whole_Blood.pos \
    --weights_dir ./WEIGHTS/ \
    --ref_ld_chr ./LDREF/1000G.EUR. \
    --chr $i \
    --out "PGZ_$i.dat"
done
ADD COMMENT

Login before adding your answer.

Traffic: 2156 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6