extract_isoform prepares representative protein sequences from FASTA and GFF/GTF files for comparative analysis.
Goal:
FASTA and annotation files may contain multiple isoforms for a gene. → extract_isoform pipelines help us to:
- Remove isoforms sharing the same gene ID
- Remove redundant isoforms carrying different gene IDs
Inputs:
- GFF or GTF files
- Protein FASTA file
Workflow:
- Remove isoforms sharing the same gene ID
- 1.0 Simplify the FASTA header
- 1.1 Create a header-ID-to-gene-ID mapping so that FASTA headers can be replaced with gene IDs.
- 1.2 Replace the header ID in fasta to Gene ID based on this list
- 1.3 Report duplicated gene IDs and their counts; skip to step 1.9 if there are no duplicates.
- 1.4 Rename duplicated gene IDs so that every sequence identifier is unique.
- 1.5 Output sequence length (with new unique gene IDs)
- 1.6 Identify the longest isoform for each gene.
- 1.7 Check if there are still duplicates
- 1.8 Use that list to subset the output from step 1.4.
- 1.9 Remove the temporary suffix after the underscore and append a species code to each header.
- AGAT process the GFF/GTF: Remove redundant isoforms carrying different gene IDs
- 2.1 Merge
- 2.2 Keep only the longest isoform
- 2.3 Extract Gene IDs (modify headers and deduplicate)
- Use the list from step 2.3 to subset the FASTA produced in step 1.9.
Example: Amphioxus (Bfl):
1. Remove isoforms sharing the same gene ID
see gff:
NC_049979.1 RefSeq region 1 35335272 . + . ID=NC_049979.1:1..35335272;Dbxref=taxon:7739;Name=1;chromosome=1;country=USA: Old Tampa Bay%2C Florida;gbkey=Src;genome=chromosome;mol_type=genomic DNA;sex=male;strain=S238N-H82;tissue-type=testes
NC_049979.1 Gnomon gene 209 23406 . - . ID=gene-LOC118414000;Dbxref=GeneID:118414000;Name=LOC118414000;end_range=23406,.;gbkey=Gene;gene=LOC118414000;gene_biotype=protein_coding;partial=true
NC_049979.1 Gnomon mRNA 209 23406 . - . ID=rna-XM_035817715.1;Parent=gene-LOC118414000;Dbxref=GeneID:118414000,Genbank:XM_035817715.1;Name=XM_035817715.1;Note=The sequence of the model RefSeq transcript was modified relative to this genomic sequence to represent the inferred CDS: added 630 bases not found in genome assembly;end_range=23406,.;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=LOC118414000;inference=similar to RNA sequence (same species):INSD:GETA01027940.1;model_evidence=Supporting evidence includes similarity to: 1 EST%2C 6 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 75 samples with support for all annotated introns;partial=true;product=serine/threonine-protein phosphatase 4 regulatory subunit 4-like%2C transcript variant X2;transcript_id=XM_035817715.1
NC_049979.1 Gnomon exon 23333 23406 . - . ID=exon-XM_035817715.1-1;Parent=rna-XM_035817715.1;Dbxref=GeneID:118414000,Genbank:XM_035817715.1;Note=The sequence of the model RefSeq transcript was modified relative to this genomic sequence to represent the inferred CDS: added 630 bases not found in genome assembly;end_range=23406,.;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=LOC118414000;inference=similar to RNA sequence (same species):INSD:GETA01027940.1;partial=true;product=serine/threonine-protein phosphatase 4 regulatory subunit 4-like%2C transcript variant X2;transcript_id=XM_035817715.1
NC_049979.1 Gnomon exon 22755 22861 . - . ID=exon-XM_035817715.1-2;Parent=rna-XM_035817715.1;Dbxref=GeneID:118414000,Genbank:XM_035817715.1;Note=The sequence of the model RefSeq transcript was modified relative to this genomic sequence to represent the inferred CDS: added 630 bases not found in genome assembly;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=LOC118414000;inference=similar to RNA sequence (same species):INSD:GETA01027940.1;partial=true;product=serine/threonine-protein phosphatase 4 regulatory subunit 4-like%2C transcript variant X2;transcript_id=XM_035817715.1
NC_049979.1 Gnomon exon 22190 22297 . - . ID=exon-XM_035817715.1-3;Parent=rna-XM_035817715.1;Dbxref=GeneID:118414000,Genbank:XM_035817715.1;Note=The sequence of the model RefSeq transcript was modified relative to this genomic sequence to represent the inferred CDS: added 630 bases not found in genome assembly;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=LOC118414000;inference=similar to RNA sequence (same species):INSD:GETA01027940.1;partial=true;product=serine/threonine-protein phosphatase 4 regulatory subunit 4-like%2C transcript variant X2;transcript_id=XM_035817715.1
NC_049979.1 Gnomon exon 21601 21722 . - . ID=exon-XM_035817715.1-4;Parent=rna-XM_035817715.1;Dbxref=GeneID:118414000,Genbank:XM_035817715.1;Note=The sequence of the model RefSeq transcript was modified relative to this genomic sequence to represent the inferred CDS: added 630 bases not found in genome assembly;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=LOC118414000;inference=similar to RNA sequence (same species):INSD:GETA01027940.1;partial=true;product=serine/threonine-protein phosphatase 4 regulatory subunit 4-like%2C transcript variant X2;transcript_id=XM_035817715.1see protein fasta:
>NP_007757.1 cytochrome c oxidase subunit I (mitochondrion) [Branchiostoma floridae]
MYITRWLFSTNHKDIGTLYFIFGAWAAMVGTAMSLLIRAELSQPGALLGDDHLYNVIVTAHAFVMIFFMVMPIMIGGFGN
WLVPMMIGAPDMAFPRMNNMSFWMLPPSFSLLLASSAVEAGVGTGWTVYPPLSSNIAHAGASVDLAIFSLHLAGVSSILG
AINFITTIHNMRASIEWNRVPLFVWSIWVTAYLLLLSLPVLAGAITMLLTDRNINTTFFDPSGGGDPILYEHLFWFFGHP
EVYILILPGFGIISHIIIHYAGKLRSFGYLGMTWAMFTIGLLGFWVWAHHMFTVGMDVDTRSYFTAATMVIAVPTGIKVF
SWLATLSGSKQLKWETPLLWALGFIFLFTVGGLTGIVLANSSLDIVLHDTYYVVAHFHYVLSMGAVFSILAGLTYWFPLF
TGFTLQESWTKIHFFVMFVGVNLTFFPQHFLGLAGMPRRYSDYPDAYTIWNVISSLGSIISLGSVVFFLFILWEAFSSQR
KAVPASHASHSAEWLMGCPPLFHTHEELPFISKKY
>NP_007758.1 cytochrome c oxidase subunit II (mitochondrion) [Branchiostoma floridae]
MATPAQLGLMDAASPVMEEMIYFHDHVMLVLILITCLIFYSMLVLISSKYIYRFLTDGHVIETVWTVIPAIILVVVALPS→ using grep in gff, we can see that only CDS sequence ($3 == “CDS”) has FASTA ID (“Name=”) and Gene ID (“gene=”)
1.0 Simplify the FASTA header
# ponytail: the POSIX collating element names a literal space and survives Notion Markdown.
sed -r 's/(>[^[. .]]*).*/\1/' processed_protein_fasta/Bfl_proteins.fa > processed_protein_fasta/Bfl_proteins_output_1.0.fa1.1 Create the header-ID-to-gene-ID mapping
- Columns here are tab-delimited, rather than separated by ordinary spaces.
awk '!/^#/'removes the initial comment rows that begin with#.
# subset only the lines with CDS in the 3rd column
awk -F '\t' '$3 == "CDS"' processed_protein_fasta/Bfl.gff | awk '!/^#/' |
# subset the Last column and use ";" as delim, and subset patter "Name=" + "gene="
awk -F '\t' '{split($9,a,";"); name=""; gene=""; for (i in a) {if (a[i] ~ /^Name=/){name=a[i]} else if (a[i] ~ /^gene=/){gene=a[i]}} print name, gene}' |
# 去除 Name= 與 gene=
sed -r 's/Name=//g' | sed -r 's/gene=//g' > processed_protein_fasta/Bfl_output_1.1.txtExample output: repeated rows are expected and do not affect this step
XP_035673608.1 LOC118414000
XP_035673608.1 LOC1184140001.2 Replace the header ID in fasta to Gene ID based on this list
awk -F ' ' 'NR==FNR{ map[">"$1] = ">"$2; next;} {if (/^>/) {$0=map[$0]} print} ' processed_protein_fasta/Bfl_output_1.1.txt processed_protein_fasta/Bfl_proteins_output_1.0.fa > processed_protein_fasta/Bfl_proteins_output_1.2.faExample output: several sequences now share the same header
>COX1
MYITRWLFSTNHKDIGTLYFIFGAWAAMVGTAMSLLIRAELSQPGALLGDDHLYNVIVTAHAFVMIFFMVMPIMIGGFGN
WLVPMMIGAPDMAFPRMNNMSFWMLPPSFSLLLASSAVEAGVGTGWTVYPPLSSNIAHAGASVDLAIFSLHLAGVSSILG
AINFITTIHNMRASIEWNRVPLFVWSIWVTAYLLLLSLPVLAGAITMLLTDRNINTTFFDPSGGGDPILYEHLFWFFGHP
EVYILILPGFGIISHIIIHYAGKLRSFGYLGMTWAMFTIGLLGFWVWAHHMFTVGMDVDTRSYFTAATMVIAVPTGIKVF
SWLATLSGSKQLKWETPLLWALGFIFLFTVGGLTGIVLANSSLDIVLHDTYYVVAHFHYVLSMGAVFSILAGLTYWFPLF
TGFTLQESWTKIHFFVMFVGVNLTFFPQHFLGLAGMPRRYSDYPDAYTIWNVISSLGSIISLGSVVFFLFILWEAFSSQR
KAVPASHASHSAEWLMGCPPLFHTHEELPFISKKY
>COX2
MATPAQLGLMDAASPVMEEMIYFHDHVMLVLILITCLIFYSMLVLISSKYIYRFLTDGHVIETVWTVIPAIILVVVALPS1.3 Report duplicate gene IDs and counts; skip to step 1.9 if none are present
grep ">" processed_protein_fasta/Bfl_proteins_output_1.2.fa | awk '{a[$0]++} {print a[$0],$0} ' | awk '{if ( $1+0 > 1 ) print $1,$2}' | sort > processed_protein_fasta/Bfl_proteins_output_1.3.txtExample output
2 >LOC118403012
2 >LOC118403089
3 >LOC1184030891.4 Rename duplicate gene IDs to make every sequence identifier unique
Append _1, _2, and so on for duplicates; append _1 even when an ID occurs once.
If the original header already contains underscores, choose another temporary delimiter, such as |1, |2, or |3, to distinguish isoforms.
awk '{ if (NR==FNR){
if (/^>/) { s[$0]++ }; next;}
if (/^>/) {k[$0]++;
if (s[$0]>1){$0=$0"_"k[$0]} else {$0=$0"_1"}}
print}' processed_protein_fasta/Bfl_proteins_output_1.2.fa processed_protein_fasta/Bfl_proteins_output_1.2.fa > processed_protein_fasta/Bfl_proteins_output_1.4.faExample output
>COX1_1
MYITRWLFSTNHKDIGTLYFIFGAWAAMVGTAMSLLIRAELSQPGALLGDDHLYNVIVTAHAFVMIFFMVMPIMIGGFGN
WLVPMMIGAPDMAFPRMNNMSFWMLPPSFSLLLASSAVEAGVGTGWTVYPPLSSNIAHAGASVDLAIFSLHLAGVSSILG
AINFITTIHNMRASIEWNRVPLFVWSIWVTAYLLLLSLPVLAGAITMLLTDRNINTTFFDPSGGGDPILYEHLFWFFGHP
EVYILILPGFGIISHIIIHYAGKLRSFGYLGMTWAMFTIGLLGFWVWAHHMFTVGMDVDTRSYFTAATMVIAVPTGIKVF
SWLATLSGSKQLKWETPLLWALGFIFLFTVGGLTGIVLANSSLDIVLHDTYYVVAHFHYVLSMGAVFSILAGLTYWFPLF
TGFTLQESWTKIHFFVMFVGVNLTFFPQHFLGLAGMPRRYSDYPDAYTIWNVISSLGSIISLGSVVFFLFILWEAFSSQR
KAVPASHASHSAEWLMGCPPLFHTHEELPFISKKY
>COX2_1
MATPAQLGLMDAASPVMEEMIYFHDHVMLVLILITCLIFYSMLVLISSKYIYRFLTDGHVIETVWTVIPAIILVVVALPS1.5 Output sequence length (with new unique gene IDs)
seqkit fx2tab --length --name --header-line processed_protein_fasta/Bfl_proteins_output_1.4.fa > processed_protein_fasta/Bfl_proteins_output_1.5.txtExample output
#name length
COX1_1 515
COX2_1 230
ATP8_1 54
ATP6_1 227
COX3_1 262
ND3_1 117
ND4L_1 91
ND4_1 452
ND5_1 5991.6 Identify the longest isoform for each gene
Inputs:
- Output from step 1.2
- Output from step 1.5
grep ">" processed_protein_fasta/Bfl_proteins_output_1.2.fa | sed 's/>//' | sort -u | while read -r gene_id; do
longest=$(grep "^${gene_id}_" processed_protein_fasta/Bfl_proteins_output_1.5.txt | sort -k 2n | tail -1 | awk '{print $1}')
echo "$longest"
done > processed_protein_fasta/Bfl_proteins_output_1.6.txtExample output
ATP6_1
ATP8_1
COX1_1
COX2_1
COX3_1
CYTB_1
LOC118403006_1
LOC118403008_1
LOC118403009_1
LOC118403010_11.7 Check if there are still duplicates
cat processed_protein_fasta/Bfl_proteins_output_1.6.txt | awk '{a[$0]++} {print a[$0],$0} ' | awk '{if ( $1+0 > 1 ) print $1,$2}'No duplicate IDs remain.
1.8 Use the selected-ID list to subset the output from step 1.4
seqtk subseq processed_protein_fasta/Bfl_proteins_output_1.4.fa processed_protein_fasta/Bfl_proteins_output_1.6.txt > processed_protein_fasta/Bfl_proteins_output_1.8.fa1.9 Remove the temporary underscore suffix and append |Bfl
sed -r 's/(>[^_]*).*/\1/' processed_protein_fasta/Bfl_proteins_output_1.8.fa | awk '{ if (/^>/) { $0=$0"|Bfl" } print }' > processed_protein_fasta/Bfl_proteins_output_1.9.faExample output
>COX1|Bfl
MYITRWLFSTNHKDIGTLYFIFGAWAAMVGTAMSLLIRAELSQPGALLGDDHLYNVIVTAHAFVMIFFMVMPIMIGGFGNWLVPMMIGAPDMAFPRMNNMSFWMLPPSFSLLLASSAVEAGVGTGWTVYPPLSSNIAHAGASVDLAIFSLHLAGVSSILGAINFITTIHNMRASIEWNRVPLFVWSIWVTAYLLLLSLPVLAGAITMLLTDRNINTTFFDPSGGGDPILYEHLFWFFGHPEVYILILPGFGIISHIIIHYAGKLRSFGYLGMTWAMFTIGLLGFWVWAHHMFTVGMDVDTRSYFTAATMVIAVPTGIKVFSWLATLSGSKQLKWETPLLWALGFIFLFTVGGLTGIVLANSSLDIVLHDTYYVVAHFHYVLSMGAVFSILAGLTYWFPLFTGFTLQESWTKIHFFVMFVGVNLTFFPQHFLGLAGMPRRYSDYPDAYTIWNVISSLGSIISLGSVVFFLFILWEAFSSQRKAVPASHASHSAEWLMGCPPLFHTHEELPFISKKY
>COX2|Bfl
MATPAQLGLMDAASPVMEEMIYFHDHVMLVLILITCLIFYSMLVLISSKYIYRFLTDGHVIETVWTVIPAIILVVVALPSLKLLYLTDELDNPQLTIKSVGHQWYWSYEYTDYYDIEFDSYMLPLGDLSKGDARLLEVDNRVVLPVDTSVRVLVTAADVIHSWTVPSLGLKMDAVPGRLNQLALQCSRVGTFYGQCSEICGANHSFMPIVIEAVPVEVFEGWCDMMLDEE
>ATP8|Bfl
MPQLNPIPWVFLFFLVWLVLGFLGLQKFTSIVTTTLDDSSEEVEVKSKEYSWPW
>ATP6|Bfl
MMVSLFSQFDSPWLLNIPLVLLALIMPWKLFVSFGPSWAGTRSSRLVYATMETLMSQVMQPLNKLGFRWVVLFSSLMLMLMTLNVIGLFPYTFTPTTQLSMNLGLAVPLWLGTVVYGFRNHPVIALAHLCPEGAPNLLVPVLVVVETLSILMRPLALGLRLTANLTAGHLLMHLISSAVLGLMELSVMLSGITLLLLVFLTMLEIAVALIQGYVFAILVTLYLDENL
>COX3|Bfl
MTGYQPHPWHLVEPSPWPLVGGSAAFTLTVGLVMWFHYNSISLMILGLVMIVATMIQWWRDVIREATFQGCHTSYVLSGLRRGMVLFIVSEVFFFLAFFWAFFHSSLAPTVELGVTWPPVGVHPLNAFAVPLLNTAVLLSSGVTVTWAHHALMEGKRTEAIQSLAITVMLGLYFTGLQAWEYYEAPFTIADSVYGSTFFVATGFHGLHVIIGSTFLMVCLGRQVFYHYTSSHHFGFEAAAWYWHFVDVVWLFLYVCIYWWGS2. AGAT process the GFF/GTF: Remove redundant isoforms carrying different gene IDs
Display the AGAT configuration file, agat_config.yaml.
agat config --exposeSet the following option in agat_config.yaml to true.
# should overlapping loci (at CDS level) be merged in a single locus
# Only one gene is kept, and the mRNA features become isoforms.
merge_loci: true2.1 Merge
agat_convert_sp_gxf2gxf.pl --gff processed_protein_fasta/Bfl.gff -o processed_protein_fasta/Bfl_output_2.1.gff2.2 Keep only the longest isoform
agat_sp_keep_longest_isoform.pl --gff processed_protein_fasta/Bfl_output_2.1.gff -o processed_protein_fasta/Bfl_output_2.2.gff2.3 Extract Gene IDs (modify headers and deduplicate)
Check the sequence type to extract, the attribute containing the gene ID, and the species code to append to the header.
awk '$3 == "CDS"' processed_protein_fasta/Bfl_output_2.2.gff | awk -F '\t' '{split($9,a,";"); gene=""; for (i in a) {if (a[i] ~ /^gene=/){gene=a[i]}} print gene}' | sed -r 's/gene=//g' | sort -u |awk '{$0=$0"|Bfl"; print }' > processed_protein_fasta/Bfl_output_2.3.txtExample output
ATP6|Pdu
COX1|Pdu
COX2|Pdu
COX3|Pdu
CYTB|Pdu
LOC118403006|Pdu
LOC118403008|Pdu
LOC118403009|Pdu
LOC118403010|Pdu
LOC118403011|Pdu3. Use the list from step 2.3 to subset the FASTA produced in step 1.9
seqtk subseq processed_protein_fasta/Bfl_proteins_output_1.9.fa processed_protein_fasta/Bfl_output_2.3.txt > processed_protein_fasta/Bfl_proteins_output_final.fa43041 → 26663 genes left