|
0
|
1 #!/bin/bash
|
|
|
2
|
|
|
3 #TV-vs-background.sh $variants $genomes ${reference.fields.crr_path} ${reference.fields.31G_var_paths} ${reference.54G_var_paths} $threshold $output_all $output_filtered
|
|
|
4
|
|
|
5 echo $@
|
|
|
6
|
|
|
7 set -- `getopt -n$0 -u -a --longoptions="variants: reference: VN_varfiles: outputfile_filtered: outputfile_all: threshold: thresholdhc:" "h:" "$@"` || usage
|
|
|
8 [ $# -eq 0 ] && usage
|
|
|
9
|
|
|
10 while [ $# -gt 0 ]
|
|
|
11 do
|
|
|
12 case "$1" in
|
|
|
13 --variants) variants=$2;shift;;
|
|
|
14 --reference) crr=$2;shift;;
|
|
|
15 --VN_varfiles) VN_varfiles_list=$2;shift;;
|
|
|
16 --outputfile_filtered) output_filtered=$2;shift;;
|
|
|
17 --outputfile_all) output_all=$2;shift;;
|
|
|
18 --threshold) threshold=$2;shift;;
|
|
|
19 --thresholdhc) thresholdhc=$2;shift;;
|
|
|
20 -h) shift;;
|
|
|
21 --) shift;break;;
|
|
|
22 -*) usage;;
|
|
|
23 *) break;;
|
|
|
24 esac
|
|
|
25 shift
|
|
|
26 done
|
|
|
27
|
|
|
28 # replace newline chars with spaces for input to testvariants
|
|
2
|
29 echo "varfiles list: $VN_varfiles_list"
|
|
|
30
|
|
0
|
31 tr '\n' ' ' < $VN_varfiles_list > VN_varfiles.txt
|
|
|
32
|
|
|
33
|
|
|
34 ### run TestVariants against 31G, 54G or 85G
|
|
|
35
|
|
|
36 echo "number of normals: $VNsetsize"
|
|
|
37 echo "list of normals: ($VN_varfiles_list)"
|
|
|
38 cat VN_varfiles.txt
|
|
|
39
|
|
|
40
|
|
|
41 echo "running TV against Virtual Normal set"
|
|
|
42 echo "command: cgatools testvariants\
|
|
|
43 --beta \
|
|
|
44 --reference $crr \
|
|
|
45 --input $variants \
|
|
|
46 --output $output_all \
|
|
|
47 --variants `cat VN_varfiles.txt`"
|
|
|
48
|
|
|
49 cgatools testvariants \
|
|
|
50 --beta \
|
|
|
51 --reference $crr \
|
|
|
52 --input $variants \
|
|
|
53 --output $output_all \
|
|
|
54 --variants `cat VN_varfiles.txt`
|
|
|
55
|
|
|
56
|
|
|
57
|
|
|
58 VNsetsize=`cat $VN_varfiles_list | wc -l`
|
|
|
59
|
|
|
60
|
|
|
61
|
|
|
62 ### filter file based on occurrence in background genomes
|
|
|
63 cp $output_all $output_filtered
|
|
|
64 cp $output_all output_expanded
|
|
|
65
|
|
|
66 ### condens file to columns with counts for all background genomes
|
|
|
67 echo "Counting..."
|
|
|
68 awk 'BEGIN{
|
|
|
69 FS="\t";
|
|
|
70 OFS="\t";
|
|
|
71 totalnormals="'"$VNsetsize"'"+0
|
|
|
72 count["00"]="0";
|
|
|
73 count["01"]="0";
|
|
|
74 count["11"]="0";
|
|
|
75 count["0N"]="0";
|
|
|
76 count["1N"]="0";
|
|
|
77 count["NN"]="0";
|
|
|
78 count["0"]="0";
|
|
|
79 count["1"]="0";
|
|
|
80 count["N"]="0";
|
|
|
81 }{
|
|
|
82 if(FNR==1) # header
|
|
|
83 print $1,$2,$3,$4,$5,$6,$7,$8,"VN_occurrences","VN_frequency","VN_fullycalled_count","VN_fullycalled_frequency","VN_00","VN_01","VN_11","VN_0N","VN_1N","VN_NN","VN_0","VN_1","VN_N"
|
|
|
84 else{
|
|
|
85 #count entries in reference genomes
|
|
|
86 for (c in count)
|
|
|
87 count[c]=0;
|
|
|
88 for (i=9; i<=NF; i++){
|
|
|
89 count[$i]++;
|
|
|
90 }
|
|
|
91 occurrences=count["11"]+count["01"]+count["1N"]+count["1"]
|
|
|
92 fullycalled=count["11"]+count["01"]+count["00"]+count["1"]+count["0"]
|
|
|
93 print $1,$2,$3,$4,$5,$6,$7,$8,occurrences,occurrences/totalnormals,fullycalled,fullycalled/totalnormals,count["00"],count["01"],count["11"],count["0N"],count["1N"],count["NN"],count["0"],count["1"],count["N"]
|
|
|
94 }
|
|
|
95 }END{
|
|
|
96
|
|
|
97
|
|
|
98 }' $output_all > "${output_all}-counted"
|
|
|
99
|
|
|
100
|
|
|
101 # this counted file is the final output file
|
|
|
102 rm $output_all
|
|
|
103 mv "${output_all}-counted" $output_all
|
|
|
104
|
|
|
105
|
|
|
106
|
|
|
107 ### filter out variants occurring in more than <threshold> of the background genomes
|
|
|
108 # if total of columns containing a 1 (01,11,1N,1) is >= threshold
|
|
|
109 awk 'BEGIN{
|
|
|
110 FS="\t";
|
|
|
111 OFS="\t";
|
|
|
112 }{
|
|
|
113 if(FNR==1){
|
|
|
114 print $0
|
|
|
115 }
|
|
|
116 if(FNR>1){
|
|
|
117 if($9 < "'"$threshold"'" )
|
|
|
118 print $0
|
|
|
119 }
|
|
|
120 }END{}' $output_all > $output_filtered
|
|
|
121
|
|
|
122
|
|
|
123 awk 'BEGIN{
|
|
|
124 FS="\t";
|
|
|
125 OFS="\t";
|
|
|
126 threshold="'"${thresholdhc}"'"+0
|
|
|
127 }{
|
|
|
128 if(FNR==1)
|
|
|
129 print $0
|
|
|
130 else if($11 >= threshold)
|
|
|
131 print $0
|
|
|
132
|
|
|
133 }END{}' $output_filtered > "output_filtered_highconf.tsv"
|
|
|
134
|
|
|
135
|
|
|
136
|
|
|
137
|
|
|
138
|
|
|
139
|
|
|
140
|
|
|
141
|
|
2
|
142
|
