Eric Coissac
195739b5f8
Realize the annotation on the normalized chromosome Former-commit-id: 73e4e016d7bfe23fc92e4071efb0a4629eb8e21f Former-commit-id: d44249c0d078d475ca80386bb428fb660289b0dd
68 lines
1.8 KiB
Bash
68 lines
1.8 KiB
Bash
#!/bin/bash
|
|
|
|
source "${THIS_DIR}/../../../scripts/bash_init.sh"
|
|
|
|
SELECTIR="${PROG_DIR}/../../normalize/lib/selectIR.py"
|
|
|
|
function lookForIR {
|
|
|
|
local QUERY="$1"
|
|
local MATCHES=$(basename ${QUERY})
|
|
MATCHES="${MATCHES/.*/}.matches"
|
|
|
|
local REPEATS="${MATCHES/.*/}.repseek"
|
|
|
|
# Blast columns:
|
|
# query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
|
|
# We keep blast matches if :
|
|
# The match is longer than 1000
|
|
# The identity is higher than 80%
|
|
#
|
|
# The match file has the following format:
|
|
# LSC/SSC begin end same_strand=1/diff_strand=0
|
|
|
|
loginfo "Locating SSC and LSC by similarity..."
|
|
blastn -db ${SCDB} \
|
|
-query ${QUERY} \
|
|
-outfmt 6 \
|
|
-max_target_seqs 10000 | \
|
|
awk '($4 > 1000) && ($3>80) { \
|
|
SAME=(($7 < $8) && ($9 < $10)) || (($7 > $8) && ($9 > $10)); \
|
|
if ($7 < $8) \
|
|
{print substr($2,1,3),$7,$8,SAME} \
|
|
else \
|
|
{print substr($2,1,3),$8,$7,SAME}}' | \
|
|
sort -nk 2 > ${MATCHES}
|
|
loginfo "Done"
|
|
|
|
|
|
loginfo "Looking for long inverted repeats..."
|
|
repseek -c -p 0.001 -i ${QUERY} 2>> /dev/null > ${REPEATS}
|
|
loginfo " --> $(wc -l ${REPEATS} | awk '{print $1}') repeats identified"
|
|
loginfo "Done"
|
|
|
|
loginfo "Marking and selecting the best inverted repeat..."
|
|
local IR=( $(${SELECTIR} ${MATCHES} ${REPEATS}) )
|
|
loginfo "Done"
|
|
|
|
loginfo " --> IR size : IRa = ${IR[5]} / IRb = ${IR[7]}"
|
|
loginfo " --> IR Score: ${IR[8]}"
|
|
|
|
let "deltaIR=${IR[5]} - ${IR[7]}"
|
|
|
|
if (( $deltaIR < -10 )) || (( $deltaIR > 10 )); then
|
|
logwarning "Differences between IR lengths ($deltaIR) is greater than 10"
|
|
fi
|
|
|
|
|
|
echo "${IR[@]}"
|
|
}
|
|
|
|
SCDB="${IR_DATA_DIR}/SC_RefDB"
|
|
|
|
if [[ ! "$1" =~ ^/ ]]; then
|
|
QUERY="${CALL_DIR}/$1"
|
|
else
|
|
QUERY="$1"
|
|
fi
|