SheepwormJM
commented
1 year ago Hi,
I too have found a similar thing. I haven't been using a gff file for the primary_transcripts.py though.
My code:
for f in *fa ; do python /Orthofinder/OrthoFinder/tools/primary_transcript.py $f ; done
For one species it worked fine. For this species, each protein sequence had the following (or very similar) header:
>HCON_00022890-00001 transcript=HCON_00022890-00001 gene=HCON_00022890
Secondary (or more) transcripts were denoted by -00002 etc. But otherwise all headers are the same.
It selected the same number of genes as I would expect to have primary transcripts.
For another species, C. elegans, it didn't work as well. It retained 25803 genes out of 28447 'accesssions', while I found only 19985 genes. The headers for this file were different, and I think they are causing an issue. I have a feeling (only really looked at the first page), that those with the header the same as the first two lines are processed correctly, and the rest may not be.
Initial file: `
2L52.1a wormpep=CE32090 gene=WBGene00007063 status=Confirmed uniprot=A4F336 insdc=CCD61130.1 2L52.1b wormpep=CE50569 gene=WBGene00007063 status=Confirmed uniprot=A0A0K3AWR5 insdc=CTQ86426.1 2RSSE.1a wormpep=CE32785 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=A4F337 insdc=CCD61138.1 product="Rho-GAP domain-containing protein" 2RSSE.1b wormpep=CE48524 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=U4PAK0 insdc=CDH92922.1 product="Rho-GAP domain-containing protein" 2RSSE.1c wormpep=CE52653 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=A0A2X0T1Z3 insdc=SPS41576.1 product="Rho-GAP domain-containing protein" 3R5.1a wormpep=CE24758 gene=WBGene00007065 locus=pot-3 status=Confirmed uniprot=G5EFG7 insdc=CAA21777.2 product="POT1PC domain-containing protein" `
Primary transcript file: `
WBGene00007063 2RSSE.1a wormpep=CE32785 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=A4F337 insdc=CCD61138.1 product="Rho-GAP domain-containing protein" 2RSSE.1b wormpep=CE48524 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=U4PAK0 insdc=CDH92922.1 product="Rho-GAP domain-containing protein" 2RSSE.1c wormpep=CE52653 gene=WBGene00007064 locus=rga-9 status=Confirmed uniprot=A0A2X0T1Z3 insdc=SPS41576.1 product="Rho-GAP domain-containing protein" 3R5.1a wormpep=CE24758 gene=WBGene00007065 locus=pot-3 status=Confirmed uniprot=G5EFG7 insdc=CAA21777.2 product="POT1PC domain-containing protein" ` Does this help?
Version is the latest (2.5.4).
All the best, Jenni
mduchoslav
Hello,
I have been investigating some odd results for an OrthoFinder run using NCBI data and new annotations from about 15 species of ducks and other fowl, and I think that I have found a bug in the primary_transcript.py script, specifically the logic used to decide if two proteins are isoforms of the same gene or not in NCBI data.
Looking at Anas platyrhynchos (assembly GCF_003850225.1 and the corresponding NCBI annotation), I see that the primary_transcript.py script returns 18,619 proteins, which should be one isoform per gene, for input to OrthoFinder. However, that particular annotation only has 15,747 genes according to the NCBI annotation report.
Digging further, I can find examples where there are multiple isoforms of the same gene that are not labeled as "isoform X" or "isoform " in the GFF file. Here is an example for the LHFPL2 gene:
For clarity I have only included part of the GFF. But if we look up those XP ids, we see that at least some of them are not indicated as isoforms, e.g.:
Based on the GFF file, these two proteins are clearly isoforms of the same gene, but the primary_transcript.py script treats them as separate genes, since they do not have 'isoform' in their name.
There does not seem to be an obvious solution to this based only on the information provided in the fasta header for protein files downloaded from NCBI; I suspect GFF parsing is the only reliable approach to accurately identify the longest protein from NCBI annotations.