<p>You can parse the GFF3 to generate a bed-like file with gene ranges and use that with <code>bedtools getfasta</code> to obtain a multi-FASTA file with gene sequences. An example case is shown below:</p>
<pre><code>## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed 
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS=&quot;\t&quot;;OFS=&quot;\t&quot;}(<span class="math-container">$3=="gene"){print $1,$4-1,$5,"name",1000,$</span>7}' &gt; genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2020-08-12T09:00:56.510Z
Thanks, I'll keep that in mind! But at the moment, I'm operating on anair-gapped server where I don't have the tools and can only use Python.
CheesyGnocchi
CC BY-SA 4.0
2020-08-12T11:54:42.363Z
In that case, you may want to check out [pybedtools](http://daler.github.io/pybedtools/) or [pyfaidx](https://github.com/mdshw5/pyfaidx). In both cases, you should be able to provide a list of ranges along with an indexed FASTA to extract sequences in multi-fasta format.
vkkodali
CC BY-SA 4.0
2020-08-11T15:46:26.017Z
2020-08-11T15:46:26.017Z
vkkodali
<p>You should do:</p>
<pre><code>if sequences_by_chrom[strand] == &quot;-&quot;:    
    rev_comp( sequences_by_chr[chromosome][start:end] )
else:
    sequences_by_chr[chromosome][start:end]
</code></pre>
<p>where</p>
<pre><code>def rev_comp(sequence):
    d = {'A':'T', 'T':'A','C':'G','G':'C'}
    return ''.join([d[i] for i in sequence][::-1])
</code></pre>

CC BY-SA 4.0
2020-08-13T18:18:09.340Z
2020-08-13T18:18:09.340Z
aerijman
<p>DNA sequence is always provided on the plus strand, if the gff tells you it is on the minus strand, you need to reverse complement it, because the start of the gene in your example will be 45290076 and the end will now be 45286331.</p>
<p>I can give you an example below, this is a potential coding sequence from yeast, no introns and on the minus strand:</p>
<pre><code>gunzip -c test.fa.gz  
&gt;1
TCATCCCCTTATATAAGATGGGAAACGAAGGTAAAAGATTAACGATAGCAT
</code></pre>
<p>If we read it in like you did, and if you apply the complement, this is not a coding sequence.</p>
<pre><code>chromosome=&quot;1&quot;
start = 1
end =51
sequences_by_chr[chromosome].complement()[(start-1):end]

Seq('AGTAGGGGAATATATTCTACCCTTTGCTTCCATTTTCTAATTGCTATCGTA', SingleLetterAlphabet())
</code></pre>
<p>Getting the sequence followed by reverse complement gives you a coding sequence, starting with ATG ending with TGA:</p>
<pre><code>sequences_by_chr[&quot;1&quot;][(start-1):end].reverse_complement()

Seq('ATGCTATCGTTAATCTTTTACCTTCGTTTCCCATCTTATATAAGGGGATGA', SingleLetterAlphabet())
</code></pre>

CC BY-SA 4.0
2020-08-13T10:27:32.137Z
2020-08-13T10:27:32.137Z
StupidWolf
<p>I've got the reference genome with Python like so:</p>
<pre><code>sequences_by_chr = {}

with gzip.open(&quot;data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz&quot;, &quot;rt&quot;) as f:
    for seq in SeqIO.parse(f, &quot;fasta&quot;):
        sequences_by_chr[seq.id] = seq.seq
</code></pre>
<p>I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation.
Suppose the following:</p>
<pre><code>chromosome = 1
strand = &quot;-&quot;
start = 45286331
end = 45290076
</code></pre>
<p>Since it's on the &quot;minus&quot; strand, should I do the following:</p>
<p><code>sequences_by_chr[chromosome].complement()[start:end]</code></p>
<p>or should I use the <code>reverse_complement()</code>?</p>
<p>I'm completely confused by this, would appreciate any hint!</p>

<fasta><python><reference-genome><gff3>
Get gene sequence based on the annotation
# Get gene sequence based on the annotation

(Asked by: username_0 on 2020-08-11)

I've got the reference genome with Python like so:
sequences_by_chr = {}

with gzip.open("data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz", "rt") as f:
    for seq in SeqIO.parse(f, "fasta"):
        sequences_by_chr[seq.id] = seq.seq

I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation.
Suppose the following:
chromosome = 1
strand = "-"
start = 45286331
end = 45290076

Since it's on the "minus" strand, should I do the following:
sequences_by_chr[chromosome].complement()[start:end]
or should I use the reverse_complement()?
I'm completely confused by this, would appreciate any hint!

**Tags:** fasta, python, reference-genome, gff3

---

## 3 answers

### Answer by username_1 on 2020-08-11. Score: 2

You can parse the GFF3 to generate a bed-like file with gene ranges and use that with bedtools getfasta to obtain a multi-FASTA file with gene sequences. An example case is shown below:
## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed 
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS="\t";OFS="\t"}($3=="gene"){print $1,$4-1,$5,"name",1000,$7}' > genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
```

> **Comments:**
> **username_0**: Thanks, I'll keep that in mind! But at the moment, I'm operating on anair-gapped server where I don't have the tools and can only use Python.  
> **username_1**: In that case, you may want to check out [pybedtools](http://daler.github.io/pybedtools/) or [pyfaidx](https://github.com/mdshw5/pyfaidx). In both cases, you should be able to provide a list of ranges along with an indexed FASTA to extract sequences in multi-fasta format.  

---

### Answer by username_2 on 2020-08-13. Score: 0

You should do:
if sequences_by_chrom[strand] == "-":    
    rev_comp( sequences_by_chr[chromosome][start:end] )
else:
    sequences_by_chr[chromosome][start:end]

where
def rev_comp(sequence):
    d = {'A':'T', 'T':'A','C':'G','G':'C'}
    return ''.join([d[i] for i in sequence][::-1])


---

### Answer by username_3 on 2020-08-13. Score: 0

DNA sequence is always provided on the plus strand, if the gff tells you it is on the minus strand, you need to reverse complement it, because the start of the gene in your example will be 45290076 and the end will now be 45286331.
I can give you an example below, this is a potential coding sequence from yeast, no introns and on the minus strand:
gunzip -c test.fa.gz  
>1
TCATCCCCTTATATAAGATGGGAAACGAAGGTAAAAGATTAACGATAGCAT

If we read it in like you did, and if you apply the complement, this is not a coding sequence.
chromosome="1"
start = 1
end =51
sequences_by_chr[chromosome].complement()[(start-1):end]

Seq('AGTAGGGGAATATATTCTACCCTTTGCTTCCATTTTCTAATTGCTATCGTA', SingleLetterAlphabet())

Getting the sequence followed by reverse complement gives you a coding sequence, starting with ATG ending with TGA:
sequences_by_chr["1"][(start-1):end].reverse_complement()

Seq('ATGCTATCGTTAATCTTTTACCTTCGTTTCCCATCTTATATAAGGGGATGA', SingleLetterAlphabet())


---

