Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-12T09:00:56.510Z
Answers.Text: 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.
Answers.UserDisplayName: CheesyGnocchi
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-12T11:54:42.363Z
Answers.Text: 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.
Answers.UserDisplayName: vkkodali
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-11T15:46:26.017Z
Answers.LastActivityDate: 2020-08-11T15:46:26.017Z
Answers.UserDisplayName: vkkodali
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-13T18:18:09.340Z
Answers.LastActivityDate: 2020-08-13T18:18:09.340Z
Answers.UserDisplayName: aerijman
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-13T10:27:32.137Z
Answers.LastActivityDate: 2020-08-13T10:27:32.137Z
Answers.UserDisplayName: StupidWolf
Body: <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>
  
Tags: <fasta><python><reference-genome><gff3>
Title: Get gene sequence based on the annotation
ThreadText: # 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())
  
  ```
  
  
  ---
  
  