Minimal DIAMOND blastx output for Roslin genome?

[yaaminiv]

I ran DIAMOND blastx a few days ago and am just getting around to checking the output. I wanted to get Uniprot accession codes for the Roslin genome so I could then get GOterms and perform a gene enrichment. My script can be found here.

I counted the lines of my output file with wc -l (found here), and I only have 209 matches! That doesn't seem right...Should I run DIAMOND blastx again with different settings?

[yaaminiv]

Realized I'm likely using a Uniprot-SwissProt database that was created with more specific taxon information (based on the readme associated with the file). I'm making my own database for blastx to see if that provides more hits.

[yaaminiv]

I made my own database but I'm still only getting 209 matches. I wasn't sure if I was making my database correctly, so I tried adding -dbtype prot to the makedb command, but that just failed.

[yaaminiv]

Here's the DIAMOND blastx output file with 209 lines I generated after I downloaded the Uniprot-SwissProt database: https://gannet.fish.washington.edu/spartina/project-oyster-oa/data/Cg-roslin/20210601-cgigas-roslin-mito-blastx.outfmt6

edited to include line numbers

[kubu4]

Can you provide a link to the input FastA? Or, change permissions (chmod 755 /gscratch/scrubbed/yaaminiv/data/cgigas_uk_roslin_v1_genomic-mito.fa) on Mox so we can view?

EDITED: Fixed markdown formatting.

[sr320]

Also, isn't this genome already annotated? Why do you need to BLAST it?

[yaaminiv]

Also, isn't this genome already annotated? Why do you need to BLAST it?

This is the format of the annotation provided:

(base) yaaminivenkataraman@big-wave ~ % head /Users/yaaminivenkataraman/Downloads/GCF_902806645.1_cgigas_uk_roslin_v1_protein.faa
>NP_001292217.1 uncharacterized LOC105326593 precursor [Crassostrea gigas]
MKVVPSLFIFIICVFISSGQRPDPFPKGCKAQIDVVFVIDGSDSISESDFQTLRTSISRIVDGFHIGSGETRMGIIVYSK
GVAFSVPLSYDPVYLKDQASIMPHPREGTNTHLGIEEMIDMFKKDKRDGVPMAGVVVTDGISKEKEKTLLQSRLARDLGI
NMFSVGVGRYTEEEELRGIASNPDQAIKVESFDELLKILSKLVQLVCPNKCMMPGVIAYPNDVSKNCRLYWKCEGEESKL
TCCPRGFSFSAPVQSCIPDPKCVEPCGDEGPICNKRPSVYQPTIYEELIEGYGWVQRSCPPGTAYDRVTCGCTITQTPPP
PRRVCRVLVHIPFDIDCVDTSGNGFLIKNHGVKFTRTGLALFEGKAKLVIPNIQRYLGSNFLVKMRYKEFPSFETQGLLS
NGDCYTPMTLQLIKDSIRHTYKIENSYRQRTGFTIPTNFKPWNEVTFSHDGNTLAGSVNGIQNEKRSNGPIQRTPYGISI
GFVDGYKHFQGYMDYVTIYKCANMNVYNKY
>NP_001292218.1 prostaglandin G/H synthase 2-like precursor [Crassostrea gigas]
MWWRAVLLVSFIAPALVSGGHNPCCSYPCENNGVCMTQGGNYYCDCANLDFYGKHCEIPTFMKKIKLFLKPSQSTLHNLL

I want to do a GOterm enrichment, but without any Uniprot Accession codes in the annotation, I don't know how to get GOterms for each gene.

[yaaminiv]

Did the SLURM email notification indicate an exit code other than "1"?

I got Failed ExitCode 134 several times, and 137 once

[kubu4]

This is the format of the annotation provided:

That's the FastA file.

Annotation files (e.g. GTF, GFF, etc) are here:

https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/806/645/GCF_902806645.1_cgigas_uk_roslin_v1/

Here's the genome GFF:

https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/806/645/GCF_902806645.1_cgigas_uk_roslin_v1/GCF_902806645.1_cgigas_uk_roslin_v1_genomic.gff.gz

[yaaminiv]

Annotation files (e.g. GTF, GFF, etc) are here

I've used those files before to create genome feature tracks, but they didn't have the information I needed for my specific task. I think the ID mapping is what I'll try, since I need Uniprot Accession codes to get GOterms, and I already have protein IDs in the FastA file.

edited to clarify goals

[sr320]

ID Mapping ... Sam has a script...

[yaaminiv]

@kubu4 Where can I find an ID mapping script?

[kubu4]

This Jupyter Notebook will allow you to get what you want:

https://github.com/RobertsLab/code/blob/master/notebooks/sam/20210601_ssal_gff-annotations.ipynb

Read through it and post any questions you have. It will require the genomic GFF file as input.

You'll likely want to modify it, to limit the final output to just the info you need/want to use downstream.