orthofinderSCO_RBH_1wayBH
This teaching workflow combines single-copy orthologues, reciprocal best hits, and one-way best hits. Before running OrthoFinder, prepare the sequence data with the extract_isoform workflow, described separately This DIAMOND workflow uses peptide sequences.
Steps:
- OrthoFinder: only get Bfl genes in SCO result
- Run DIAMOND reciprocal best hits (RBH) between human TFs and the Bfl genes not assigned to single-copy orthogroups.
- Run DIAMOND one-way best hits between human TFs and the Bfl genes still unassigned after RBH.
- merge all 3 outputs (*Bfl_TF_list.txt)
Tutorial 1 (GOAL: Predict the TFs in Bfl using a TFs FASTA of Human from AnimalTFDB4.0)
Work directory
/data/core/jordan/Shared/518_projects/workshop_20231213/orthofinderSCO_RBH_1wayBH-Bfl
conda activate /data/core/jordan/CondaEnv/orthofinder1. OrthoFinder: only get Bfl genes in SCO result
Use more permissive MMseqs2 settings: -s 6 -e 1e-3 -k 5 --num-iterations 5.
vim ~/CondaEnv/orthofinder/bin/scripts_of/config.json
orthofinder -f orthofinder_input/ -S mmseqs -t 150OrthoFinder assigned 23208 genes (81.9% of total) to 2421 orthogroups. Fifty percent of all genes were in orthogroups with 18 or more genes (G50 was 18) and were contained in the largest 203 orthogroups (O50 was 203). There were 407 orthogroups with all species present and 179 of these consisted entirely of single-copy genes.grep -wf orthofinder_input/OrthoFinder/Results_Dec07/Orthogroups/Orthogroups_SingleCopyOrthologues.txt orthofinder_input/OrthoFinder/Results_Dec07/Phylogenetic_Hierarchical_Orthogroups/N0.tsv > ./orthofinderSCO.txt- -w: exact match
awk -F "\t" '{print $4}' orthofinderSCO.txt > orthofinderSCO_Bfl_TF_list.txt→ 179 Bfl genes
LOC118403119|Bfl
LOC118403141|Bfl
LOC118403181|Bfl
LOC118403236|Bfl
LOC118403522|Bfl
LOC118403701|Bfl
LOC118403813|Bfl
LOC118403816|Bfl
LOC118404176|Bfl
LOC118404577|Bfl2. Run DIAMOND reciprocal best hits (RBH) between human TFs and the Bfl genes not assigned to single-copy orthogroups.
2.1 List the Bfl genes not present in orthofinderSCO_Bfl_TF_list.txt.
venn except <(grep ">" orthofinder_input/Bfl_proteins_output_final.fa) <(awk '{print $0=">"$0}' orthofinderSCO_Bfl_TF_list.txt) | sed -r 's/>(.*)/\1/g' > Bfl_not_assigned_in_SCO_gene_list.txt→ 26484 genes 2.2 Subset these genes from the FASTA and use them for the reciprocal best-hit search.
seqtk subseq orthofinder_input/Bfl_proteins_output_final.fa Bfl_not_assigned_in_SCO_gene_list.txt > Bfl_for_DIAMOND_RBH.faDIAMOND
conda activate /data/core/jordan/CondaEnv/DIAMONDdiamond makedb --in orthofinder_input/Homo_only_TF_protein_output_final.fasta -d Homo_TF
diamond makedb --in Bfl_for_DIAMOND_RBH.fa -d Bfl_for_DIAMOND_RBH
diamond blastp -q Bfl_for_DIAMOND_RBH.fa -d Homo_TF -o Bfl-vs-HomoTF.tsv --ultra-sensitive --matrix BLOSUM90 --iterate
diamond blastp -q orthofinder_input/Homo_only_TF_protein_output_final.fasta -d Bfl_for_DIAMOND_RBH -o HomoTF-vs-Bfl.tsv --ultra-sensitive --matrix BLOSUM90 --iterateSelect the best alignment using bit_score.
awk '{print $1}' Bfl-vs-HomoTF.tsv | sort -u | while read -r gene_id; do
best_hit=$(grep "^${gene_id}" Bfl-vs-HomoTF.tsv | sort -k 12n | tail -1 | awk '{print $1, $2}' OFS="\t")
echo "$best_hit"
done > Bfl_one_way_BH_for_RBH.txt
awk '{print $1}' HomoTF-vs-Bfl.tsv | sort -u | while read -r gene_id; do
best_hit=$(grep "^${gene_id}" HomoTF-vs-Bfl.tsv | sort -k 12n | tail -1 | awk '{print $2, $1}' OFS="\t") # 這邊刻意把 1st 2nd column 對調，後續比較好跟 Bfl_BH.txt 比較
echo "$best_hit"
done > Homo_one_way_BH_for_RBH.txtBfl_one_way_BH_for_RBH.txt: 2219 Bfl genes vs. 553 Homo genes Homo_one_way_BH_for_RBH.txt: 504 Bfl genes vs. 1234 Homo genes RBH
venn intersection Bfl_one_way_BH_for_RBH.txt Homo_one_way_BH_for_RBH.txt > Bfl_Homo_RBH.txtLOC118405535|Bfl ENSP00000262032.5|Homo
LOC118411066|Bfl ENSP00000360017.4|Homo
LOC118414201|Bfl ENSP00000296302.7|Homo
LOC118423077|Bfl ENSP00000290399.6|Homo
LOC118426390|Bfl ENSP00000384018.1|Homo
LOC118409569|Bfl ENSP00000377721.3|Homo
LOC118403621|Bfl ENSP00000507485.1|Homo
LOC118418362|Bfl ENSP00000326272.3|Homo
LOC118404628|Bfl ENSP00000357644.4|Homo
LOC118403563|Bfl ENSP00000257915.5|HomoThe reciprocal search identifies 351 RBH pairs.
awk '{print $1}' Bfl_Homo_RBH.txt > RBH_Bfl_TF_list.txt3. Run DIAMOND one-way best hits between human TFs and the Bfl genes still unassigned after RBH.
2.1 From Bfl_for_DIAMOND_RBH.fa, list genes not present in RBH_Bfl_TF_list.txt.
venn except <(grep ">" Bfl_for_DIAMOND_RBH.fa) <(awk '{print $0=">"$0}' RBH_Bfl_TF_list.txt) | sed -r 's/>(.*)/\1/g' > Bfl_not_assigned_in_SCO+RBH_gene_list.txt→ 26133 genes 2.2 Subset these genes from the FASTA for the one-way best-hit search.
seqtk subseq Bfl_for_DIAMOND_RBH.fa Bfl_not_assigned_in_SCO+RBH_gene_list.txt > Bfl_for_DIAMOND_one_way_BH.faDIAMOND
conda activate /data/core/jordan/CondaEnv/DIAMONDdiamond makedb --in orthofinder_input/Homo_only_TF_protein_output_final.fasta -d Homo_TF
diamond blastp -q Bfl_for_DIAMOND_one_way_BH.fa -d Homo_TF -o Bfl-vs-HomoTF_for_one_way_BH.tsv --ultra-sensitive --matrix BLOSUM90 --iterateBfl-vs-HomoTF_for_one_way_BH.txt: 1868 Bfl genes vs. 731 Homo genes Apply the filtering criteria shown below.
awk -F "\t" '$11<=1e-4' Bfl-vs-HomoTF_for_one_way_BH.tsv | awk -F "\t" '$3>=30' | awk '{print $1}' | sort -u > one_way_BH_Bfl_TF_list.txt→ 1346 Bfl genes
4. merge all 3 outputs (*Bfl_TF_list.txt)
orthofinderSCO_Bfl_TF_list.txt : 179 Bfl genes
RBH_Bfl_TF_list.txt : 351 Bfl genes
one_way_BH_Bfl_TF_list.txt : 1346 Bfl genes
cat orthofinderSCO_Bfl_TF_list.txt RBH_Bfl_TF_list.txt one_way_BH_Bfl_TF_list.txt > pred_Bfl_TF_list.txt→ 1876 predicted TFs in Bfl
Comparison with a collaborator’s Bfl TF predictions (putative_Bfl_TF_output_1.9.fa; 1,928 TFs)
venn intersection <(grep ">" ~/Shared/R145_SchneiderSQ/TO-GCN/TF-OrthoFinder/processed_protein_fasta/putative_Bfl_TF_output_1.9.fa) <(awk '{print $0=">"$0}' pred_Bfl_TF_list.txt) | wc -l- 1,372 genes overlap: 1,372 / 1,876 = 73.13% of these predictions are also in the comparison list.
- 504 genes occur only in the
orthofinderSCO_RBH_1wayBHresult. - 556 genes occur only in the collaborator’s predictions.
Tutorial 2 (GOAL: Predict the TFs in GRCh38 (human) using a TFs FASTA of Human from AnimalTFDB4.0)
Work directory
/data/core/jordan/Shared/518_projects/workshop_20231213/orthofinderSCO_RBH_1wayBH-GRCh38
conda activate /data/core/jordan/CondaEnv/orthofinder
# 為了避免搞混，先更改名稱
sed -r -i '/>/ s/\|Homo/\|GRCh38/g' orthofinder_input/GRCh38_proteins_output_final.fa 1. OrthoFinder: only get GRCh38 genes in SCO result
Use more permissive MMseqs2 settings: -s 6 -e 1e-3 -k 5 --num-iterations 5.
vim ~/CondaEnv/orthofinder/bin/scripts_of/config.json
orthofinder -f orthofinder_input/ -S mmseqs -t 150OrthoFinder assigned 17733 genes (82.6% of total) to 3658 orthogroups. Fifty percent of all genes were in orthogroups with 4 or more genes (G50 was 4) and were contained in the largest 619 orthogroups (O50 was 619). There were 1591 orthogroups with all species present and 1474 of these consisted entirely of single-copy genes.grep -wf orthofinder_input/OrthoFinder/Results_Dec08/Orthogroups/Orthogroups_SingleCopyOrthologues.txt orthofinder_input/OrthoFinder/Results_Dec08/Phylogenetic_Hierarchical_Orthogroups/N0.tsv > ./orthofinderSCO.txt- -w: exact match
awk -F "\t" '{print $4}' orthofinderSCO.txt > orthofinderSCO_GRCh38_TF_list.txt→ 1474 GRCh38 genes
2. DIAMOND RBH for GRCh38 genes not assigned to single-copy orthogroups
2.1 List the GRCh38 genes not present in orthofinderSCO_GRCh38_TF_list.txt.
venn except <(grep ">" orthofinder_input/GRCh38_proteins_output_final.fa) <(awk '{print $0=">"$0}' orthofinderSCO_GRCh38_TF_list.txt) | sed -r 's/>(.*)/\1/g' > GRCh38_not_assigned_in_SCO_gene_list.txt→ 18348 genes 2.2 Subset these genes from the FASTA and use them for the reciprocal best-hit search.
seqtk subseq orthofinder_input/GRCh38_proteins_output_final.fa GRCh38_not_assigned_in_SCO_gene_list.txt > GRCh38_for_DIAMOND_RBH.faDIAMOND
conda activate /data/core/jordan/CondaEnv/DIAMONDdiamond makedb --in orthofinder_input/Homo_only_TF_protein_output_final.fasta -d Homo_TF
diamond makedb --in GRCh38_for_DIAMOND_RBH.fa -d GRCh38_for_DIAMOND_RBH
diamond blastp -q GRCh38_for_DIAMOND_RBH.fa -d Homo_TF -o GRCh38-vs-HomoTF.tsv --ultra-sensitive --matrix BLOSUM90 --iterate
diamond blastp -q orthofinder_input/Homo_only_TF_protein_output_final.fasta -d GRCh38_for_DIAMOND_RBH -o HomoTF-vs-GRCh38.tsv --ultra-sensitive --matrix BLOSUM90 --iterateSelect the best alignment using bit_score.
awk '{print $1}' GRCh38-vs-HomoTF.tsv | sort -u | while read -r gene_id; do
best_hit=$(grep "^${gene_id}" GRCh38-vs-HomoTF.tsv | sort -k 12n | tail -1 | awk '{print $1, $2}' OFS="\t")
echo "$best_hit"
done > GRCh38_one_way_BH_for_RBH.txt
awk '{print $1}' HomoTF-vs-GRCh38.tsv | sort -u | while read -r gene_id; do
best_hit=$(grep "^${gene_id}" HomoTF-vs-GRCh38.tsv | sort -k 12n | tail -1 | awk '{print $2, $1}' OFS="\t") # 這邊刻意把 1st 2nd column 對調，後續比較好跟 Bfl_BH.txt 比較
echo "$best_hit"
done > Homo_one_way_BH_for_RBH.txtGRCh38_one_way_BH_for_RBH.txt: 1229 GRCh38 genes vs. 274 Homo genes Homo_one_way_BH_for_RBH.txt: 237 GRCh38 genes vs. 1160 Homo genes RBH
venn intersection GRCh38_one_way_BH_for_RBH.txt Homo_one_way_BH_for_RBH.txt > GRCh38_Homo_RBH.txtZNF799|GRCh38 ENSP00000411084.2|Homo
HSFX4|GRCh38 ENSP00000489814.2|Homo
LOC105379522|GRCh38 ENSP00000498738.1|Homo
FOXP1|GRCh38 ENSP00000498110.1|Homo
ZNF146|GRCh38 ENSP00000400391.1|Homo
PRDM9|GRCh38 ENSP00000296682.4|Homo
SOX7|GRCh38 ENSP00000301921.1|Homo
GTF2IRD2|GRCh38 ENSP00000498563.1|Homo
PLSCR2|GRCh38 ENSP00000345494.4|Homo
SMARCA4|GRCh38 ENSP00000296302.7|HomoThe reciprocal search identifies 201 RBH pairs.
awk '{print $1}' GRCh38_Homo_RBH.txt > RBH_GRCh38_TF_list.txt3. DIAMOND one-way best hits for GRCh38 genes still unassigned after RBH
2.1 From GRCh38_for_DIAMOND_RBH.fa, list genes not present in RBH_GRCh38_TF_list.txt.
venn except <(grep ">" GRCh38_for_DIAMOND_RBH.fa) <(awk '{print $0=">"$0}' RBH_GRCh38_TF_list.txt) | sed -r 's/>(.*)/\1/g' > GRCh38_not_assigned_in_SCO+RBH_gene_list.txt→ 18147 genes 2.2 Subset these genes from the FASTA for the one-way best-hit search.
seqtk subseq GRCh38_for_DIAMOND_RBH.fa GRCh38_not_assigned_in_SCO+RBH_gene_list.txt > GRCh38_for_DIAMOND_one_way_BH.faDIAMOND
conda activate /data/core/jordan/CondaEnv/DIAMOND#diamond makedb --in orthofinder_input/Homo_only_TF_protein_output_final.fasta -d Homo_TF
diamond blastp -q GRCh38_for_DIAMOND_one_way_BH.fa -d Homo_TF -o GRCh38-vs-HomoTF_for_one_way_BH.tsv --ultra-sensitive --matrix BLOSUM90 --iterateGRCh38-vs-HomoTF_for_one_way_BH.tsv: 1028 GRCh38 genes vs. 268 Homo genes Apply the filtering criteria shown below.
awk -F "\t" '$11<=1e-4' GRCh38-vs-HomoTF_for_one_way_BH.tsv | awk -F "\t" '$3>=30' | awk '{print $1}' | sort -u > one_way_BH_GRCh38_TF_list.txt→ 750 GRCh38 genes (the original line used “Bfl” in this human example).
4. merge all 3 outputs (*GRCh38_TF_list.txt)
orthofinderSCO_GRCh38_TF_list.txt : 1474 GRCh38 genes
RBH_GRCh38_TF_list.txt : 201 GRCh38 genes
one_way_BH_GRCh38_TF_list.txt : 750 GRCh38 genes
cat orthofinderSCO_Bfl_TF_list.txt RBH_Bfl_TF_list.txt one_way_BH_Bfl_TF_list.txt > pred_GRCh38_TF_list.txt→ 2425 predicted TFs in GRCh38
Comparison with the human TF list from AnimalTFDB 4.0
venn intersection <(sed -r '/>/ s/\|Homo/\|GRCh38/g' ~/Shared/R145_SchneiderSQ/TO-GCN/TF-AnimalTFDB4.0/Homo_only_TF_protein_output_final_symbol.fasta) <(awk '{print $0=">"$0 }' pred_GRCh38_TF_list.txt) | wc -l- 1,613 genes overlap: 1,613 / 2,425 = 66.52% of the predictions are in the reference TF list.
- The recorded comparison lists 811 genes only in the
orthofinderSCO_RBH_1wayBHresult. This count differs by one from 2,425 − 1,613 = 812; retain that discrepancy when interpreting the record. - 31 genes occur only in AnimalTFDB 4.0. Using the AnimalTFDB 4.0 human TF list as a reference, the proportions of predicted genes present in that list are:
- AnimalTFDB4.0-1wayBH Intersection: 92.73%
- AnimalTFDB4.0-1wayBH Union: 64.88%
- orthofinder-RBH-1wayBH: 66.52%