2
|
1 ## This is a sample PBS script for profiling STR from reference genome using STR-FM
|
|
2 ##
|
|
3 ##requirement
|
|
4 ##1 STR error rates (can be downloaded from https://usegalaxy.org/u/guru%40psu.edu/h/error-rates-files) --> errorrate.bymajorallele
|
|
5 ##
|
|
6 echo " "
|
|
7 echo " "
|
|
8 echo "Job started on `hostname` at `date`"
|
|
9 cd /working/directory/
|
|
10 echo ${MOTIF}
|
|
11 echo ${OUTPUT}
|
|
12 echo " "
|
|
13 echo "Generate all possible combination of STR length profile" ## See detail in profilegenerator.xml on https://github.com/Arkarachai/STR-FM
|
|
14 python profilegenerator.py errorrate.bymajorallele ${MOTIF} 30 > ${OUTPUT}.30
|
|
15
|
|
16 echo "remove duplicated profiles"
|
|
17 cat ${OUTPUT}.30 | sort | uniq > ${OUTPUT}.30.sort
|
|
18
|
|
19 echo "genotyping using error correction model" ## See detail in GenotypingSTR.xml on https://github.com/Arkarachai/STR-FM
|
|
20 python GenotypeTRcorrection.py ${OUTPUT}.30.sort errorrate.bymajorallele ${OUTPUT}.30.prob 0.5
|
|
21
|
|
22 echo "select only full motif different --> need to replace 4 with motif size (1-6)"
|
|
23 cat ${OUTPUT}.30.prob | grep hetero | awk '(($7-$8)==4) || (($8-$7)==4) {print $0}' > ${OUTPUT}.30.prob.screen
|
|
24
|
|
25 echo "Evaluate the probability of the allele combination to generate read profile" ## See detail in probvalueforhetero.xml on https://github.com/Arkarachai/STR-FM
|
|
26 python heteroprob.py ${OUTPUT}.30.prob.screen ${INPUT} > ${OUTPUT}.30.bino
|
|
27
|
|
28 echo "formatting"
|
|
29 cat ${OUTPUT}.30.bino | sort -k 12n,12 -k 6n,6 > ${OUTPUT}.30.bino.sort
|
|
30
|
|
31 echo "Combine read profile probabilities" ## See detail in combineprobforallelecombination.xml on https://github.com/Arkarachai/STR-FM
|
|
32 python combinedprobforallelecombination.py ${OUTPUT}.30.bino.sort > ${OUTPUT}.30.bino.sort.plot
|
|
33
|
|
34
|
|
35 echo "Job end on `hostname` at `date`"
|