Answers.Body: <p>I proceeded as follows:</p>
  <pre><code>import gzip
  import sys
  
  qtls=['sqtls','eqtls']
  tissues=['Adipose_Subcutaneous', 'Adipose_Visceral_Omentum', 'Adrenal_Gland', 'Artery_Aorta', 'Artery_Coronary', 'Artery_Tibial', 'Brain_Amygdala', 'Brain_Anterior_cingulate_cortex_BA24', 'Brain_Caudate_basal_ganglia', 'Brain_Cerebellar_Hemisphere', 'Brain_Cerebellum', 'Brain_Cortex', 'Brain_Frontal_Cortex_BA9', 'Brain_Hippocampus', 'Brain_Hypothalamus', 'Brain_Nucleus_accumbens_basal_ganglia', 'Brain_Putamen_basal_ganglia', 'Brain_Spinal_cord_cervical_c-1', 'Brain_Substantia_nigra', 'Breast_Mammary_Tissue', 'Cells_Cultured_fibroblasts', 'Cells_EBV-transformed_lymphocytes', 'Colon_Sigmoid', 'Colon_Transverse', 'Esophagus_Gastroesophageal_Junction', 'Esophagus_Mucosa', 'Esophagus_Muscularis', 'Heart_Atrial_Appendage', 'Heart_Left_Ventricle', 'Kidney_Cortex', 'Liver', 'Lung', 'Minor_Salivary_Gland', 'Muscle_Skeletal', 'Nerve_Tibial', 'Ovary', 'Pancreas', 'Pituitary', 'Prostate', 'Skin_Not_Sun_Exposed_Suprapubic', 'Skin_Sun_Exposed_Lower_leg', 'Small_Intestine_Terminal_Ileum', 'Spleen', 'Stomach', 'Testis', 'Thyroid', 'Uterus', 'Vagina', 'Whole_Blood']
  
  gene_old=''
  
  
  def index(tissue,qtl,gene_old):
          counter_old,counter=1,1
          if qtl=='sqtls':
                  path='/path/to/file/'+tissue+'.v8.sqtl_allpairs.txt.gz'
          else:
                  path='/path/to/file/'+tissue+'.allpairs.txt.gz'
          f = gzip.open(path, 'r')
          next(f)
          for line in f:
                  try:
                          line=line.decode()
                          line=line[:line.find('\t')]
                          gene=line[line.rfind('ENS'):].split('.')[0]
                          #print(gene,gene_old,counter_old,counter)
                          with open(qtl+'_'+tissue+'.idx','a') as indexed:          
                                   indexed.write(str(gene_old)+'\t'+str(counter_old)+'\t'+str(counter)+'\n')
                          gene_old=gene
                          counter +=1
                          counter_old=counter
                  except IndexError:
                          print('tissue,qtl,gene_old,counter')
                          counter+=1
          f.close()
  
  
  for qtl in qtls:
          for tissue in tissues:
                 qtl=row.qtls
                 tissue=row.tissues
                 print(qtl, tissue)
                 index(tissue,qtl,gene_old)
  </code></pre>
  <p>The above script writes an index for each gene which should look as follows (took 2 days to finish all files):</p>
  <pre><code>ensid    row_start    row_end
  ......
  ENSG00000230337     1526391 1534329
  ENSG00000171819 1534330 1542321
  ENSG00000198793    1542322 1550178
  ENSG00000120942  1550179 1558010
  .....
  </code></pre>
  <p>And to get a gene of interest you will have to look up the coordinates for an ensemblid/gene and for example:</p>
  <pre><code>zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' &gt; ENSG00000120942.txt
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-02T13:23:41.480Z
Answers.LastActivityDate: 2020-09-02T14:39:46.240Z
Answers.LastEditDate: 2020-09-02T14:39:46.240Z
Answers.UserDisplayName: which_command
Answers.Body: <p>Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-31T08:12:18.343Z
Answers.Text: Hi @Colin thanks for getting to me- I just want to be able to take out all variants inside a particular gene- for example if the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to `zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt`
  Hence I thought tabix would be the right one for this
Answers.UserDisplayName: which_command
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-31T15:48:25.117Z
Answers.Text: Here is an example of a tabix file that has been hacked up to allow "free text query". http://target.wustl.edu/ENCFF000ARO_ALL.iteres.loci.gz  So they can query `tabix file.txt "A(n)"` and because there are many lines starting with A(n) it will return all of them
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-31T16:08:17.487Z
Answers.Text: Can also check out other options for text indexing https://unix.stackexchange.com/questions/3086/command-line-friendly-full-text-indexing
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-31T06:03:53.403Z
Answers.LastActivityDate: 2020-08-31T06:03:53.403Z
Answers.UserDisplayName: Colin D
Body: <p>I have a gtex variant file- the head of which looks as follows:</p>
  <pre><code>phenotype_id                                    variant_id     
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...
  </code></pre>
  <p>Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file:
  E.g. If the variants in gene &quot;ENSG00000148481- MINDY3&quot; are on line number 87528225 to 87536766- effectively what I would like is the equivalent to <code>zcat file.gz | sed -n '87528225,87536766p' &gt; MINDY3.txt</code>. However <code>zgrep ENSG00000148481 file.gz</code> is just as fast as the above...</p>
  <p>Hence I thought tabix would be the right tool for this-</p>
  <p>I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:</p>
  <pre><code>zcat gtex.txt.gz | bgzip &gt; gtex.txt.bgz
  </code></pre>
  <p>However I am not quite sure how to proceed from there given that the data is not tab-delimited.</p>
  <p>As a trial I tried the first 1000 lines:</p>
  <pre><code>zcat gtex.txt.gz | head -n 1000 | bgzip &gt; gtex_1000.txt.bgz
  
  ./tabix -p bed gtex_1000.gz #index as a bed file
  [get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......
  
  ./tabix -p vcf gtex_1000.gz #index as a vcf file
  Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:
  
  ./tabix test.gz chr1:15000:17000
  </code></pre>
  <p>It returns nothing.</p>
  <p>I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time
  ... Does anyone know of a trick to index the files with unconventional delimiting?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-29T02:21:28.870Z
Comments.Text: Suppose that you could successfully creat an index file for your txt file, which program would be used to read it?
Comments.UserDisplayName: Phoenix Mu
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-31T16:25:54.800Z
Comments.Text: @PhoenixMu see the way that I have used `sed` in the updated question... this was my plan... although this was the same speed as `grep`
Comments.UserDisplayName: which_command
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-01T17:13:59.130Z
Comments.Text: I think tabix only works for bed and vcf format. Your file is not in bed format?
Comments.UserDisplayName: Phoenix Mu
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-02T09:06:44.663Z
Comments.Text: Yes it's not in bed format as the coordinate columns are : delimited and the file is too big to justify producing a copy with properly delimited columns- I am in the process of making an index file myself and will post an answer if it works
Comments.UserDisplayName: which_command
Tags: <parsing><tabix><gtex>
Title: tabix indexing a GTEx variant file- ':' delimited
ThreadText: # tabix indexing a GTEx variant file- ':' delimited
  
  (Asked by: username_0 on 2020-08-27)
  
  I have a gtex variant file- the head of which looks as follows:
  
  ```
  phenotype_id                                    variant_id     
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
  chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...
  
  ```
  
  Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file:
  E.g. If the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to `zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt`. However `zgrep ENSG00000148481 file.gz` is just as fast as the above...
  
  Hence I thought tabix would be the right tool for this-
  
  I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:
  
  ```
  zcat gtex.txt.gz | bgzip > gtex.txt.bgz
  
  ```
  
  However I am not quite sure how to proceed from there given that the data is not tab-delimited.
  
  As a trial I tried the first 1000 lines:
  
  ```
  zcat gtex.txt.gz | head -n 1000 | bgzip > gtex_1000.txt.bgz
  
  ./tabix -p bed gtex_1000.gz #index as a bed file
  [get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......
  
  ./tabix -p vcf gtex_1000.gz #index as a vcf file
  Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:
  
  ./tabix test.gz chr1:15000:17000
  
  ```
  
  It returns nothing.
  
  I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time
  ... Does anyone know of a trick to index the files with unconventional delimiting?
  
  **Tags:** parsing, tabix, gtex
  
  > **Comments:**
  > **username_1**: Suppose that you could successfully creat an index file for your txt file, which program would be used to read it?  
  > **username_0**: @username_1 see the way that I have used `sed` in the updated question... this was my plan... although this was the same speed as `grep`  
  > **username_1**: I think tabix only works for bed and vcf format. Your file is not in bed format?  
  > **username_0**: Yes it's not in bed format as the coordinate columns are : delimited and the file is too big to justify producing a copy with properly delimited columns- I am in the process of making an index file myself and will post an answer if it works  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2020-09-02. Score: 0 | Accepted answer
  
  I proceeded as follows:
  
  ```
  import gzip
  import sys
  
  qtls=['sqtls','eqtls']
  tissues=['Adipose_Subcutaneous', 'Adipose_Visceral_Omentum', 'Adrenal_Gland', 'Artery_Aorta', 'Artery_Coronary', 'Artery_Tibial', 'Brain_Amygdala', 'Brain_Anterior_cingulate_cortex_BA24', 'Brain_Caudate_basal_ganglia', 'Brain_Cerebellar_Hemisphere', 'Brain_Cerebellum', 'Brain_Cortex', 'Brain_Frontal_Cortex_BA9', 'Brain_Hippocampus', 'Brain_Hypothalamus', 'Brain_Nucleus_accumbens_basal_ganglia', 'Brain_Putamen_basal_ganglia', 'Brain_Spinal_cord_cervical_c-1', 'Brain_Substantia_nigra', 'Breast_Mammary_Tissue', 'Cells_Cultured_fibroblasts', 'Cells_EBV-transformed_lymphocytes', 'Colon_Sigmoid', 'Colon_Transverse', 'Esophagus_Gastroesophageal_Junction', 'Esophagus_Mucosa', 'Esophagus_Muscularis', 'Heart_Atrial_Appendage', 'Heart_Left_Ventricle', 'Kidney_Cortex', 'Liver', 'Lung', 'Minor_Salivary_Gland', 'Muscle_Skeletal', 'Nerve_Tibial', 'Ovary', 'Pancreas', 'Pituitary', 'Prostate', 'Skin_Not_Sun_Exposed_Suprapubic', 'Skin_Sun_Exposed_Lower_leg', 'Small_Intestine_Terminal_Ileum', 'Spleen', 'Stomach', 'Testis', 'Thyroid', 'Uterus', 'Vagina', 'Whole_Blood']
  
  gene_old=''
  
  def index(tissue,qtl,gene_old):
          counter_old,counter=1,1
          if qtl=='sqtls':
                  path='/path/to/file/'+tissue+'.v8.sqtl_allpairs.txt.gz'
          else:
                  path='/path/to/file/'+tissue+'.allpairs.txt.gz'
          f = gzip.open(path, 'r')
          next(f)
          for line in f:
                  try:
                          line=line.decode()
                          line=line[:line.find('\t')]
                          gene=line[line.rfind('ENS'):].split('.')[0]
                          #print(gene,gene_old,counter_old,counter)
                          with open(qtl+'_'+tissue+'.idx','a') as indexed:          
                                   indexed.write(str(gene_old)+'\t'+str(counter_old)+'\t'+str(counter)+'\n')
                          gene_old=gene
                          counter +=1
                          counter_old=counter
                  except IndexError:
                          print('tissue,qtl,gene_old,counter')
                          counter+=1
          f.close()
  
  for qtl in qtls:
          for tissue in tissues:
                 qtl=row.qtls
                 tissue=row.tissues
                 print(qtl, tissue)
                 index(tissue,qtl,gene_old)
  
  ```
  
  The above script writes an index for each gene which should look as follows (took 2 days to finish all files):
  
  ```
  ensid    row_start    row_end
  ......
  ENSG00000230337     1526391 1534329
  ENSG00000171819 1534330 1542321
  ENSG00000198793    1542322 1550178
  ENSG00000120942  1550179 1558010
  .....
  
  ```
  
  And to get a gene of interest you will have to look up the coordinates for an ensemblid/gene and for example:
  
  ```
  zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' > ENSG00000120942.txt
  
  ```
  
  
  ---
  
  ### Answer by username_2 on 2020-08-31. Score: 0
  
  Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?
  
  > **Comments:**
  > **username_0**: Hi @username_2 thanks for getting to me- I just want to be able to take out all variants inside a particular gene- for example if the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to `zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt` Hence I thought tabix would be the right one for this  
  > **username_2**: Here is an example of a tabix file that has been hacked up to allow "free text query". http://target.wustl.edu/ENCFF000ARO_ALL.iteres.loci.gz  So they can query `tabix file.txt "A(n)"` and because there are many lines starting with A(n) it will return all of them  
  > **username_2**: Can also check out other options for text indexing https://unix.stackexchange.com/questions/3086/command-line-friendly-full-text-indexing  
  
  ---
  
  