text: user3186: I want to find the length of a list of genes, of Homo sapiens, that is reported in the GEO database. I have gathered the ensembl id's of those genes.
  I understand this information can be parsed from the start and end position given in GTF file available in ensembl.org. 
  I would like to ask for advice on how to programmatically download the gtf files with ensembl id's.
  On the ensembl webpage a couple of options are listed, API and ftp access. I would like to know if sample scripts are available.
  Update: I found a script here, it shows how to parse the data from gtf file. Could someone help in understanding how the code can be modified to obtain gene length from start and end position?
  Edit 1: 
  By start and end position, I refer to the start and end values given in the "fields" section here and the file from which I want to parse the values is the second last file available in this link.
  For instance, the end - start value of the gene with ensembl id ENSG00000131482,
  gene :G6PC
  Edit 2: I get the gene names and the expression values from the SOFT files documented in GEO, link
  The following is a sample output from SOFT file. It gives me the probe_set(ID_REF) , Gene name and expression values from samples (GSMXXXXXX). Here we see fro G6PC, values are reported from two probe_sets. Can we get the transcript id's from the probe_sets?
  Obtaining the transcript id will help me in finding the length of the specific transcript. 
  > Table(gds)[2449,0:5]
             ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
  2449 1555612_s_at       G6PC   4.39472   4.52235   4.26869
  > Table(gds)[16399,0:5]
           ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
  16399 206952_at       G6PC   2.76693   3.45107    2.9815
  
  
  user298: If you really just want the start and end positions of all genes in a list of genes from Homo_sapiens.GRCh38.93.gtf.gz, you could do this in awk. For example, if your list of genes is:
  $ cat genes
  ENSG00000131482
  ENSG00000223972
  ENSG00000228794
  ENSG00000187961
  
  You can do:
  $ awk '{ 
          if(NR==FNR){
              wantedGenes[$1]++ 
          } 
          else{
              if($3=="gene"){
                  gsub(/";*/,"",$10); 
                  if($10 in wantedGenes){
                      print $10,$4,$5
                  } 
              }
          }
        }' genes Homo_sapiens.GRCh38.93.gtf 
  ENSG00000223972 11869 14409
  ENSG00000228794 825138 859446
  ENSG00000187961 960587 965715
  ENSG00000131482 42900797 42913369
  
  Explanation:
  
  if(NR==FNR){ wantedGenes[$1]++ } : NR is the input line number and FNR is the current line number of this input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array wantedGenes.
  if($3=="gene"){ : if this line is a 'gene' entry.
  gsub(/";*/,"",$10); : the gene names in the GTF file are quoted and have an ; after the quote. For example: "ENSG00000187961";. Since we don't have this in the list, we need to remove it so we can match the gene name. This gsub command will remove all " followed by 0 or more ; (the ;* means 0 or more occurrences of ;) from the 10th field of each line. So, it will replace "ENSG00000187961"; with ENSG00000187961.
  if($10 in wantedGenes){print $10,$4,$5} : if the 10th field (the modified gene name from above) is in the wantedGenes array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).
  
  You can easily modify this to print the length instead of the start and end positions:
  $ awk '{ 
          if(NR==FNR){
              wantedGenes[$1]++ 
          } 
          else{
              if($3=="gene"){
                  gsub(/";*/,"",$10); 
                  if($10 in wantedGenes){
                      print $10,$5-$4
                  } 
              }
          }
        }' genes Homo_sapiens.GRCh38.93.gtf 
  ENSG00000223972 2540
  ENSG00000228794 34308
  ENSG00000187961 5128
  ENSG00000131482 12572
  
  As requested in chat, here's a python version of the same thing:
  #!/usr/bin/env python3
  import re
  import sys
  
  wantedGenes = []
  with open(sys.argv[1], mode='r') as f:
      wantedGenes = f.readlines()
  
  wantedGenes = [ x.rstrip() for x in wantedGenes]
  
  print('# Gene\tstart\tend\tlength')
  with open(sys.argv[2], mode='r') as f:
      for line in f:
          if not line.startswith('#'): 
              fields = line.split('\t')
              annotations = fields[8].split(' ')
              geneName = re.sub('[";]', '', annotations[1])
              if geneName in wantedGenes and fields[2] == 'gene':
                  print('%s\t%s\t%s\t%d' %
                        (geneName, fields[3],fields[4],
                         int( int(fields[4]) - int(fields[3]) )
                        ))
  
  Save the file as foo.py and run with the file names as arguments:
  python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf
  
  
  user383: You could also do this in R with the GenomicFeatures library.
  library(GenomicFeatures)
  
  ## make TxDb from GTF file 
  txdb <- makeTxDbFromGFF('Homo_sapiens.GRCh38.93.gtf')
  
  ## get gene information
  all.genes <- genes(txdb)
  
  ## import your list of gene names
  my.genes <- c('ENSG00000141510','ENSG00000184571','ENSG00000011007')
  
  ## get the length of each of those genes
  my.genes.lengths <- width(all.genes[my.genes])
  ## put the names back on the lengths
  names(my.genes.lengths) <- my.genes
  
  ## print lengths
  print(my.genes.lengths)
  
  ENSG00000141510 ENSG00000184571 ENSG00000011007
        25772           55687           18905
  
  
original_text: 3186: <p>I want to find the length of a list of genes, of Homo sapiens, that is reported in the GEO database. I have gathered the ensembl id's of those genes.</p>
  
  <p>I understand this information can be parsed from the start and end position given in GTF file available in ensembl.org. </p>
  
  <p>I would like to ask for advice on how to programmatically download the gtf files with ensembl id's.</p>
  
  <p>On the ensembl webpage a couple of options are listed, API and ftp access. I would like to know if sample scripts are available.</p>
  
  <p>Update: I found a script <a href="https://github.com/dpryan79/Answers/blob/master/SEQanswers_42420/GTF2LengthGC.R" rel="noreferrer">here</a>, it shows how to parse the data from gtf file. Could someone help in understanding how the code can be modified to obtain gene length from start and end position?</p>
  
  <p>Edit 1: 
  By start and end position, I refer to the start and end values given in the "fields" section <a href="https://asia.ensembl.org/info/website/upload/gff.html" rel="noreferrer">here</a> and the file from which I want to parse the values is the second last file available in this <a href="ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/" rel="noreferrer">link</a>.<br>
  For instance, the end - start value of the gene with ensembl id ENSG00000131482,
  gene :G6PC</p>
  
  <p>Edit 2: I get the gene names and the expression values from the SOFT files documented in GEO, <a href="https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS4027" rel="noreferrer">link</a></p>
  
  <p>The following is a sample output from SOFT file. It gives me the probe_set(ID_REF) , Gene name and expression values from samples (GSMXXXXXX). Here we see fro G6PC, values are reported from two probe_sets. Can we get the transcript id's from the probe_sets?</p>
  
  <p>Obtaining the transcript id will help me in finding the length of the specific transcript. </p>
  
  <pre><code>&gt; Table(gds)[2449,0:5]
             ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
  2449 1555612_s_at       G6PC   4.39472   4.52235   4.26869
  &gt; Table(gds)[16399,0:5]
           ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
  16399 206952_at       G6PC   2.76693   3.45107    2.9815
  </code></pre>
   298: <p>If you really just want the start and end positions of all genes in a list of genes from <a href="ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.93.gtf.gz" rel="nofollow noreferrer">Homo_sapiens.GRCh38.93.gtf.gz</a>, you could do this in <code>awk</code>. For example, if your list of genes is:</p>
  
  <pre><code>$ cat genes
  ENSG00000131482
  ENSG00000223972
  ENSG00000228794
  ENSG00000187961
  </code></pre>
  
  <p>You can do:</p>
  
  <pre><code>$ awk '{ 
          if(NR==FNR){
              wantedGenes[$1]++ 
          } 
          else{
              if($3=="gene"){
                  gsub(/";*/,"",$10); 
                  if($10 in wantedGenes){
                      print $10,$4,$5
                  } 
              }
          }
        }' genes Homo_sapiens.GRCh38.93.gtf 
  ENSG00000223972 11869 14409
  ENSG00000228794 825138 859446
  ENSG00000187961 960587 965715
  ENSG00000131482 42900797 42913369
  </code></pre>
  
  <h3>Explanation:</h3>
  
  <ul>
  <li><code>if(NR==FNR){ wantedGenes[$1]++ }</code> : <code>NR</code> is the input line number and <code>FNR</code> is the current line number of <em>this</em> input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array <code>wantedGenes</code>.</li>
  <li><code>if($3=="gene"){</code> : if this line is a 'gene' entry.</li>
  <li><code>gsub(/";*/,"",$10);</code> : the gene names in the GTF file are quoted and have an <code>;</code> after the quote. For example: <code>"ENSG00000187961";</code>. Since we don't have this in the list, we need to remove it so we can match the gene name. This <code>gsub</code> command will remove all <code>"</code> followed by 0 or more <code>;</code> (the <code>;*</code> means 0 or more occurrences of <code>;</code>) from the 10th field of each line. So, it will replace <code>"ENSG00000187961";</code> with <code>ENSG00000187961</code>.</li>
  <li><code>if($10 in wantedGenes){print $10,$4,$5}</code> : if the 10th field (the modified gene name from above) is in the <code>wantedGenes</code> array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).</li>
  </ul>
  
  <p>You can easily modify this to print the length instead of the start and end positions:</p>
  
  <pre><code>$ awk '{ 
          if(NR==FNR){
              wantedGenes[$1]++ 
          } 
          else{
              if($3=="gene"){
                  gsub(/";*/,"",$10); 
                  if($10 in wantedGenes){
                      print $10,$5-$4
                  } 
              }
          }
        }' genes Homo_sapiens.GRCh38.93.gtf 
  ENSG00000223972 2540
  ENSG00000228794 34308
  ENSG00000187961 5128
  ENSG00000131482 12572
  </code></pre>
  
  <p>As requested in chat, here's a python version of the same thing:</p>
  
  <pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python3
  import re
  import sys
  
  wantedGenes = []
  with open(sys.argv[1], mode='r') as f:
      wantedGenes = f.readlines()
  
  wantedGenes = [ x.rstrip() for x in wantedGenes]
  
  print('# Gene\tstart\tend\tlength')
  with open(sys.argv[2], mode='r') as f:
      for line in f:
          if not line.startswith('#'): 
              fields = line.split('\t')
              annotations = fields[8].split(' ')
              geneName = re.sub('[";]', '', annotations[1])
              if geneName in wantedGenes and fields[2] == 'gene':
                  print('%s\t%s\t%s\t%d' %
                        (geneName, fields[3],fields[4],
                         int( int(fields[4]) - int(fields[3]) )
                        ))
  </code></pre>
  
  <p>Save the file as <code>foo.py</code> and run with the file names as arguments:</p>
  
  <pre><code>python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf
  </code></pre>
   383: <p>You could also do this in R with the GenomicFeatures library.</p>
  
  <pre><code>library(GenomicFeatures)
  
  ## make TxDb from GTF file 
  txdb &lt;- makeTxDbFromGFF('Homo_sapiens.GRCh38.93.gtf')
  
  ## get gene information
  all.genes &lt;- genes(txdb)
  
  ## import your list of gene names
  my.genes &lt;- c('ENSG00000141510','ENSG00000184571','ENSG00000011007')
  
  ## get the length of each of those genes
  my.genes.lengths &lt;- width(all.genes[my.genes])
  ## put the names back on the lengths
  names(my.genes.lengths) &lt;- my.genes
  
  ## print lengths
  print(my.genes.lengths)
  
  ENSG00000141510 ENSG00000184571 ENSG00000011007
        25772           55687           18905
  </code></pre>
  