I have a reference genome ref_genome.fna I got by downloading via datasets download genome via its accession number.
GenBank assembly and RefSeq assembly are identical.
When I use RefSeq accession number to find it in RefSeq database, I can see features that were annotated, and those are Gene; mRNA; CDS; ncRNA.
When I go to it's annotation release, link here, I can see that 25,293 genes are protein coding. Likewise, annotation products are available on ftp site, link here.
What I want to do is extract the locations of those 25,293 protein coding genes (alongside any gene identifier) in a single file, eg.
chrom chromStart chromEnd geneID
chr2 100000 155000 RK031
... ... ... ...
For that purpose, what file do I need to download, and what tool do I need to use?
1 answer
Hi Ivan,
To create a table like you described in your post, you can use grep and awk to extract the fields from the feature table. Here's how I did it:
grep "protein_coding" GCF_000372685.2_Astyanax_mexicanus-2.0_feature_table.txt | awk -F '\t' '{print substr($5,0,5)$6,$8,$9,$16}' > GCF_000372685.2_Astyanax_mexicanus_protein_coding.txt
I used grep to extract all lines containing "protein_coding"; I piped the output to awk, and extracted the fields 5 (seq_type, and only the first five characters, so you have chrom), 6 (chromosome number), 8 (start), 9 (end) and 16 (GeneID). The fields 5 and 6 are concatenated, so you have only one column with the chromosome info. Those entries without chromosome info (unplaced scaffolds) will appear as unpla.
I got only 25,117 protein coding from the feature table. Not sure about the discrepancy though.
I hope it helps :)
Log in to answer this question.
You can simply download the transcript and protein sequences for Astyanax mexicanus from NCBI.
Otherwise you can use
AGAT(Extracting genomic feature sequences from GTF/GFF files with AGAT ) for extracting this information from genome file.