Setting a batch analysis for single alignment-wide ω

[desmodus1984]

Hi,
I have identified orthologues with Orthofinder, and I am interested in identifying genes under selection with hyphy and if possible doing the analysis in parallel.

I installed HyPhy with conda, and I am trying to understand how to set the parameters for the analysis.
Since I have duplications and potential paralogues I will restrict the analysis to alignments with single-copy orthologues.

I was reading the tutorial on the website and it says: "This tutorial assumes you are specifically using the HyPhy executable HYPHYMP. If you have installed a different executable (e.g. HYPHYMPI), you may need to alter some commands."
Since I installed with conda, I'd appreciate some help with the code.

The website indicates to use:

# Path to files of interest
FILEPATH=/absolute/path/to/directory/with/all/your/files/

# Run BUSTED on all files $FILEPATH/*.fas
for FILE in `ls $FILEPATH/*.fas`; do
    echo "Running BUSTED on $FILE"
    echo `(echo "10"; echo "4"; echo "1"; echo $FILE; echo "Y"; echo "1"; echo "d") | HYPHYMP`
done

I don't have HYPHYMP. First, I want to do a quick single alignment-wide ω to find genes with signal of selection and then do BUSTED with one or two species of interest, and since I have access to a cluster I would like to submit several jobs with HYPHYMPI.

The code I thought would be:

#!/bin/bash

#PBS -l nodes=1:ppn=32
#PBS -l walltime=03:00:00
#PBS -l mem=100gb
#PBS -t 1-100
#PBS -k oe
#PBS -N single-msa-w
#PBS -o ${PBS_JOBNAME}_${PBS_JOBID}.log
#PBS -e ${PBS_JOBNAME}_${PBS_JOBID}.err
#PBS -l naccesspolicy=singlejob

# Load bash environment
source ~/.bashrc

# Init conda
source /data/common/MMA/mamba.init.sh
conda activate NGS-forge

cd /data/common/MMA/data/genomes/

# Save starting time
start_time=$(date +%s)
echo "Start ${PBS_JOBNAME}: Job started: $(date)"

# Get the sample name for this array job
base=$(sed -n "${PBS_ARRAYID}p" dn.dn.sps.txt)

#Calculate single alignment-wide ω

echo "Starting single alignment-wide ω on ${base} at $(date)"

HYPHYMPI CPU=32  ...

I can't find the code for running the single alignment-wide ω.
Could you please indicate the code to run the single alignment-wide ω?

Then, to run busted on one branch out of 11 (species), so I assume it would be:

HYPHYMPI CPU=32  busted-ph --alignment path/to/${base}.fasta --tree my_tree.nwk --branches Hydlep --output ${base}.json

As mentioned before, I ran Orthofinder to identify orthologues and orthogroups. The issue/ perhaps advantage, is that it adds the file name as prefix to the protein fasta headeres, so the multiple sequence alignments look like:
Ailmel_rna-XM_034655643.1
Calurs_rna-XM_073887989.1
Canlup_rna-XM_038532347.1
Hydlep_GNX-026372-mRNA
Lepwed_LOC_00017029-mRNA-10
Mirang_rna-XM_064596364.1
Mirleo_rna-XM_035014847.1
Musnig_rna-XM_059416383.1
Neosch_rna-XM_044916755.1
Ursarc_rna-XM_026498632.3

In my case, the tree is:
((((Calurs,(Halgry,((((Hydlep,Lepwed),(Mirang,Mirleo))),Neosch))),Musnig),(Ursarc,Ailmel)),Canlup);

Can I reconcile the alignment and the tree with trim-tree, next either parse or save the reconciled files and then run the analysis?

Orthofinder uses protein data. I looked at the alignments and many seem to need trimming due to issues like gaps. Is there any program that you recommend for doing protein-alignment trimming? I am planning to use trimAI but I wanted to ask you first.

Lastly, after trimming the aa, can HyPhy run the single alignment-wide ω/ busted analysis with the protein alignment (set input data) or must I backtranslate to codon sequences?

Thank you.

[spond]

Dear @desmodus1984,

Sorry for not responding sooner.

1. HyPhy requires codon data, so you MUST backtranslate. dN/dS estimation requires both synonymous and non-synonymous signal, and with amino-acid data you obviously lose the former.
2. trim-tree does not perform partial name matching. I think your best bet is to run a little Python script (requires Biopython); attached

reconcile.py.txt

Best,
Sergei

PS One small remark: busted-ph does not use MPI at all, so you should run the multi-threaded hyphy CPU=N, where N is the number of threads.

[desmodus1984]

Dear @spond
Thank you very much for the response. I would appreciate if you don't mind asking a few more questions.

I figure out a way to take advantage of the output format of Orthofinder and make the alignment and tree compatible with a very easy and quick fix. Since Orthofinder adds the prefix of the protein files, I just have to add the species prefix to the CDS to then retrieve them for each orthogroup.

Here is my first issue. I am analyzing marine mammals and I added terrestrial species: dog, panda, and brown bear. When I tried backtranslating with pal2nal, weirdly, annoyingly, and funnily, it had an error with the panda sequences. I extracted the CDS and proteins with gffread, I don't think the issue is with gffread but with the panda genome/annotation. I am sure that I will encounter errors with other species.
I aligned the CDS sequences with MACSE, and the panda sequence didn't look that bad.

The species with the long stretch is the dog!
I thought that since I have the CDS alignment I could use it to move forward.
I read that it is recommended to trim the alignments. I would like to know which software do you recommend for trimming, or do you recommend other program for doing the backtranslation?
I was thinking of using ClipKit, but you might know a better one for trimming alignments.

Sorry for being dumb, I want to calculate single alignment-wide ω, but I couldn't find which command option for doing this. Would you please indicate it. Thank you.
I have searched for the command/option but the only page mentioning it uses CLI site and not command line.

My idea is to first calculate ω for all single-copy orthogroups, and then on the ones with selection signal do RELAX, BUSTED, and BUSTED-PH- analysis. Does this workflow look reasonable for you?

Do you have any recommendations for how distant species should be, and how many (minimal sample size) are enough for calculating ω?
I work with two sister species with different diets and different diving capacities (major interest of my PI). I expect that some genes will have signal of purifying selection and other of relaxed selection. I have them, and I added in ladder-wise manner, one species from upper clades following this figure 1 (tribe>subfamily>family>Superfamily) until reaching dog, total: 11 species.

Thank you very much for the help.