First version of if inverted repeat annotation tools
Former-commit-id: 9d13deae459b086f11ba4f0f82aac899b4ff0fd9 Former-commit-id: 74a7f07d3fb26f6ce12808b8fe2dc7f176181285
This commit is contained in:
@ -3,104 +3,37 @@
|
||||
# NORMALISATION D'UN PLASTIDE
|
||||
#
|
||||
#========================================================================================
|
||||
# Ce programme dispose de 4 fonctions pour traiter les donnees fasta issues de genbank
|
||||
# - seqlength : compte le nombre de paire de base du fichier
|
||||
# ex : seqlength $1
|
||||
#
|
||||
# - cutseq : permet de couper un morceau de la sequence
|
||||
# cutseq [x] [y]
|
||||
# [x] : coordonne du debut de la sequence a couper
|
||||
# [y] : coordonne de la fin de la sequence a couper
|
||||
# ex : cutseq $1 10 100
|
||||
# Normalize the way the chloroplaste genome sequence is linearized in the fasta file
|
||||
# The normalized sequence is:
|
||||
#
|
||||
# LSC + IRB + SSC + IRA
|
||||
#
|
||||
# - revcomp : donne le brin reverse
|
||||
# ex : $1 | revcomp
|
||||
# The SSC and LSC are approximatively mapped by similarity with a reference database
|
||||
# Inverted repeats (IRs) are identified for maximizing the segregation between
|
||||
# LSC and SSC match
|
||||
#
|
||||
#
|
||||
# - formatfasta : permet de coller a la suite plusieurs morceaux de sequence au moment de
|
||||
# la reecriture
|
||||
# - joinfasta : enleve les titres au moment de la reecriture du fichier et renvoie les
|
||||
# informations dans la fonction formatfasta
|
||||
# ex : joinfasta $1
|
||||
# go_normalize.sh <FASTAFILE>
|
||||
#
|
||||
#========================================================================================
|
||||
# Pour lancer le programme, utiliser les commandes :
|
||||
# chmod +x normalize_plastid.sh
|
||||
#./normalize_plastid.sh [fichier].fasta
|
||||
# - <FASTAFILE> : The fasta file containing the genome to normalize
|
||||
#
|
||||
# ex : seqlength $1
|
||||
# Results are printed to the standart output
|
||||
#
|
||||
# cutseq $1 [x] [y]
|
||||
# [x]:coordonne du debut [y]:coordonne de la fin de la sequence a couper
|
||||
# ex : cutseq $1 10 100
|
||||
#
|
||||
# ex : $1 | revcomp
|
||||
#
|
||||
# ex : joinfasta $1
|
||||
#========================================================================================
|
||||
|
||||
# -- CAUTION -- Works as long than the script
|
||||
# is not called through a symlink
|
||||
SCRIPT_DIR="$(dirname ${BASH_SOURCE[0]})"
|
||||
source "${SCRIPT_DIR}/../../../scripts/bash_init.sh"
|
||||
source ${SCRIPT_DIR}/../lib/lookforIR.lib.sh
|
||||
|
||||
function lookForIR {
|
||||
|
||||
local QUERY="$1"
|
||||
local MATCHES=$(basename ${QUERY})
|
||||
MATCHES="${MATCHES/.*/.matches}"
|
||||
|
||||
local REPEATS="${MATCHES/.*/.repseek}"
|
||||
|
||||
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=( $(${PROG_DIR}/selectIR.py ${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[@]}"
|
||||
}
|
||||
|
||||
pushTmpDir ORG.normalize
|
||||
|
||||
SCDB="${IR_DATA_DIR}/SC_RefDB"
|
||||
QUERY="${CALL_DIR}/$1"
|
||||
MATCHES="${1/.*/.matches}"
|
||||
REPEATS="${1/.*/.repseek}"
|
||||
|
||||
tmpfasta1="tmp_$$_1.fasta"
|
||||
tmpfasta2="tmp_$$_2.fasta"
|
||||
|
||||
|
||||
openLogFile "${QUERY/.*/.log}"
|
||||
|
||||
loginfo "Computing the genome size..."
|
||||
genome_length=$(seqlength $QUERY)
|
||||
loginfo " --> $genome_length bp"
|
||||
@ -136,7 +69,7 @@ pushTmpDir ORG.normalize
|
||||
rm -f ${tmpfasta1}
|
||||
QUERY=${tmpfasta2}
|
||||
|
||||
loginfo "Recompute location of the IR..."
|
||||
loginfo "Recomputing location of the IR..."
|
||||
declare -a IR=( $(lookForIR ${QUERY}) )
|
||||
loginfo "Done"
|
||||
|
||||
@ -224,7 +157,7 @@ pushTmpDir ORG.normalize
|
||||
fi
|
||||
|
||||
# Merges the four parts of the genome.
|
||||
cat ${tmpSSC} ${tmpIR1} ${tmpLSC} ${tmpIR2} | joinfasta
|
||||
cat ${tmpLSC} ${tmpIR2} ${tmpSSC} ${tmpIR1} | joinfasta
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,135 +0,0 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import sys
|
||||
|
||||
data = open(sys.argv[1])
|
||||
repeats = open(sys.argv[2])
|
||||
|
||||
chloro = {'LSC' : [], 'SSC' : [] }
|
||||
chlorosize =0
|
||||
|
||||
for line in data:
|
||||
parts = line.strip().split()
|
||||
if len(parts) >= 4:
|
||||
single = parts[0]
|
||||
begin = int(parts[1])
|
||||
end = int(parts[2])
|
||||
direction = int(parts[3])
|
||||
|
||||
|
||||
if direction==0:
|
||||
direction=-1
|
||||
|
||||
if end > chlorosize:
|
||||
extsize = end - chlorosize
|
||||
chloro['LSC'].extend([0] * extsize)
|
||||
chloro['SSC'].extend([0] * extsize)
|
||||
chlorosize=len(chloro['LSC'])
|
||||
|
||||
begin-=1
|
||||
|
||||
chr = chloro[single]
|
||||
|
||||
for p in range(begin,end):
|
||||
chr[p]+=direction
|
||||
|
||||
maxSSC = float(max(abs(n) for n in chloro['SSC']))
|
||||
maxLSC = float(max(abs(n) for n in chloro['LSC']))
|
||||
|
||||
chloro['SSC']=[n / maxSSC for n in chloro['SSC']]
|
||||
chloro['LSC']=[n / maxLSC for n in chloro['LSC']]
|
||||
|
||||
|
||||
scoreMax=0
|
||||
imax = len(chloro['LSC'])
|
||||
|
||||
for line in repeats:
|
||||
parts = line.strip().split()
|
||||
|
||||
pos1 = int(parts[1]) -1
|
||||
len1 = int(parts[3])
|
||||
|
||||
pos2 = int(parts[2]) -1
|
||||
len2 = int(parts[4])
|
||||
|
||||
c_begin = min(pos1 + len1,imax)
|
||||
c_end = min(pos2,imax)
|
||||
o_max = min(pos1 ,imax)
|
||||
o_min = min(pos2 + len2, imax)
|
||||
|
||||
c_lsc = sum(abs(chloro['LSC'][n]) for n in range(c_begin,c_end))
|
||||
c_ssc = sum(abs(chloro['SSC'][n]) for n in range(c_begin,c_end))
|
||||
|
||||
o_lsc = sum(abs(chloro['LSC'][n]) for n in range(0,o_max))
|
||||
o_ssc = sum(abs(chloro['SSC'][n]) for n in range(0,o_max))
|
||||
|
||||
o_lsc += sum(abs(chloro['LSC'][n]) for n in range(o_min,len(chloro['LSC'])))
|
||||
o_ssc += sum(abs(chloro['SSC'][n]) for n in range(o_min,len(chloro['SSC'])))
|
||||
|
||||
c = float(c_lsc + c_ssc)
|
||||
o = float(o_lsc + o_ssc)
|
||||
|
||||
if c > 0:
|
||||
c_lsc /= c
|
||||
c_ssc /= c
|
||||
|
||||
if o > 0:
|
||||
o_lsc /= o
|
||||
o_ssc /= o
|
||||
|
||||
score = ((c_lsc - c_ssc) ** 2 + (o_lsc - o_ssc) ** 2) / 2.0
|
||||
|
||||
# print >>sys.stderr,"c.lsc = %f c.ssc = %f o.lsc = %f o.ssc = %f score = %6.4f (len=%d)" % (c_lsc,c_ssc,o_lsc,o_ssc,score,len1)
|
||||
|
||||
if (score > scoreMax):
|
||||
scoreMax = score
|
||||
pos1Max = pos1
|
||||
pos2Max = pos2
|
||||
len1Max = len1
|
||||
len2Max = len2
|
||||
|
||||
c_begin = min(pos1Max + len1Max,imax)
|
||||
c_end = min(pos2Max,imax)
|
||||
o_max = min(pos1Max,imax)
|
||||
o_min = min(pos2Max + len2Max,imax)
|
||||
|
||||
c_lsc = sum(chloro['LSC'][n] for n in range(c_begin,c_end))
|
||||
c_ssc = sum(chloro['SSC'][n] for n in range(c_begin,c_end))
|
||||
|
||||
o_lsc = sum(chloro['LSC'][n] for n in range(0,o_max))
|
||||
o_ssc = sum(chloro['SSC'][n] for n in range(0,o_max))
|
||||
|
||||
o_lsc += sum(chloro['LSC'][n] for n in range(o_min,len(chloro['LSC'])))
|
||||
o_ssc += sum(chloro['SSC'][n] for n in range(o_min,len(chloro['SSC'])))
|
||||
|
||||
if abs(c_lsc) > abs(c_ssc):
|
||||
center = "LSC"
|
||||
dcenter= "+" if c_lsc > 0 else "-"
|
||||
else:
|
||||
center = "SSC"
|
||||
dcenter= "+" if c_ssc > 0 else "-"
|
||||
|
||||
if abs(o_lsc) > abs(o_ssc):
|
||||
out = "LSC"
|
||||
dout = "+" if o_lsc > 0 else "-"
|
||||
else:
|
||||
out = "SSC"
|
||||
dout = "+" if o_ssc > 0 else "-"
|
||||
|
||||
|
||||
|
||||
sys.stdout.write("%s %s %s %s %d %d %d %d %6.5f\n" % (center,
|
||||
dcenter,
|
||||
out,
|
||||
dout,
|
||||
pos1Max + 1,
|
||||
len1Max,
|
||||
pos2Max + 1,
|
||||
len2Max,
|
||||
scoreMax))
|
||||
|
||||
|
||||
|
||||
#for p in range(chlorosize):
|
||||
# sys.stdout.write("%d %d %d\n" % (p,chloro['SSC'][p],chloro['LSC'][p]))
|
||||
|
||||
Reference in New Issue
Block a user