tabix indiziert eine GTEx-Variantendatei - ':' getrennt

Question

Ich habe eine gtex-Variantendatei, deren Kopf wie folgt aussieht:

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

Effektiv möchte ich in der Lage sein, die Varianten innerhalb eines bestimmten Gens nachzuschlagen (die Datei ist sortiert) und sie in eine temporäre Datei zu legen: ZB Wenn die Varianten im Gen "ENSG00000148481- MINDY3" in der Zeilennummer 87528225 bis 87536766- effektiv sind, was Ich möchte, ist das Äquivalent zu zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt. Ist zgrep ENSG00000148481 file.gzaber genauso schnell wie oben ...

Daher dachte ich, Tabix wäre das richtige Werkzeug dafür.

Ich möchte es tabix indizieren, um die Suche zu beschleunigen. Es ist gzip komprimiert und ich muss zuerst tun:

zcat gtex.txt.gz | bgzip > gtex.txt.bgz

Ich bin mir jedoch nicht ganz sicher, wie ich von dort aus vorgehen soll, da die Daten nicht durch Tabulatoren getrennt sind.

Als Versuch habe ich die ersten 1000 Zeilen ausprobiert:

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

Es gibt nichts zurück.

Ich fange an zu denken, dass ich nur ein Skript schreiben muss, das sich in das ':' aufteilt und die Daten in eine neue Datei schreibt ... und dann diese Datei indiziert - was sehr viel Zeit in Anspruch nimmt ... Kennt jemand einen Trick, um die Dateien mit unkonventioneller Begrenzung zu indizieren?

Answer 1

Tabix-Dateien benötigen mindestens einen Namen und eine Genomkoordinate. Tabix kann (ab) verwendet werden, indem eine gefälschte Genomkoordinate erstellt wird. Anschließend werden nur Namen indiziert, aber ich bin mir nicht sicher, ob ich diese Verwendung empfehlen würde. Was ist Ihre tatsächlich erwartete Anwendung davon?

Answer 2

Ich ging wie folgt vor:

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)

Das obige Skript schreibt einen Index für jedes Gen, der wie folgt aussehen sollte (es dauerte 2 Tage, bis alle Dateien fertig waren):

ensid    row_start    row_end
......
ENSG00000230337     1526391 1534329
ENSG00000171819 1534330 1542321
ENSG00000198793    1542322 1550178
ENSG00000120942  1550179 1558010
.....

Und um ein interessierendes Gen zu erhalten, müssen Sie die Koordinaten für ein Ensemblid / Gen nachschlagen und zum Beispiel:

zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' > ENSG00000120942.txt