text: user9361: 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!
  
  user3577: 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
  ```
  
  
original_text: 9361: <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>
   3577: <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>
  