<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>

CC BY-SA 4.0
2020-09-02T13:23:41.480Z
2020-09-02T14:39:46.240Z
2020-09-02T14:39:46.240Z
which_command
<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>

CC BY-SA 4.0
2020-08-31T08:12:18.343Z
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
which_command
CC BY-SA 4.0
2020-08-31T15:48:25.117Z
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
Colin D
CC BY-SA 4.0
2020-08-31T16:08:17.487Z
Can also check out other options for text indexing https://unix.stackexchange.com/questions/3086/command-line-friendly-full-text-indexing
Colin D
CC BY-SA 4.0
2020-08-31T06:03:53.403Z
2020-08-31T06:03:53.403Z
Colin D
<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>

CC BY-SA 4.0
2020-08-29T02:21:28.870Z
Suppose that you could successfully creat an index file for your txt file, which program would be used to read it?
Phoenix Mu
CC BY-SA 4.0
2020-08-31T16:25:54.800Z
@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`
which_command
CC BY-SA 4.0
2020-09-01T17:13:59.130Z
I think tabix only works for bed and vcf format. Your file is not in bed format?
Phoenix Mu
CC BY-SA 4.0
2020-09-02T09:06:44.663Z
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
which_command
<parsing><tabix><gtex>
tabix indexing a GTEx variant file- ':' delimited
# 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**: @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`  
> **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 @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  
> **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  

---

