Utilities/inversions_VCF.sh at master · pmorrell/Utilities · 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
#!/bin/bash -l

# Peter L. Morrell - 27 September 2024 - Falcon Heights, MN

set -e
set -o pipefail

# Usage: Provide a list of inversion positions in a BED file. The
# other required input is the named of a VCF file
# indexed by bcftools. The third should be a path to the output directory.

module load bcftools/1.16-gcc-8.2.0-5d4xg4y

BED_FILE=$1
VCF=$2
OUTDIR=$3
OUT_PREFIX=WBDC_GBS
SAMPLE=WBDC-355
MIN_INFORM=5


# Declare an array to store genomic positions
declare -a GENOMIC_POSITIONS


while read -r line; do
  # Split the line into fields
  IFS=$'\t' read -r chrom start end <<< "$line"
  # Add 1 to the start and end positions
  new_start=$((start + 1))
  new_end=$((end + 1))
  # Combine chromosome, start, and end into a single string
  genomic_position="${chrom}:${new_start}-${new_end}"
  # Append the genomic position to the array
  GENOMIC_POSITIONS+=("$genomic_position")
done < "$BED_FILE"


# Print the genomic positions
for position in "${GENOMIC_POSITIONS[@]}"; do
  echo "$position"
done


if [ ! -d "$OUTDIR" ]; then
    mkdir -p "$OUTDIR";
fi;


for i in "${GENOMIC_POSITIONS[@]}"
    do
    	INFORM_SNPS=$(bcftools query -r "$i" -s "$SAMPLE" -f '[%GT]\n' \
            "$VCF" | grep -E -c '1/1|0/1|1/0')
    		if "$INFORM_SNPS" gt "$MIN_INFORM"; then
        bcftools view -r "$i" "$VCF" -o "$OUTDIR"/"$OUT_PREFIX"_{"$i"}.vcf.gz -Oz
            else
            echo "There are no SNPs in the interval $i"
        fi
    done