Skip to content

C
CBBC-paralogs
Commit 1f09d480 by Edoardo SARTI
Options

table1 script updated

parent 530735fc
Showing with 92 additions and 124 deletions
scripts/annexes/CB_proteins.txt 0 → 100644
FBA
GAPDH
PGK
PRK
RPE
RPI
RubiscoLS
RubiscoSS
TPI
TRK
FBPase
SBPase
scripts/annexes/CB_proteins_an.txt 0 → 100644
FBA
GAPDH
PGK
PRK
RPE
RPI
TPI
TRK
FBPase
SBPase
scripts/annexes/FBPase_domains.txt 0 → 100644
FBPase_PF00316
FBPase_PF18913
scripts/annexes/GAPDH_domains.txt 0 → 100644
GAPDH_PF00044
GAPDH_PF02800
scripts/annexes/RubiscoLS_domains.txt 0 → 100644
RubiscoLS_PF00016
RubiscoLS_PF02788
scripts/annexes/SBPase_domains.txt 0 → 100644
SBPase_PF00316
SBPase_PF18913
scripts/annexes/TRK_domains.txt 0 → 100644
TRK_PF00456
TRK_PF02779
TRK_PF02780
scripts/annexes/x.tmp.txt 0 → 100644
SBPase_PF00316
SBPase_PF18913
scripts/get_all_up_subtree.py 0 → 100644
import sys
import ete3
tree_fn = sys.argv[1]
node_name = sys.argv[2]
t = ete3.Tree(tree_fn)
if node_name == "root":
subtree_root = t
else:
subtree_root = t.search_nodes(name=node_name)[0]
for l in [x for x in subtree_root]:
n = l.name.split()[0]
up = l.name.split("|")[1]
print(n)
scripts/table_S1.sh
for CBP in `cat CB_proteins.txt` INPUT_DIR="PV_trees/input_sequences_for_trees"
TREES_DIR="PV_trees/PV_annotated_trees"
SCRIPTS_DIR="scripts"
ANNEXES_DIR="scripts/annexes"
TMP_DIR="tmp/"
GROUND_TRUTH_DIR="PV_trees/PV_annotated_trees/ground_truth"
for CBP in `cat ${ANNEXES_DIR}/CB_proteins_an.txt`
do do
grep ${CBP} PF_domains > ${CBP}_domains.txt cat ${ANNEXES_DIR}/${CBP}_domains.txt > ${TMP_DIR}/x.tmp.txt
if [[ "${CBP}" == "GAPDH" ]]
then
for x in "GAPDH_PF00044" "GAPDH_PF02800"
do
grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus2_exp.fa | awk 'length($0)>2' | sort | uniq > f0.txt
P_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
P_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
P_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all paralogs predicted by Profileview
grep ${CBP} PF_domains > ${CBP}_domains.txt
python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_monodomain/${x}/${x}_plus2.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus2_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
PR_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
PR_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
PR_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all predicted paralogs having functional annotations
grep -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
F_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
F_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}"
done
else
for x in `cat ${CBP}_domains.txt`
do
# Count all paralogs in UniProtKB sequences
#grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | awk 'substr($0,1,1)==">"{h=$0} substr($0,1,1)!=">"{a[$1]=h} END{for (i in a) {print a[i]}}' | sort | uniq > f0.txt
grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | sort | uniq > f0.txt
P_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
P_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
P_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all paralogs predicted by Profileview
grep ${CBP} PF_domains > ${CBP}_domains.txt
#python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_monodomain/${x}/${x}_plus.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | awk 'substr($0,1,1)==">"{h=$0} substr($0,1,1)!=">"{a[$1]=h} END{for (i in a) {print a[i]}}' | sort | uniq > f0.txt
python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_monodomain/${x}/${x}_plus.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
#python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_monodomain/${x}/${x}_plus.nhx root > f1.txt; grep "CHLRE\|VOLCA\|ARATH" f1.txt | awk 'length($0)>2' | sort | uniq > f0.txt
PR_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
PR_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
PR_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all predicted paralogs having functional annotations
# grep -A1 -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | awk 'substr($0,1,1)==">"{h=$0} substr($0,1,1)!=">"{a[$1]=h} END{for (i in a) {print a[i]}}' > f0.txt
grep -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
F_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
F_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}"
done
if [[ "${CBP}" == "PGK" ]] if [[ "${CBP}" == "PGK" ]]
then then
x="PGK_PGAbinding" echo "PGK_PGAbinding" >> ${TMP_DIR}/x.tmp.txt
grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | sort | uniq > f0.txt elif [[ "${CBP}" == "TRK" ]]
P_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
P_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
P_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all paralogs predicted by Profileview
grep ${CBP} PF_domains > ${CBP}_domains.txt
python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_mixed/PGK_PGAbinding.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
PR_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
PR_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
PR_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
grep -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
F_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
F_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}"
fi
if [[ "${CBP}" == "TRK" ]]
then then
x="TRK_multi" echo "TRK_multi" >> ${TMP_DIR}/x.tmp.txt
grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | sort | uniq > f0.txt elif [[ "${CBP}" == "RPI" ]]
P_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
P_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
P_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
# Count all paralogs predicted by Profileview
grep ${CBP} PF_domains > ${CBP}_domains.txt
python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_mixed/TRK_multi.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
PR_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
PR_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
PR_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
grep -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
F_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'`
F_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'`
echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}"
fi
if [[ "${CBP}" == "RPI" ]]
then then
for x in "RPI_isomerase" "RPI_dimerase" "RPI_isomerase_dimerase" echo "RPI_isomerase" >> ${TMP_DIR}/x.tmp.txt
echo "RPI_dimerase" >> ${TMP_DIR}/x.tmp.txt
echo "RPI_isomerase_dimerase" >> ${TMP_DIR}/x.tmp.txt
fi
for x in `cat ${TMP_DIR}/x.tmp.txt`
do do
grep -A1 "CHLRE\|VOLCA\|ARATH" /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa | awk 'length($0)>2' | sort | uniq > f0.txt grep -A1 "CHLRE\|VOLCA\|ARATH" ${INPUT_DIR}/${CBP}/${CBP}.fa | awk 'length($0)>2' | sort | uniq > ${TMP_DIR}/f0.txt
P_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'` P_CHLRE=`grep "CHLRE" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
P_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'` P_VOLCA=`grep "VOLCA" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
P_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'` P_ARATH=`grep "ARATH" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
# Count all paralogs predicted by Profileview # Count all paralogs predicted by Profileview
grep ${CBP} PF_domains > ${CBP}_domains.txt python3 ${SCRIPTS_DIR}/get_all_up_subtree.py ${TREES_DIR}/${x}/${x}.nhx root > ${TMP_DIR}/f1.txt; grep -A1 -f ${TMP_DIR}/f1.txt ${INPUT_DIR}/${CBP}/${CBP}.fa > ${TMP_DIR}/f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" ${TMP_DIR}/f2.txt | awk 'length($0)>2{print $1}' | sort | uniq > ${TMP_DIR}/f0.txt
python3 get_all_up_subtree.py /workspaces/2024_PV_on_CB/data/PV_trees/from_mixed/${x}.nhx root > f1.txt; grep -A1 -f f1.txt /workspaces/2024_PV_on_CB/data/CB_sequences/${CBP}/${CBP}_plus_exp.fa > f2.txt; grep -A1 "CHLRE\|VOLCA\|ARATH" f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt
PR_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'` PR_CHLRE=`grep "CHLRE" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
PR_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'` PR_VOLCA=`grep "VOLCA" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
PR_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'` PR_ARATH=`grep "ARATH" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
grep -f /workspaces/2024_PV_on_CB/data/experimental_evidence/${CBP}_class_all_primaryacc_names.txt f2.txt | awk 'length($0)>2' | sort | uniq > f0.txt # Count all predicted paralogs having functional annotations
F_CHLRE=`grep "CHLRE" f0.txt | wc -l | awk '{print $1}'` grep -f ${GROUND_TRUTH_DIR}/${CBP}_class_all_primaryacc_names.txt ${TMP_DIR}/f2.txt | awk 'length($0)>2{print $1}' | sort | uniq > ${TMP_DIR}/f0.txt
F_VOLCA=`grep "VOLCA" f0.txt | wc -l | awk '{print $1}'` F_CHLRE=`grep "CHLRE" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" f0.txt | wc -l | awk '{print $1}'` F_VOLCA=`grep "VOLCA" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
F_ARATH=`grep "ARATH" ${TMP_DIR}/f0.txt | wc -l | awk '{print $1}'`
echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}" echo "${x} | ${P_CHLRE} ${PR_CHLRE} ${F_CHLRE} | ${P_VOLCA} ${PR_VOLCA} ${F_VOLCA} | ${P_ARATH} ${PR_ARATH} ${F_ARATH}"
done done
fi
fi
done done
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment