2014-11-18 2 views
0

저는 몇 달 동안 프로그래밍을 해왔으므로 전문가는 아닙니다. 두 개의 거대한 텍스트 파일 (옴니 ~ 20GB, ~ 2.5M 라인, dbSNP, ~ 10GB, ~ 60M 라인)이 있습니다. 그것들은 처음 몇 줄을 가지고 있는데, 탭으로 구분 된 것은 아니지만 "#"(머리글)로 시작하고 나머지 줄은 탭으로 구분 된 열 (실제 데이터)로 구성됩니다.거대한 텍스트 파일의 데이터를 가져 와서 다른 거대한 텍스트 파일의 데이터를 효율적으로 대체하십시오 (Python)

각 줄의 처음 두 열은 염색체 번호와 염색체상의 위치를 ​​포함하며 세 번째 열은 식별 코드를 포함합니다. "옴니"파일에는 ID가 없으므로 dbSNP 파일 (데이터베이스)에서 위치를 찾고 ID로 완료된 첫 번째 파일의 복사본을 만들어야합니다.

메모리 제한으로 인해 두 파일을 한 줄씩 읽고 마지막 줄에서 다시 읽기로 결정했습니다. 코드의 효율성에 만족하지 않습니다. 느릴 수 있기 때문입니다. 나는 경험이 부족하기 때문에 내 잘못이라고 확신한다. 파이썬을 사용하여 더 빠르게 만들 수있는 방법이 있습니까? 문제는 파일을 열거 나 닫을 수 있습니까?

나는 보통 그놈 터미널이 같은 (파이썬 2.7.6, 14.04 우분투)에서 스크립트를 실행

:

파이썬 -u Replace_ID.py> Replace.log 2>

Replace.err

미리 감사드립니다. 무

(Omni example)

...
#CHROM POS ID REF ALT ...
1 534,247. CT ...
...

dbSNP (dbSNP example)

...
#CHROM POS ID REF의 ALT ...
1 10019 rs376643643 TA T ... 출력은 있지만 positi 후의 RS의 ID와, 옴니 파일과 완전히 동일해야


... 에.

코드 :

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." 
+0

테스트를 위해 자세한 정보를 제공해 주시겠습니까? 예를 들어 각 파일에 5 줄을 제공하고 교체 프로세스가 완료되었음을 나타내는 5 줄 출력을 제공합니까? 아마 당신이 제공하는 코드를 더 잘 이해할 수있을 것이며 또한 약간의 수정을 테스트하고 일관성을 검증 할 수있을 것입니다. –

+0

@ ArthurVaïsese 귀하의 제안에 감사드립니다. 원본 메시지를 편집했습니다. 두 개의 링크 만 추가 할 수 있으므로 다음은 출력 파일 예제입니다. lucabetti.altervista.org/file/Output_example.vcf – 4Kinesis

답변

0

여기 문제에 대한 나의 기여. 우선, 여기에 내가 당신의 문제에 대해 이해하고있는 것, 나는 정확하다는 것을 확인한다 : 두 개의 파일이 있는데, 각각은 표로 구분 된 값 파일이다. 첫 번째 dbSNP는 데이터를 포함하며 세 번째 열은 유전자의 염색체 번호 (열 1)와 염색체상의 유전자 위치 (열 2)에 해당하는 식별자입니다.

옴니 파일을 취하고 dbNSP 파일 (염색체 번호와 유전자 위치를 기반으로 함)에서 오는 모든 값으로 ID 열을 채우는 작업으로 구성됩니다.

문제는 파일 크기에서 비롯된 것입니다. 모든 dbnsp 파일을 메모리에 두지 않으려면 각 행의 파일 위치에서 찾기를 수행하고 양호한 행으로 바로 이동하려고했습니다. 이 방법은 여러 파일 열기로 인해 사용자에게 충분하지 않으므로 여기에 제가 제안한 내용이 있습니다.

dbNSP 파일을 한 번 구문 분석하여 필수 정보 (예 : 커플 (number,position):ID) 만 유지합니다. 에 해당하는 귀하의 예제에서 : 당신이 2GB 이하 적은 메모리에로드 때문에 20GB의 파일에서 시작, 메모리의 용어에서 파일 크기의 10 % 미만에

1 534247 rs201475892 
1 569624 rs6594035 
1 689186 rs374789455 

이 대응은, 아마 저렴 (전에 시도한로드 중 어떤 종류인지는 알 수 없습니다.)

그래서 여기에 제 코드가 있습니다. 당신과 달리 설명을 요청하는 것을 망설이지 말고, 나는 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 파일을 한 번 구문 분석하고 사전에있는 모든 ID를 수집하는 것입니다. 그런 다음 omnifile을 한 행씩 읽고 결과를 출력하십시오.

당신은이 같은 스크립트를 호출 할 수 있습니다

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

내가 사용하고 가져 오기 매개 변수의 설명을

python replace.py -h 

또는 함께 사용할 수 있도록 매개 변수는 자동 도움을 제출 처리하는 argparse 모듈

python replace.py --help 

저는 자신의 접근 방식이 d 파일을 한 번 실행하고 나중에 RAM에서 작업하고 테스트 해 보도록 권유합니다.

NB : 당신이 Object programming에 익숙한 지 모르겠다. 그래서 여기서는 모든 클래스가 스택 오버 플로우에 게시하기 위해 같은 파일에 있다고 언급해야한다. 실제 사용 사례에서 모든 클래스를 "Chromosome.py", "ChromosomeIdentifierIndex.py"및 "DefaultWritter.py"와 같은 별도의 파일에 넣은 다음 "replace.py"클래스로 가져 오는 것이 좋습니다. "파일.

+0

해답을 가져 주셔서 감사합니다. 나는 그것을 테스트했고 실제로 훨씬 더 빠릅니다. 또한 Python에서 Object 프로그래밍에 대해 많은 것을 배웠습니다. 이것이 가장 중요한 것입니다. 안녕 – 4Kinesis

관련 문제