PhyloSkims/annotate
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Initial commit of ORG.Annotate

Browse Source 
Former-commit-id: 8c4c1d0391ad06ed67e748e0f837399e9b3de7ad
Former-commit-id: 5af9ed909a6f6bed7bae86c77bd770ae4b5fc14e
This commit is contained in:
Eric Coissac
2015-05-07 16:41:43 +02:00
commit 7eb55b5d0b
9 changed files with 1892 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

96
annotate_plastid.sh Executable file
View File

@ -0,0 +1,96 @@
#!/bin/bash
#
#
export ORGANNOT_HOME=`dirname $0`
REPSEEK=${ORGANNOT_HOME}/repseek
SUMATRA=${ORGANNOT_HOME}/sumatra
ARAGORN=${ORGANNOT_HOME}/aragorn
WRAPARAGORN=${ORGANNOT_HOME}/aragorn_wrapper.awk
ECOFIND=${ORGANNOT_HOME}/ecofind
function annotateCAU {
QUERY="$$.query.fasta"
echo $1 | sed 's/&/ /' | tr '@' '\n' > ${QUERY}
${SUMATRA} -d -n ${QUERY} $2 2> /dev/null | \
awk ' {n[$2]+=1;d[$2]+=$3} \
END {for (i in n) \
print i, n[i],d[i], d[i]/n[i]\
}' | \
sort -rnk4 | \
egrep '^trn(I|M|fM)' | \
tail -1 | \
awk '{print $1,$NF}'
rm -rf ${QUERY}
}
function gffTRNA {
${ARAGORN} -w -io -seq $3 | awk -v gid=${1} -f ${WRAPARAGORN}
}
# s'alimente avec un fichier.fasta
# $3 : nb de caractere du fichier, t : nb de caractere du titre,
# $1+1 : nb de retour chariot du fichier
function seqlength {
cat $1 | \
wc |\
awk -v t="`head -1 $1 | wc -c`" '{print $3 - t - $1 + 1}'
}
# recupere les informations issues du programme repseek avec l'origine des deux
# IR et leur taille
function lookforIR {
${REPSEEK} -c -p 0.001 $1 | \
grep 'Distant.inv' | \
sort -n -k4 | \
tail -1 | \
awk '{print $7}' | \
sed 's/-/ /g'
}
# recupere le nom de la sequence analyse
function seqName {
head -n1 $1| \
awk '{print $1}' | \
sed 's/^>//' | \
sed -E 's/.*\|([^|]+)\|/\1/'
}
# cree un resume du fichier analyse au format gff
# ex : GFF (NC_*** Repseek IR1 start end . + . )
function gffIR {
lseq=$2
nom=$1
lookforIR $3 | \
awk -v nom="$nom" -v lseq="$lseq" \
'BEGIN {OFS="\t"} \
{ startIR1=$1; \
startIR2=$2; \
endIR1=startIR1 + $3 -1; \
endIR2=startIR2 + $3 -1; \
startSSC=1; \
endSSC=startIR1-1; \
startLSC=endIR1+1; \
endLSC=startIR2-1; \
\
print nom,"RepSeek","misc_feature",startSSC,endSSC,"\.","+","\.","ID=SSC;note=small single copy region";\
print nom,"RepSeek","repeat_region",startIR1,endIR1,"\.","+","\.","ID=IRA;note=inverted repeat A";\
print nom,"RepSeek","misc_feature",startLSC,endLSC,"\.","+","\.","ID=LSC;note=large single copy region";\
print nom,"RepSeek","repeat_region",startIR2,endIR2,"\.","-","\.","ID=IRB;note=inverted repeat B";\
}'
}
echo "##gff-version 3"
genome=$1
genome_name=`seqName $1`
genome_length=`seqlength $1`
gffIR ${genome_name} ${genome_length} ${genome}| grep -v '^ *$'
gffTRNA ${genome_name} ${genome_length} ${genome}| grep -v '^ *$'

229
aragorn_wrapper.awk Executable file
View File

@ -0,0 +1,229 @@
#!/usr/bin/awk -f
function genomeid() {
if (gid=="") {
gid="XXXXXXX";
}
return gid;
}
function home() {
"echo $ORGANNOT_HOME" | getline homedir;
return homedir;
}
function prog(program) {
return home() "/" program;
}
function trnalib(prognam) {
return home() "/lib/trnaCAU.ref.fasta";
}
function awkPID() {
"bash -c 'echo $PPID'" | getline pid
return pid
}
function tmpName(prefix,suffix) {
return prefix awkPID() suffix
}
function rm(filename) {
system("rm -f " filename);
}
function printfasta(id,seq,filename) {
if (filename=="")
filename = tmpName(id "_",".fasta");
seqlen=length(seq);
print ">" id > filename;
for (i=1; i <= seqlen; i+=60)
print substr(seq,i,60) >> filename;
close(filename);
return filename;
}
function epissage(intron,seq) {
if (intron != "") {
l=length(intron);
intron=substr(intron,2,l-2);
split(intron,intronpos,",");
lseq=length(seq);
lintron2=lseq - intronpos[2] + 1;
seq = substr(seq,1,intronpos[1]) substr(seq,intronpos[2],lintron2);
}
return seq;
}
function patchtRNA(anticodon,trna,seq) {
if (anticodon == "cat") {
file=printfasta(trna "_" anticodon,seq,"");
command= prog("sumatra") " -d -n " file " " trnalib();
while ((command | getline output) > 0) {
split(output,field," ");
n[field[2]]++;
d[field[2]]=field[3];
}
close(command)
dmin=1;
for (i in n) {
dist=d[i]/n[i];
if (dist < dmin) {
dmin=dist;
trna=i;
}
}
}
else {
trna="trn" AA1[substr(trna,6,3)];
}
return trna;
}
function gene2product(gene) {
return "tRNA-" AA3[substr(gene,4,1)];
}
function gffTRNA(geneid,trna,loc,anti,intron,seq) {
if (loc ~ /^c/) {
complement="-";
loc=substr(loc,3,length(loc)-3);
}
else {
complement="+";
loc=substr(loc,2,length(loc)-2);
}
split(loc,pos,",");
anti=toupper(anti);
gsub("T","U",anti);
product=gene2product(trna);
printf("%s\taragorn\tgene\t%d\t%d\t.\t%s\t.\tID=genetrn%d;gbkey=Gene;gene=%s\n",
genomeid(),pos[1],pos[2],complement,geneid,trna);
printf("%s\taragorn\ttRNA\t%d\t%d\t.\t%s\t.\tID=trn%d;parent=genetrn%d;gbkey=tRNA;gene=%s;Note=anticodon: %s;product=%s\n",
genomeid(),pos[1],pos[2],complement,geneid,geneid,trna,anti,product);
if (intron=="") {
printf("%s\taragorn\texon\t%d\t%d\t.\t%s\t.\tID=exontrn%d;parent=trn%d;gbkey=tRNA;gene=%s;Note=anticodon: %s;product=%s\n",
genomeid(),pos[1],pos[2],complement,geneid,geneid,trna,anti,product);
}
else {
l=length(intron);
intron=substr(intron,2,l-2);
split(intron,intronpos,",");
lseq=pos[2]-pos[1]+1;
lintron2=lseq - intronpos[2] + 1;
seq = substr(seq,1,intronpos[1]) substr(seq,intronpos[2],lintron2);
printf("%s\taragorn\texon\t%d\t%d\t.\t%s\t.\tID=exontrn%da;parent=trn%d;gbkey=tRNA;gene=%s;Note=anticodon: %s;product=%s\n",
genomeid(),pos[1],pos[1]+intronpos[1]-2,complement,geneid,geneid,trna,anti,product);
printf("%s\taragorn\tintron\t%d\t%d\t.\t%s\t.\tID=exontrn%di;parent=trn%d;gbkey=intron;gene=%s;Note=anticodon: %s;product=%s\n",
genomeid(),pos[1]+intronpos[1]-1,pos[1]+intronpos[2]-2,complement,geneid,geneid,trna,anti,product);
printf("%s\taragorn\texon\t%d\t%d\t.\t%s\t.\tID=exontrn%db;parent=trn%d;gbkey=tRNA;gene=%s;Note=anticodon: %s;product=%s\n",
genomeid(),pos[1]+intronpos[2]-1,pos[2],complement,geneid,geneid,trna,anti,product);
}
}
BEGIN {
print ARGV[1];
AA1["Ala"]="A";
AA1["Cys"]="C";
AA1["Asp"]="D";
AA1["Glu"]="E";
AA1["Phe"]="F";
AA1["Gly"]="G";
AA1["His"]="H";
AA1["Ile"]="I";
AA1["Lys"]="K";
AA1["Leu"]="L";
AA1["Met"]="M";
AA1["Asn"]="N";
AA1["Pyl"]="O";
AA1["Pro"]="P";
AA1["Gln"]="Q";
AA1["Arg"]="R";
AA1["Ser"]="S";
AA1["Thr"]="T";
AA1["Sec"]="U";
AA1["Val"]="V";
AA1["Trp"]="W";
AA1["Tyr"]="Y";
AA3["A"]="Ala";
AA3["C"]="Cys";
AA3["D"]="Asp";
AA3["E"]="Glu";
AA3["F"]="Phe";
AA3["G"]="Gly";
AA3["H"]="His";
AA3["I"]="Ile";
AA3["K"]="Lys";
AA3["L"]="Leu";
AA3["M"]="Met";
AA3["N"]="Asn";
AA3["O"]="Pyl";
AA3["P"]="Pro";
AA3["Q"]="Gln";
AA3["R"]="Arg";
AA3["S"]="Ser";
AA3["T"]="Thr";
AA3["U"]="Sec";
AA3["V"]="Val";
AA3["W"]="Trp";
AA3["Y"]="Tyr";
AA3["f"]="fMet";
}
/^>/ { id = substr($1,2);
}
/^[0-9]+ +genes found/ \
{ nbgene=$1
}
((geneid != "") && /^[0-9]+/ && ! /genes found/) \
{ seq=epissage(intron,seq);
trna=patchtRNA(anti,trna,seq);
# print geneid,trna,loc,anti,"'"intron"'",seq;
gffTRNA(geneid,trna,loc,anti,intron,seq);
seq=""
}
/^[0-9]+/ && ! /genes found/ \
{ geneid=$1;
trna =$2;
loc =$3;
lseq =$4;
x=$5;
split($5,intron_desc,"i");
anti =substr(intron_desc[1],2,3);
intron=intron_desc[2];
}
/^[^>0-9]/ \
{ seq=seq $1
}
END { seq=epissage(intron,seq);
trna=patchtRNA(anti,trna,seq);
# print geneid,trna,loc,anti,"'"intron"'",seq;
gffTRNA(geneid,trna,loc,anti,intron,seq);
}

BIN
lib/ncbi20150303.adx Normal file
View File

Binary file not shown.

1
lib/ncbi20150303.ndx.REMOVED.git-id Normal file
View File

@ -0,0 +1 @@
0b9d49e3330cd8e6abe0c8b010e9931f7b2f22f3

BIN
lib/ncbi20150303.rdx Normal file
View File

Binary file not shown.

1
lib/ncbi20150303.tdx.REMOVED.git-id Normal file
View File

@ -0,0 +1 @@
87cabbb1913329a9e693acad81654beb01f0e9bc

1263
lib/trnaCAU.ref.fasta Normal file
View File

File diff suppressed because it is too large Load Diff

278
normalize_plastid.sh Executable file
View File

@ -0,0 +1,278 @@
#!/bin/bash
#
# 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
#
#
# - revcomp : donne le brin reverse
# ex : $1 | revcomp
#
# - 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
#
#========================================================================================
# Pour lancer le programme, utiliser les commandes :
# chmod +x normalize_plastid.sh
#./normalize_plastid.sh [fichier].fasta
#
# ex : seqlength $1
#
# 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
#========================================================================================
FINDCDS=`dirname $0`/findcds
REPSEEK=`dirname $0`/repseek
# s'alimente avec un fichier.fasta
# $3 : nb de caractere du fichier, t : nb de caractere du titre,
# $1+1 : nb de retour chariot du fichier
function seqlength {
cat $1 | \
wc |\
awk -v t="`head -1 $1 | wc -c`" '{print $3 - t - $1 + 1}'
}
# selectionne une sequence parmi le fichier
# $2 : debut de la sequence a couper, $3 : fin de la sequence a couper
function cutseq {
awk -v from=$2 -v end=$3 'function printfasta(seq) { \
seqlen=length(seq); \
for (i=1; i <= seqlen; i+=60) \
print substr(seq,i,60); \
} \
\
/^>/ {print $0} \
! /^>/ {seq=seq$0} \
END {printfasta(substr(seq,from,end-from+1))}' $1
}
# donne le brin reverse de la sequence :
# la sous-fonction comp reecrit la sequence a l'envers
# la sous-fonction rev remplace les bases par leurs bases associe
# la sous-fonction revcomp reprend les deux precedente
function revcomp {
awk 'function printfasta(seq) { \
seqlen=length(seq); \
for (i=1; i <= seqlen; i+=60) \
print substr(seq,i,60); \
} \
function comp(seq) { \
"echo "seq" | tr acgtACGT tgcaTGCA " | getline res; \
return res; \
} \
function rev(seq) { \
"echo "seq" | rev " | getline res; \
return res; \
} \
function revcomp(seq) { \
res=rev(comp(seq)); \
return res; \
} \
\
/^>/ {print $0} \
! /^>/ {seq=seq$0} \
END {printfasta(revcomp(seq))}' $1
}
function formatfasta {
awk 'function printfasta(seq) { \
seqlen=length(seq); \
for (i=1; i <= seqlen; i+=60) \
print substr(seq,i,60); \
} \
/^>/ { print $0 } \
! /^>/ { seq=seq $0 } \
END { printfasta(seq)}' $1
}
# colle bout a bout deux sequence en mettant le meme nombre de paire de base par ligne
# sur le fichier, ici regle a 60
# enleve les titres intermediaire entre deux sequences recollees si il y en a
function joinfasta {
awk '(NR==1 && /^>/) {print $0} \
! /^>/ {print $0}' $1 | \
formatfasta
}
# recupere les informations issues du programme repseek avec l'origine des deux
# IR et leur taille
function lookforIR {
repseek -c -p 0.001 $1| \
grep 'Distant.inv'| \
sort -n -k4 | \
tail -1 | \
awk '{print $7}' | \
sed 's/-/ /g'
}
# determine si le fragment analyse doit etre recolle en forward ou reverse dans le
# nouveau .fasta
function maxCDS {
${FINDCDS} -F $1 -c -l 150 | \
awk '/^[^#]/ && ($2 == "Watson") { Watson+=$5-$4+1} \
/^[^#]/ && ($2 == "Crick") {Crick+=$5-$4+1} \
END {print Watson - Crick}'
}
#
# Exemple results from repseek
#
# $1 $2 $3 $4 $5 $6 $7 $8 $9 $10 $11 $12
# Class pos_r1 pos_r2 len_r1 len_r2 Delta Seed ident score Rmean Rmod frac
# Distant.inv 86608 130934 25319 25319 19007 86608-130934-25319-2.01 100.000 25262.43 2.01 2 0.99
# Test where the sequence is cut
# Definie les variables utilisé : le debut, fin et taille des deux IR
genome=$1
genome_length=`seqlength $1`
IRS=(`lookforIR ${genome}`)
posIR1=${IRS[0]}
posIR2=${IRS[1]}
lenIR=${IRS[2]}
let "endIR2=$posIR2 + $lenIR - 1"
let "endIR1=$posIR1 + $lenIR - 1"
# Defini la coupe a adopter en fonction de :
# - la position de la fin de l'IR2 par rapport a la sequence total, pour identifier une
# coupe au sein d'une IR
# Le programme repseek considere toujours que la position maximal de la fin d'IR2 ne peut
# pas depasser celle de la sequence, que l'IR2 soit coupe ou non. Donc dans tous les cas
# on choisit de recouper la sequence a mi-distance entre la fin de l'IR1 et le debut de
# l'IR2
if (( endIR2 == genome_length )) ; then
tmpfasta1="tmp_$$_1.fasta"
tmpfasta2="tmp_$$_2.fasta"
# defini la localisation de la coupure entre les deux IR
let "posCut=($endIR1+$posIR2)/2"
# realise la coupure du fichier d'entre du nucleotide calcule jusqu'a la fin de la sequence
cutseq ${genome} ${posCut} ${genome_length} > ${tmpfasta1}
let "posCut=$posCut-1"
# realise la coupure du fichier d'entre du debut de la sequence jusqu'au nucleotide calcule
cutseq ${genome} 1 ${posCut} >> ${tmpfasta1}
# ces deux fragment sont rassembles dans un fichier temporaire
joinfasta ${tmpfasta1} > ${tmpfasta2}
rm -f ${tmpfasta1}
genome=${tmpfasta2}
# recalcul la nouvelle position des IR
IRS=(`lookforIR ${genome}`)
posIR1=${IRS[0]}
posIR2=${IRS[1]}
lenIR=${IRS[2]}
let "endIR2=$posIR2 + $lenIR - 1"
let "endIR1=$posIR1 + $lenIR - 1"
fi
tmpIR1="tmp_$$_IR1.fasta"
tmpIR2="tmp_$$_IR2.fasta"
#enregistre les deux fragments IR1 et IR2 complet
cutseq ${genome} ${posIR1} ${endIR1} > ${tmpIR1}
cutseq ${genome} ${posIR2} ${endIR2} > ${tmpIR2}
let "lenSC1=$posIR1 -1 + ($genome_length - endIR2)"
let "lenSC2=$posIR2 - $endIR1"
tmpLSC="tmp_$$_LSC.fasta"
tmpSSC="tmp_$$_SSC.fasta"
# Defini la coupe a adopter en fonction de :
# - la taille de la SC1 par rapport a la taille de la SC2, pour identifier une
# coupe au sein d'une SC. La coupe a lieu au sein de la SC1, le but est d'identifier la
# LSC et la SSC parmis les SC1 et la SC2
# si la SC1 est plus grande que la SC2, alors la SC1 est la LSC et la coupe a eu lieu
# dans la LSC
if (( lenSC1 > lenSC2 )); then
# defini le debut de la LSC
let "beginLSC=$endIR2+1"
cutseq ${genome} ${beginLSC} ${genome_length} > ${tmpLSC}
# defini la fin de la LSC
let "endLSC=$posIR1-1"
cutseq ${genome} 1 ${endLSC} >> ${tmpLSC}
tmpfasta1="tmp_$$_1.fasta"
# rejoint les deux morceaux pour former la LSC
joinfasta ${tmpLSC} > ${tmpfasta1}
mv ${tmpfasta1} ${tmpLSC}
# donc la SC2 est la SSC,
# definit l'origine et la fin a couper pour avoir le fragment SSC
let "beginSSC=$endIR1+1"
let "endSSC=$posIR2-1"
cutseq ${genome} ${beginSSC} ${endSSC} > ${tmpSSC}
tmp=${tmpIR1}
tmpIR1=${tmpIR2}
tmpIR2=${tmp}
else
# sinon la SC2 est la LSC, et la coupe a eu lieu dans la SSC
# definit l'origine et la fin a couper pour avoir le fragment LSC
let "beginLSC=$endIR1+1"
let "endLSC=$posIR2-1"
cutseq ${genome} ${beginLSC} ${endLSC} > ${tmpLSC}
# definit le debut de la SSC et coupe la premiere partie
let "beginSSC=$endIR2+1"
cutseq ${genome} ${beginSSC} ${genome_length} > ${tmpSSC}
# definit la fin de la SSC et coupe la seconde partie
let "endSSC=$posIR1-1"
cutseq ${genome} 1 ${endSSC} >> ${tmpSSC}
# joint les deux parties afin de reformer la SSC
joinfasta ${tmpSSC} > ${tmpfasta1}
mv ${tmpfasta1} ${tmpSSC}
fi
if [[ ! -z $tmpfasta2 ]]; then
rm -f $tmpfasta2
fi
# determine si les fragments doivent etre recolle en reverse ou forward
maxSSC=`maxCDS ${tmpSSC}`
# si maxSSC est negatif, le rapport Watson - Crick est negatif, le fragment est
# donc reverse
if (( maxSSC < 0 )); then
revcomp ${tmpSSC} > ${tmpfasta1}
mv ${tmpfasta1} ${tmpSSC}
fi
maxLSC=`maxCDS ${tmpLSC}`
# si maxLSC est negatif, le rapport Watson - Crick est negatif, le fragment est
# donc reverse
if (( maxLSC < 0 )); then
revcomp ${tmpLSC} > ${tmpfasta1}
mv ${tmpfasta1} ${tmpLSC}
fi
# Les quatre fragments sont recolle ensemble sans erreur de coupure et dans un ordre connus.
cat ${tmpSSC} ${tmpIR1} ${tmpLSC} ${tmpIR2} | joinfasta
exit 0

24
reannoteCAURefDB.sh Executable file
View File

@ -0,0 +1,24 @@
#!/bin/bash
SUMATRA=`dirname $0`/sumatra
function annotateCAU {
QUERY="$$.query.fasta"
echo $1 | sed 's/&/ /' | tr '@' '\n' > ${QUERY}
${SUMATRA} -d -n ${QUERY} $2 2> /dev/null | \
awk ' {n[$2]+=1;d[$2]+=$3} \
END {for (i in n) \
print i, n[i],d[i], d[i]/n[i]\
}' | \
sort -rnk4 | \
egrep '^trn(I|M|fM)' | \
tail -1 | \
awk '{print $1,$NF}'
rm -rf ${QUERY}
}
for seq in `awk '(on==1) && /^>/ {print "";on=0} /^>/ {printf("%s@",$0)} ! /^>/ {on=1;printf($0)}' $1 | tr ' ' '&'`; do
new=(`annotateCAU ${seq} $1`)
echo $seq | sed 's/&/ /g' | sed -E 's/>([^ ]+) />'${new[0]}' /' | tr '@' '\n'
done