VariantFilterScripts/xCompoundHetFilters.sh at master · AshKernow/VariantFilterScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#!/bin/bash
#$ -cwd -l mem=1G,time=1:: -N FilterCHet


## use "0" for missing parent
ArrNum=1
#get arguments
while getopts v:t:n:a:p: opt; do
    case "$opt" in
        v) VcfFil="$OPTARG";;
        t) FamFil="$OPTARG";;
        n) DirPre="$OPTARG";;
        a) ArrNum="$OPTARG";;
        p) AddPrm="$OPTARG";;
        #H) echo "$usage"; exit;;
    esac
done

echo "Starting"

FiltScrDir="/home/local/ARCS/ads2202/scripts/Filtering_scripts"

VcfFil=`readlink -f $VcfFil`
FamFil=`readlink -f $FamFil`

FamNam=`cut -f 1 $FamFil | head -n $ArrNum | tail -n 1`
ModNam=`cut -f 2 $FamFil | head -n $ArrNum | tail -n 1`
PatNam=`cut -f 3 $FamFil | head -n $ArrNum | tail -n 1`
MatNam=`cut -f 4 $FamFil | head -n $ArrNum | tail -n 1`
FilPrm=`cut -f 5 $FamFil | head -n $ArrNum | tail -n 1`

if [[ $FamNam == [0-9]* ]]; then FamNam=Fam$FamNam; fi

DirNam=$FamNam
if [[ -n $DirPre ]]; then DirNam=$DirPre"_"$DirNam; fi
mkdir -p $DirNam
cd $DirNam

LogFil=$FamNam.compound_heterozygous.log

echo "Start :`date`"
echo "VCF file: $VcfFil"
echo "Filter file: $FamFil"
if [[ ! -e $VcfFil ]] | [[ ! -e $FamFil ]]; then "Missing/Incorrect required arguments"; exit; fi


#get paternally interited allele
PatPrm=$FilPrm
if [[ "$PatNam" != "0" ]]; then
    PatPrm=`echo $PatPrm | sed s/"--het "/"--het $PatNam,"/`
fi
if [[ "$MatNam" != "0" ]]; then
    if [[ "$PatPrm" == *--ref* ]]; then
        PatPrm=`echo $PatPrm | sed s/"--ref "/"--ref $MatNam,"/`
    else
        PatPrm=$PatPrm" --ref "$MatNam
    fi
fi
CMD="$FiltScrDir/ExmFilt.CustomGenotype.py -v $VcfFil -o $FamNam.tempheppat $PatPrm -P -f 0.03"
if [[ ! -z $AddPrm ]]; then CMD=$CMD" $AddPrm"; fi
echo $CMD
echo `date`
eval $CMD


MatPrm=$FilPrm
if [[ "$MatNam" != "0" ]]; then
    MatPrm=`echo $MatPrm | sed s/"--het "/"--het $MatNam,"/`
fi
if [[ "$PatNam" != "0" ]]; then
    if [[ "$MatPrm" == *--ref* ]]; then
        MatPrm=`echo $MatPrm | sed s/"--ref "/"--ref $PatNam,"/`
    else
        MatPrm=$MatPrm" --ref "$PatNam
    fi
fi
CMD="$FiltScrDir/ExmFilt.CustomGenotype.py -v $VcfFil -o $FamNam.temphepmat $MatPrm -P -f 0.03"
if [[ ! -z $AddPrm ]]; then CMD=$CMD" $AddPrm"; fi
echo $CMD
echo `date`
eval $CMD

echo "Run R script to find compound het"
echo `date`

R --vanilla <<RSCRIPT
options(stringsAsFactors=F)

mathet <- read.delim("$FamNam.temphepmat.tsv")
pathet <- read.delim("$FamNam.tempheppat.tsv")

mathet <- mathet[mathet[,"Gene"]%in%pathet[,"Gene"],]
pathet <- pathet[pathet[,"Gene"]%in%mathet[,"Gene"],]

comphet <- rbind(mathet, pathet)
comphet <- comphet[order(comphet[,"Chromosome"], comphet[,"Position"]),]
write.table(comphet, "$FamNam.compound_heterozygous.tsv", col.names=T, row.names=F, quote=F, sep="\t")
RSCRIPT

echo "Paternal Filtering" > $LogFil
echo "---------------------------------------------------" >> $LogFil
cat $FamNam.tempheppat.log >> $LogFil
echo "" >> $LogFil
echo "===================================================" >> $LogFil
echo "" >> $LogFil

echo "Maternal Filtering" >> $LogFil
echo "---------------------------------------------------" >> $LogFil
cat $FamNam.temphepmat.log >> $LogFil
echo "" >> $LogFil
echo "===================================================" >> $LogFil
echo "" >> $LogFil


rm -rf *temphep*
LEN=`cat $FamNam.compound_heterozygous.tsv | wc -l`
if [[ $LEN -gt 1 ]]; then
    CMD="nohup $FiltScrDir/xAnnotateVariantTSV.sh -i $FamNam.compound_heterozygous.tsv &"
    echo $CMD >> $LogFil
    eval $CMD
fi

GENLEN=`cut -f 6 $FamNam.compound_heterozygous.tsv | tail -n +2 | sort | uniq | wc -l`


echo "Final Number of variants: $LEN in $GENLEN genes" >> $LogFil
echo "" >> $LogFil
echo "===================================================" >> $LogFil

echo "Finished `date`"