Retrieve the FASTA transcript sequences of a list of NCBI genes

Question

[Sibs]

Hello everyone,

I want to design a microarray probe set for Coccomyxa subellipsoidea C-169 which has been fully sequenced. I need to make a target sequence data file in FASTA format or a TDT file of GenBank accessions. I got a list of 10091 potential genes to design the microarray (NCBI-->Gene). How can I make this file? Do you know any step by step guideline that I can use? Thanks

[pawlowac]

Do you need to pull the genes themselves using identifiers or just get the transcript sequence from your fasta file?

Answer 1

[umer.zeeshan.ijaz]

Say you are able to extract the identifiers with one identifier/line (perhaps in a file). I am just using the echo to generate the list:

$ echo -e "A0QTF8\nA0QTG1\nA0QU63"
A0QTF8
A0QTG1
A0QU63

Then you can use the following one-liner to get the sequences:

$ echo -e "A0QTF8\nA0QTG1\nA0QU63" | while read l; do echo -e ">"$l"\n"$(curl -s http://www.uniprot.org/uniprot/$l.txt | awk '/SQ/,/\/\//{if ($0!~/^\/\// && $0!~/^SQ/) {gsub(" ","",$0); printf $0}}'); done
>A0QTF8
MDLINGMGTSPGYWRTPREPGNDHRRARLDVMAQRIVITGAGGMVGRVLADQAAAKGHTVLALTSSQCDITDEDAVRRFVANGDVVINCAAYTQVDKAEDEPERAHAVNAVGPGNLAKACAAVDAGLIHISTDYVFGAVDRDTPYEVDDETGPVNIYGRTKLAGEQAVLAAKPDAYVVRTAWVYRGGDGSDFVATMRRLAAGDGAIDVVADQVGSPTYTGDLVGALLQIVDGGVEPGILHAANAGVASRFDQARATFEAVGADPERVRPCGSDRHPRPAPRPSYTVLSSQRSAQAGLTPLRDWREALQDAVAAVVGATTDGPLPSTP
>A0QTG1
MSAAANAEHGAADRVEILPVPGLPEFRPGDDLVGSLAEAAPWLRDGDVLVVTSKVVSKCEGRIVAAPSDPEERDTLRRKLIDDEAVRVLARKGRTLITENAIGLVQAAAGVDGSNVGSTELALLPVDPDRSAATLREGLRERLGVTVGVVITDTMGRAWRTGQTDFAIGASGLTVLQGYAGSRDRHGNELVVTEVAVADEIAAAADLVKGKLTAIPVAVVRGLRLPDDGSTAHRLVRAGEDDLFWLGTAEAIELGRRQAQLLRRSVRRFSAEPVPHDAIEAAVGEALTAPAPHHTRPVRFVWVQDSETRTRLLDRMKEQWRADLTADGLDADAVDRRVARGQILYDAPELVIPFLVPDGAHSYPDDARTAAEHTMFTVAVGAAVQGLLVALAVRDIGSCWIGSTIFAADLVRAELELPDDWEPLGAIAIGYPEQTPQPLGPRDPVPTDELLVRK
>A0QU63
MTKKSASSNNKVVATNRKARHNYTILDTYEAGIVLMGTEVKSLREGQASLADAFATVDDGEIWLRNVHIAEYHHGTWTNHAPRRNRKLLLHRKQIDNLIGKIRDGNLTLVPLSIYFTDGKVKVELALARGKQAHDKRQDLARRDAQREVIRELGRRAKGKI

Best Wishes,
Umer

Answer 2

[Dan D]

Step 1: Gene a track from the UCSC table browser, which contains the following for the hg19 genome:

chromosome
txStart (transcription start)
txEnd (transcription end)
gene name

Step 2:

With that data in hand, and with a local copy of the hg19 reference genome, use samtools faidx along with the chromosome, start, and end positions to get the FASTA sequence for each gene. By simply launching a separate samtools process for each gene and capturing the output, you can programmatically build a FASTA containing all of the gene sequences

If you need clarification or coding help for any of these steps, just let me know what specific questions you have.

[David Langenberger]

You can directly tell the UCSC Genome Browser to give you the sequences instead of the genomic loci.

Just select 'sequence' as 'output format'.

Answer 3

[st.ph.n]

This python snippet, using BioPython, will pull the sequence record in genbank format from NCBI using an input file containing a list of accession numbers.

import sys
from Bio import Entrez
from Bio import SeqIO

Entrez.email = #put your email address here

input_list = sys.argv[1]

with open(input_list, "r") as f:
        idlist = [line.strip() for line in f]

handle = Entrez.efetch(db="nucleotide", rettype="gb", retmode="text",id=idlist, email=Entrez.email)

seq_rec_list = []
for seq_record in SeqIO.parse(handle, "gb"):
    seq_rec_list.append(seq_record)

out_handle = open("output.fasta", "w")
SeqIO.write(seq_rec_list, out_handle, "fasta")
out_handle.close()