Получите данные из огромного текстового файла, чтобы эффективно заменить данные в другом огромном текстовом файле (Python)

Question

Я программировал в течение нескольких месяцев, поэтому я не эксперт. У меня есть два огромных текстовых файла (omni, ~ 20 GB, ~ 2.5M строк, dbSNP, ~ 10 GB, ~ 60M строк). У них есть первые несколько строк, не обязательно с разделителями табуляции, начиная с «#» (заголовок), а остальные строки организованы в столбцах с разделителями табуляции (фактические данные).

Первые две колонки каждой строки содержат номер хромосомы и положение на хромосоме, а третий столбец содержит идентификационный код. В «omni» файле у меня нет идентификатора, поэтому мне нужно найти позицию в файле dbSNP (база данных) и создать копию первого файла, заполненного идентификаторами.

Из-за ограничений памяти я решил прочитать два файла по очереди и перезапустить их из последней строки. Меня не удовлетворяет эффективность моего кода, потому что я чувствую, что он медленнее, чем может быть. Я почти уверен, что это моя вина из-за отсутствия опыта. Есть ли способ сделать это быстрее с помощью Python? Может возникнуть проблема с открытием и закрытием файлов?

Я обычно запустить скрипт в GNOME Terminal (Python 2.7.6, Ubuntu 14,04), как это:

питон -u Replace_ID.py> Replace.log 2> Replace.err

Спасибо вам большое заблаговременно.

всенаправленный (Omni example):

...
#CHROM POS ID REF ALT ...
1 534247. КТ ...
...

dbSNP (dbSNP example):

...
#CHROM POS-ID Ref ALT ...
1 10019 rs376643643 ТА Т ...
...

выход должен быть точно таким же, как Omni файла, но с Rs ID после positi на.

Код:

SNPline = 0 #line in dbSNP file 
SNPline2 = 0 #temporary copy 
omniline = 0 #line in omni file 
line_offset = [] #beginnings of every line in dbSNP file (stackoverflow.com/a/620492) 
offset = 0 
with open("dbSNP_build_141.vcf") as dbSNP: #database 
    for line in dbSNP: 
     line_offset.append(offset) 
     offset += len(line) 
    dbSNP.seek(0) 
with open("Omni_replaced.vcf", "w") as outfile:  
    outfile.write("")  
with open("Omni25_genotypes_2141_samples.b37.v2.vcf") as omni: 
    for line in omni:   
     omniline += 1 
     print str(omniline) #log 
     if line[0] == "#":  #if line is header 
      with open("Omni_replaced.vcf", "a") as outfile: 
       outfile.write(line) #write as it is 
     else: 
      split_omni = line.split('\t') #tab-delimited columns 
      with open("dbSNP_build_141.vcf") as dbSNP: 
       SNPline2 = SNPline   #restart from last line found 
       dbSNP.seek(line_offset[SNPline])  
       for line in dbSNP: 
        SNPline2 = SNPline2 + 1 
        split_dbSNP = line.split('\t') 
        if line[0] == "#": 
         print str(omniline) + "#" + str(SNPline2) #keep track of what's happening. 
         rs_found = 0 #it does not contain the rs ID 
        else: 
         if split_omni[0] + split_omni[1] == split_dbSNP[0] + split_dbSNP[1]: #if chromosome and position match 
          print str(omniline) + "." + str(SNPline2) #log 
          SNPline = SNPline2 - 1 
          with open("Omni_replaced.vcf", "a") as outfile: 
           split_omni[2] = split_dbSNP[2] #replace the ID 
           outfile.write("\t".join(split_omni)) 
          rs_found = 1 #ID found 
          break   
         else: 
          rs_found = 0 #ID not found 
       if rs_found == 0: #if ID was not found in dbSNP, then: 
        with open("Omni_replaced.vcf", "a") as outfile: 
         outfile.write("\t".join(split_omni)) #keep the line unedited 
       else: #if ID was found: 
        pass #no need to do anything, line already written 
    print "End." 

Answer 1

вот мой вклад в вашу проблему. Прежде всего, вот что я понимаю о вашей проблеме, просто чтобы проверить, что я прав: У вас есть два файла, каждый из которых представляет собой файл с разделителями таблиц. Первый, dbSNP, содержит данные, третьи столбцы которых являются идентификаторами, соответствующими номеру хромосомы гена (столбец 1) и положению гена на хромосоме (колонка 2).

Задача состоит в том, чтобы взять omni-файл и заполнить столбец ID всеми значениями, поступающими из файла dbNSP (на основе номера хромосомы и положения гена).

Проблема возникает из-за размера файлов. Вы пытались сохранить в позиции файла каждую строку, чтобы сделать поиск и перейти непосредственно к хорошей строке, чтобы избежать поместить весь файл dbnsp в память. Поскольку этот подход недостаточно для вас из-за множественного открытия файла, вот что я предлагаю.

Разберите файл dbNSP один раз, чтобы сохранить только основную информацию, то есть пары (number,position):ID. Из вашего примера, что соответствует:

1 534247 rs201475892 
1 569624 rs6594035 
1 689186 rs374789455 

это соответствует менее чем на 10% от размера файла в срок памяти, поэтому, начиная с файла 20GB будет загружаться в память меньше, чем 2 Гб, это, вероятно, доступным (не знаю, какую загрузку вы пытались раньше).

Итак, вот мой код для этого. Не стесняйтесь спрашивать объяснение, как в отличие от вас, я использую Object Programming.

import argparse 

#description of this script 
__description__ = "This script parse a Database file in order to find the genes identifiers and provide them to a other file.[description to correct]\nTake the IDs from databaseFile and output the targetFile content enriched with IDs" 

# -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- 
#classes used to handle and operate on data 
class ChromosomeIdentifierIndex(): 
    def __init__(self): 
     self.chromosomes = {} 

    def register(self, chromosomeNumber, positionOnChromosome, identifier): 
     if not chromosomeNumber in self.chromosomes: 
      self.chromosomes[chromosomeNumber] = {} 

     self.chromosomes[chromosomeNumber][positionOnChromosome] = identifier 

    def __setitem__(self, ref, ID): 
     """ Allows to use alternative syntax to chrsIndex.register(number, position, id) : chrsIndex[number, position] = id """ 
     chromosomeNumber, positionOnChromosome = ref[0],ref[1] 
     self.register(chromosomeNumber, positionOnChromosome, ID) 

    def __getitem__(self, ref): 
     """ Allows to get IDs using the syntax: chromosomeIDindex[chromosomenumber,positionOnChromosome] """ 
     chromosomeNumber, positionOnChromosome = ref[0],ref[1] 
     try: 
      return self.chromosomes[chromosomeNumber][positionOnChromosome] 
     except: 
      return "." 

    def __repr__(self): 
     for chrs in self.chromosomes.keys(): 
      print "Chromosome : ", chrs 
      for position in self.chromosomes[chrs].keys(): 
       print "\t", position, "\t", self.chromosomes[chrs][position] 

class Chromosome(): 
    def __init__(self, string): 
     self.values = string.split("\t") 
     self.chrs  = self.values[0] 
     self.position = self.values[1] 
     self.ID  = self.values[2] 

    def __str__(self): 
     return "\t".join(self.values) 

    def setID(self, ID): 
     self.ID = ID 
     self.values[2] = ID 

class DefaultWritter(): 
    """ Use to print if no output is specified """ 
    def __init__(self): 
     pass 
    def write(self, string): 
     print string 
    def close(self): 
     pass 

# -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- -*- 
#The code executed when the scrip is called 
if __name__ == "__main__": 

    #initialisation 
    parser = argparse.ArgumentParser(description = __description__) 
    parser.add_argument("databaseFile" , help="A batch file that contains many informations, including the IDs.") 
    parser.add_argument("targetFile" , help="A file that contains informations, but miss the IDs.") 
    parser.add_argument("-o", "--output", help="The output file of the script. If no output is specified, the output will be printed on the screen.") 
    parser.add_argument("-l", "--logs" , help="The log file of the script. If no log file is specified, the logs will be printed on the screen.") 
    args = parser.parse_args() 

    output = None 
    if args.output == None: 
     output = DefaultWritter() 
    else: 
     output = open(args.output, 'w') 

    logger = None 
    if args.logs == None: 
     logger = DefaultWritter() 
    else: 
     logger = open(args.logs, 'w') 

    #start of the process 

    idIndex = ChromosomeIdentifierIndex() 

    #build index by reading the database file. 
    with open(args.databaseFile, 'r') as database: 
     for line in database: 
      if not line.startswith("#"): 
       chromosome = Chromosome(line) 
       idIndex[chromosome.chrs, chromosome.position] = chromosome.ID 

    #read the target, replace the ID and output the result 
    with open(args.targetFile, 'r') as target: 
     for line in target: 
      if not line.startswith("#"): 
       chromosome = Chromosome(line) 
       chromosome.setID(idIndex[chromosome.chrs, chromosome.position]) 
       output.write(str(chromosome)) 
      else: 
       output.write(line) 

    output.close() 
    logger.close()   

Основная идея разборе dbNSP файл один раз и собирает все идентификаторы в словаре. Затем прокручивайте строку за строкой и выводите результат.

Вы можете вызвать скрипт так:

python replace.py ./dbSNP_example.vcf ./Omni_example.vcf -o output.vcf 

argparse модуль я использовал и импорта для обработки параметров также предоставить автоматическую помощь, поэтому описание параметра доступно с

python replace.py -h 

или

python replace.py --help 

Я считаю, что его подход будет быстрее, чем ваш, поскольку я знаю d файл один раз и просто работать с ОЗУ позже, и я приглашаю вас протестировать его.

NB: Я не знаю, знакомы ли вы программе Object, поэтому я должен упомянуть, что здесь все классы находятся в одном файле для публикации при переполнении стека. В практическом примере использования хорошей практикой было бы поместить весь класс в отдельный файл, например «Chromosome.py», «ChromosomeIdentifierIndex.py» и «DefaultWritter.py», а затем импортировать их в «replace.py " файл.