title: How do I get the microarray data?
body: <p>thank you for your help. I want to use the following python code, to read and process data from an affymetrix microarray data set. I want to elucidate differential gene expression in disease conditions of Crohn's disease and Ulcerative colitis, in mononuclear cells. The code runs perfectly, but when I try to see the content of X, I get an empty array at the output (like this : array([], dtype=float64)), which of course is not useful. Here is a link to the raw data set : <a href="https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS1615" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS1615</a>
  I have tried long to figure out why I have an empty and unprocessable output, but to no avail. Here is the code:</p>
  
  <pre><code>import gzip
  import numpy as np
  
  """
  Read in a SOFT format data file.  The following values can be exported:
  
  GID : A list of gene identifiers of length d
  SID : A list of sample identifiers of length n
  STP : A list of sample descriptions of length d
  X   : A dxn array of gene expression values
  """
  #path to the data file
  fname = "../data/GDS1615_full.soft.gz"
  
  ## Open the data file directly as a gzip file
  with gzip.open(fname) as fid:
      SIF = {}
      for line in fid:
          if line.startswith(line, len("!dataset_table_begin")):
              break
          elif line.startswith(line, len("!subject_description")):
              subset_description = line.split("=")[1].strip()
          elif line.startswith(line, len("!subset_sample_id")):
              subset_ids = [x.strip() for x in subset_ids]
              for k in subset_ids:
                  SIF[k] = subset_description
      ## Next line is the column headers (sample id's)
      SID = next(fid).split("\t")
  
      ## The column indices that contain gene expression data
      I = [i for i,x in enumerate(SID) if x.startswith("GSM")]
  
      ## Restrict the column headers to those that we keep
      SID = [SID[i] for i in I]
  
      ## Get a list of sample labels
      STP = [SIF[k] for k in SID]
  
      ## Read the gene expression data as a list of lists, also get the gene
      ## identifiers
      GID,X = [],[]
      for line in fid:
  
          ## This is what signals the end of the gene expression data
          ## section in the file
          if line.startswith("!dataset_table_end"):
              break
  
          V = line.split("\t")
  
          ## Extract the values that correspond to gene expression measures
          ## and convert the strings to numbers
          x = [float(V[i]) for i in I]
  
          X.append(x)
          GID.append(V[0] + ";" + V[1])
  X = np.array(X)
  
  ## The indices of samples for the ulcerative colitis group
  UC = [i for i,x in enumerate(STP) if x == "ulcerative colitis"]
  
  ## The indices of samples for the Crohn's disease group
  CD = [i for i,x in enumerate(STP) if x == "Crohn's disease"]
  </code></pre>
  
  <p>At the console, I get such output:
  X
  Out[94]: array([], dtype=float64)</p>
  
  <p>X.shape
  Out[95]: (0,)</p>
  
  <p>Thank you once more for your suggestions.</p>