This is a test version of Biostars. For the public version, visit https://www.biostars.org.
struggle to get fasta files from ucsc goldenPath

Hi all, I am very interested in get fasta seq (mrna) from a given genes IDs. My first approach has been to download from here https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ the refMrna.fa.gz file and then use bedtools getfasta -rna to retrieve it:

e.g: bedtools getfasta -fi refMrna.fa -rna -bed test.bed -fo test.fa.out

(bed file is eg: chr8 9999999 9999998)

However, after getting the error: 'WARNING. chromosome (chr8) was not found in the FASTA file. Skipping.' I decided to explore the ref mrna file and I realised the gene IDs are in format refseq format:

e.g: grep -o -E "^>\w+" refMrna.fa | tr -d ">" > NM_001383, NR_164501 ..

Now I am wondering if there is a smartest way to get the mrna fasta seq or if I should transform the gene ID using the UCSC tool (RefSeq all) which provides an output with lots of refseqs for the given gene and use that file to be passed as a bed file.

Any feedbacks on this issue?

Thanks!

ucsc getfasta fasta
• 2,745 views • link
written by Lila M Lila M

If your intervals are for main chromosomes why are you using the mRNA file instead of the entire genome? If you have the gene ID's you could also use the table browser to get the sequence.

• 0 views • link
written by GenoMax GenoMax

It is because I want to pull the mRNA fasta files and not DNA (which I guess is the case if I pull it from the entire genome)

• 0 views • link
written by Lila M Lila M

mRNA file is not going to have the same co-ordinates as the main chromosomes. Can you clarify what coordinates are in your BED file?

• 0 views • link
written by GenoMax GenoMax

So this is my approach: I want to interrogate gene PVT1, I know the coordinates for it is chr8:127794532-128101252 In UCSC table I selected the following options:

group: mRNA and EST
track: Human mRNAs
table: all_mrna

This prints out a big file, but I can select the coordinates of interest (my BED file <test.bed> ), which are:

chr8    127794532   128101253 LF385466
chr8    127794532   128101253 LF385467
chr8    127794532   128101253 MA621043
chr8    127794532   128101253 MA621044

Using that 'portion' of the file I can then do:

bedtools getfasta -fi mrna.fa -rna -bed test.bed  -fo test.fa.out

Is this approach not correct? What is your suggestion then?

• 0 views • link
written by Lila M Lila M

If coordinates in your BED file refer to chromosomal locations then you need to use the whole genome file and get those sections by the method you noted above.

That said you should be able to get the sequence using table browser directly.

• 0 views • link
written by GenoMax GenoMax

Could you please elaborate how to retrieve the sequences directly from the table browser? Thanks!

• 0 views • link
written by Lila M Lila M

Can you select the output format as sequence instead of BED?

• 0 views • link
written by GenoMax GenoMax

Of course ... But as 'expected' the range range=chr8:127794533-128101253 in the output sequence file belongs to: LF385466, LF385467, MA621043 and MA621044, which exactly the same sequences as my approach using the mrna.fa file.

• 0 views • link
written by Lila M Lila M

Not sure if I am missing something but you want 4 copies of the same sequence with a different fasta header?

• 0 views • link
written by GenoMax GenoMax

They are +ve and -ve strands, but identical two by two ... I can filter that out later on. This is probably a very naive question, but what the header (LF385466, LF385467 ...) means? I know the are the qnames, but why they can be identical sequences and different headers?

• 0 views • link
written by Lila M Lila M

Yes! I've just realised about it! Thank you so much for following this issue!

• 0 views • link
written by Lila M Lila M

0 answers

No answers yet.

Log in to answer this question.

More posts like this