Example (Hofstenia miamia):

REQUIRE:

  • Hof Peptide FASTA: 22632 genes
  • Hof Nucleotide FASTA: 20966 genes

1.1 Convert the peptide FASTA line endings with dos2unix if needed

sudo dos2unix Isodiametra_pulchra/Isodiametra_aa.fa

1.2 Modify the header

awk '{ if (/^>/) { $0=$0"|Hof" } print }' Hofstenia_miamia/Hofstenia_aa.faa > Hofstenia_miamia/processed/Hof_output_1.2.faa

1.3 MMseqs2 RBH

Match the peptide and nucleotide FASTA identifiers using a translated search (TBLASTN).

mmseqs createdb Hofstenia_miamia/processed/Hof_output_1.2.faa Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_pep;
mmseqs createdb Hofstenia_miamia/Hof.fasta Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_nucl

TBLASTN

mmseqs rbh \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_pep \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_nucl \
Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl \
Hofstenia_miamia/processed/mmseqs2_RBH/tmp --threads 100 -e 1e-5 -s 7;

mmseqs convertalis \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_pep \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_nucl \
Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl \
Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_result

awk -F"\t" '$3 >= 0.8' Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_result > Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.3

→ 12723 sequence pairs ⇒ *filter out alignment identity<0.8 *⇒ 12033 sequence pairs

  • 11852 uniq nucleotide ID
  • 12022 uniq peptide ID 1.4 Remove matches with repeated peptide IDs.
grep -v -f <(awk '{print $1}' Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.3 | sort | uniq -c | awk '$1 > 1 {print $2}') \
Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.3 > Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.4

→ 12011 gene pairs

1.5 Remove the 12,011 matched peptide IDs from the peptide FASTA, then run a one-way best-hit search on the remaining sequences

seqtk seq Hofstenia_miamia/processed/Hof_output_1.2.faa | paste - - | \
grep -v -f <(awk -F"\t" '{print $1}' Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.4 | sort -u) | \
tr "\t" "\n" > Hofstenia_miamia/processed/mmseqs2_RBH/Hof_output_1.5.faa
  • seqtk seq places each sequence on one line.
  • paste - - combines a header and sequence into one tab-delimited line.
  • tr "\t" "\n" converts the tabs back into line breaks. → 10621 unnamed peptide IDs

1.6 1way BH (peptide to nucleotide)

mmseqs createdb Hofstenia_miamia/processed/mmseqs2_RBH/Hof_output_1.5.faa Hofstenia_miamia/processed/mmseqs2_BH/HofDB_output_1.6;

mmseqs search \
Hofstenia_miamia/processed/mmseqs2_BH/HofDB_output_1.6 \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_nucl \
Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_pep.to.nucl \
Hofstenia_miamia/processed/mmseqs2_BH/tmp --threads 100 -e 1e-5 -s 7;


mmseqs convertalis \
Hofstenia_miamia/processed/mmseqs2_BH/HofDB_output_1.6 \
Hofstenia_miamia/processed/mmseqs2_RBH/HofDB_nucl \
Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_pep.to.nucl \
Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_result


awk -F"\t" '$3 >= 0.8' Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_result > Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_output_1.6

→ 5456 gene pairs (3517 unique peptide IDs)

1.7 For each peptide ID we only keep the best gene pair

  • highest bitscore and identity
awk -F"\t" '{print $1}' Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_output_1.6 | sort -u | while read -r pep_id; do

desire_pep_nucl_pair=$(grep -a "^${pep_id}" Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_output_1.6 | sort -k 12n -k 3n | tail -1 | awk '{print}' OFS="\t" )

echo "$desire_pep_nucl_pair"

done > Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_output_1.7

→ 3517 gene pairs

1.8 Concatenate the 3,517 additional gene pairs with the previous 12,011 pairs and check for repeated peptide IDs

cat Hofstenia_miamia/processed/mmseqs2_BH/Hof_BH_output_1.7 Hofstenia_miamia/processed/mmseqs2_RBH/Hof_RBH_pep.vs.nucl_output_1.4 | awk '{print $1, $2}' OFS="\t" > Hofstenia_miamia/processed/Hof_pep.vs.nucl_output_1.8

→ 15528 gene pairs

  • 15528 uniq peptide ID
  • 12701 uniq nucleotide ID

1.9 Rename the peptide FASTA identifiers using this mapping

Keep the original peptide ID when no corresponding nucleotide ID is available.

awk -F ' ' 'NR==FNR{ map[">"$1] = ">"$2; next;} {if (/^>/) {$0=(map[$0]) ? map[$0] : $0} print} ' Hofstenia_miamia/processed/Hof_pep.vs.nucl_output_1.8 \
Hofstenia_miamia/processed/Hof_output_1.2.faa > Hofstenia_miamia/processed/Hof_output_1.9.faa

The awk conditional expression uses this form: (condition) ? (expression_if_true) : (expression_if_false)

1.10 Rename duplicated gene IDs to make all sequence identifiers unique, as in the extract_isoform workflow

Append _1, _2, and so on for duplicates; append _1 even when an ID occurs once. If the original header already contains _ or |, use another temporary separator, such as .1, .2, or .3.

awk '{ if (NR==FNR){ 
			if (/^>/) { s[$0]++ }; next;} 
			if (/^>/) {k[$0]++; 
								 if (s[$0]>1){$0=$0"_"k[$0]} else {$0=$0"_1"}} 
print}' Hofstenia_miamia/processed/Hof_output_1.9.faa Hofstenia_miamia/processed/Hof_output_1.9.faa > Hofstenia_miamia/processed/Hof_output_1.10.faa

1.11 Output sequence length (with new unique gene IDs)

seqkit fx2tab --length --name --header-line Hofstenia_miamia/processed/Hof_output_1.10.faa > Hofstenia_miamia/processed/Hof_output_1.11.txt

1.12 Identify the longest isoform for each gene. Inputs:

  • Output from step 1.9
  • Output from step 1.11
grep -a ">" Hofstenia_miamia/processed/Hof_output_1.9.faa  | sed 's/>//' | sort -u | while read -r gene_id; do

longest=$(grep -a "^${gene_id}_" Hofstenia_miamia/processed/Hof_output_1.11.txt | sort -k 2n | tail -1 | awk '{print $1}')

echo "$longest"

done > Hofstenia_miamia/processed/Hof_output_1.12.txt

→ 19805 gene left 1.13 Check if there are still duplicates

cat Hofstenia_miamia/processed/Hof_output_1.12.txt | awk '{a[$0]++} {print a[$0],$0} ' | awk '{if ( $1+0 > 1 ) print $1,$2}'

No duplicate IDs remain. 1.14 Use this list to subset the output from step 1.10.

seqtk subseq Hofstenia_miamia/processed/Hof_output_1.10.faa Hofstenia_miamia/processed/Hof_output_1.12.txt > Hofstenia_miamia/processed/Hof_output_1.14.fa

→ 19805 genes 1.15 Remove the temporary underscore suffix from the headers.

sed -r 's/_.*$//' Hofstenia_miamia/processed/Hof_output_1.14.fa > Hofstenia_miamia/processed/Hof_output_1.15.fa
← Back to projects