파이썬은 슬라이딩 윈도우 분석의 반복마다 스텝 크기를 줄입니다.

Question

저는 염색체를 가로 질러 일련의 슬라이딩 윈도우에서 통계치 Tajima의 D를 추정하는 프로그램을 만들고 있습니다. 염색체 자체는 기능적 중요성을 가진 여러 다른 영역으로 나뉘어져있다. 슬라이딩 윈도우 분석은 각 지역의 내 스크립트로 수행됩니다.

프로그램 시작시 슬라이딩 창의 크기와 한 창에서 다음 창으로 이동하는 단계의 크기를 정의합니다. 각 다른 염색체 지역에 대한 좌표를 포함하는 파일을 가져오고 작업중인 모든 SNP 데이터를 포함하는 다른 파일을 가져옵니다 (이것은 큰 파일이므로 줄 단위로 읽음). 프로그램은 염색체 위치 목록을 반복합니다. 각 위치에 대해 분석을위한 단계 및 창 인덱스를 생성하고 SNP 데이터를 출력 파일 (단계에 해당)로 분할하고 각 단계 파일의 주요 통계를 계산하며 이러한 통계를 결합하여 각 창에 대한 타지 마의 D를 추정합니다.

이 프로그램은 SNP 데이터의 작은 파일에 적합합니다. 또한 첫 번째 염색체 중단 점 이상의 첫 번째 반복에도 잘 작동합니다. 그러나 SNP 데이터의 대용량 파일의 경우 프로그램이 각 염색체 영역을 반복 할 때 분석 단계 크기가 알 수 없게 줄어 듭니다. 첫 번째 염색체 영역의 경우 단계 크기는 2500 개의 뉴클레오타이드 (이것은 예상되는 것입니다)입니다. 그러나 두 번째 염색체 세그먼트의 경우 단계 크기는 1966이고 세 번째 값은 732입니다.

이유가 궁금한 사람이 있으면 알려 주시기 바랍니다. 나는이 프로그램이 작은 파일에는 크기로 작동하지만 큰 파일에는 작동하지 않는 것처럼 특히 곤란하다.

내 코드는 다음과 같습니다 : 그래서,

stepStart = int(L[1]) 
    stepStop = int(L[2]) 
    stepSize = int(stepStop-(stepStart-1)) 

stepStop 및 stepStart 파일 '내용에 의존 나타납니다

import sys 
import math 
import fileinput 
import shlex 
import string 
windowSize = int(500) 
stepSize = int(250) 
n = int(50)  #number of individuals in the anaysis 
SNP_file = open("SNPs-1.txt",'r') 
SNP_file.readline() 
breakpoints = open("C:/Users/gwilymh/Desktop/Python/Breakpoint coordinates.txt", 'r') 
breakpoints = list(breakpoints) 
numSegments = len(breakpoints) 
# Open a file to store the Tajima's D results: 
outputFile = open("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/Tajima's D estimates.txt", 'a') 
outputFile.write(str("segmentNumber\tchrSegmentName\tsegmentStart\tsegmentStop\twindowNumber\twindowStart\twindowStop\tWindowSize\tnSNPs\tS\tD\n")) 

#Calculating parameters a1, a2, b1, b2, c1 and c2 
numPairwiseComparisons=n*((n-1)/2) 
b1=(n+1)/(3*(n-1)) 
b2=(2*(n**2+n+3))/(9*n*(n-1)) 
num=list(range(1,n))    # n-1 values as a list 
i=0 
a1=0 
for i in num: 
    a1=a1+(1/i) 
    i=i+1 
j=0 
a2=0 
for j in num: 
    a2=a2+(1/j**2) 
    j=j+1 
c1=(b1/a1)-(1/a1**2) 
c2=(1/(a1**2+a2))*(b2 - ((n+2)/(a1*n))+ (a2/a1**2)) 

counter6=0 
#For each segment, assign a number and identify the start and stop coodrinates and the segment name 
for counter6 in range(counter6,numSegments): 
    segment = shlex.shlex(breakpoints[counter6],posix = True) 
    segment.whitespace += '\t' 
    segment.whitespace_split = True 
    segment = list(segment) 
    segmentName = segment[0] 
    segmentNumber = int(counter6+1) 
    segmentStartPos = int(segment[1]) 
    segmentStopPos = int(segment[2]) 
    outputFile1 = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_Count of SNPs and mismatches per step.txt")%(str(segmentNumber),str(segmentName))), 'a') 

#Make output files to index the lcoations of each window within each segment 
    windowFileIndex = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_windowFileIndex.txt")%(str(segmentNumber),str(segmentName))), 'a') 
    k = segmentStartPos - 1 
    windowNumber = 0 
    while (k+1) <=segmentStopPos: 
     windowStart = k+1 
     windowNumber = windowNumber+1 
     windowStop = k + windowSize 
     if windowStop > segmentStopPos: 
      windowStop = segmentStopPos 
     windowFileIndex.write(("%s\t%s\t%s\n")%(str(windowNumber),str(windowStart),str(windowStop))) 
     k=k+stepSize 
    windowFileIndex.close() 

# Make output files for each step to export the corresponding SNP data into + an index of these output files 
    stepFileIndex = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_stepFileIndex.txt")%(str(segmentNumber),str(segmentName))), 'a') 
    i = segmentStartPos-1 
    stepNumber = 0 
    while (i+1) <= segmentStopPos: 
     stepStart = i+1 
     stepNumber = stepNumber+1 
     stepStop = i+stepSize 
     if stepStop > segmentStopPos: 
      stepStop = segmentStopPos 
     stepFile = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_step_%s.txt")%(str(segmentNumber),str(segmentName),str(stepNumber))), 'a') 
     stepFileIndex.write(("%s\t%s\t%s\n")%(str(stepNumber),str(stepStart),str(stepStop))) 
     i=i+stepSize 
    stepFile.close() 
    stepFileIndex.close() 

# Open the index file for each step in current chromosomal segment 
    stepFileIndex = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_stepFileIndex.txt")%(str(segmentNumber),str(segmentName))), 'r') 
    stepFileIndex = list(stepFileIndex) 
    numSteps = len(stepFileIndex) 

    while 1: 
     currentSNP = SNP_file.readline() 
     if not currentSNP: break 
     currentSNP = shlex.shlex(currentSNP,posix=True) 
     currentSNP.whitespace += '\t' 
     currentSNP.whitespace_split = True 
     currentSNP = list(currentSNP) 
     SNPlocation = int(currentSNP[0]) 
     if SNPlocation > segmentStopPos:break 
     stepIndexBin = int(((SNPlocation-segmentStartPos-1)/stepSize)+1) 
     #print(SNPlocation, stepIndexBin) 
     writeFile = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_step_%s.txt")%(str(segmentNumber),str(segmentName),str(stepIndexBin))), 'a') 
     writeFile.write((("%s\n")%(str(currentSNP[:])))) 
     writeFile.close() 

    counter3=0 
    for counter3 in range(counter3,numSteps): 
# open up each step in the list of steps across the chromosomal segment: 
     L=shlex.shlex(stepFileIndex[counter3],posix=True) 
     L.whitespace += '\t' 
     L.whitespace_split = True 
     L=list(L) 
     #print(L) 
     stepNumber = int(L[0]) 
     stepStart = int(L[1]) 
     stepStop = int(L[2]) 
     stepSize = int(stepStop-(stepStart-1)) 
#Now open the file of SNPs corresponding with the window in question and convert it into a list: 
     currentStepFile = open(("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_step_%s.txt")%(str(segmentNumber),str(segmentName),str(counter3+1)),'r') 
     currentStepFile = list(currentStepFile) 
     nSNPsInCurrentStepFile = len(currentStepFile) 
     print("number of SNPs in this step is:", nSNPsInCurrentStepFile) 
     #print(currentStepFile) 
     if nSNPsInCurrentStepFile == 0: 
      mismatchesPerSiteList = [0] 
     else: 
# For each line of the file, estimate the per site parameters relevent to Tajima's D 
      mismatchesPerSiteList = list() 
      counter4=0 
      for counter4 in range(counter4,nSNPsInCurrentStepFile): 
       CountA=0 
       CountG=0 
       CountC=0 
       CountT=0 
       x = counter4 
       lineOfData = currentStepFile[x] 
       counter5=0 
       for counter5 in range(0,len(lineOfData)): 
        if lineOfData[counter5]==("A" or "a"): CountA=CountA+1 
        elif lineOfData[counter5]==("G" or "g"): CountG=CountG+1 
        elif lineOfData[counter5]==("C" or "c"): CountC=CountC+1 
        elif lineOfData[counter5]==("T" or "t"): CountT=CountT+1 
        else: continue 
       AxG=CountA*CountG 
       AxC=CountA*CountC 
       AxT=CountA*CountT 
       GxC=CountG*CountC 
       GxT=CountG*CountT 
       CxT=CountC*CountT 
       NumberMismatches = AxG+AxC+AxT+GxC+GxT+CxT 
       mismatchesPerSiteList=mismatchesPerSiteList+[NumberMismatches] 
     outputFile1.write(str(("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n")%(segmentNumber, segmentName,stepNumber,stepStart,stepStop,stepSize,nSNPsInCurrentStepFile,sum(mismatchesPerSiteList)))) 
    outputFile1.close() 

    windowFileIndex = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_windowFileIndex.txt")%(str(segmentNumber),str(segmentName))), 'r') 
    windowFileIndex = list(windowFileIndex) 
    numberOfWindows = len(windowFileIndex) 
    stepData = open((("C:/Users/gwilymh/Desktop/Python/Sliding Window Analyses-2/%s_%s_Count of SNPs and mismatches per step.txt")%(str(segmentNumber),str(segmentName))), 'r') 
    stepData = list(stepData) 
    numberOfSteps = len(stepData) 

    counter = 0 
    for counter in range(counter, numberOfWindows): 
     window = shlex.shlex(windowFileIndex[counter], posix = True) 
     window.whitespace += "\t" 
     window.whitespace_split = True 
     window = list(window) 
     windowNumber = int(window[0]) 
     firstCoordinateInCurrentWindow = int(window[1]) 
     lastCoordinateInCurrentWindow = int(window[2]) 
     currentWindowSize = lastCoordinateInCurrentWindow - firstCoordinateInCurrentWindow +1 
     nSNPsInThisWindow = 0 
     nMismatchesInThisWindow = 0 

     counter2 = 0 
     for counter2 in range(counter2,numberOfSteps): 
      step = shlex.shlex(stepData[counter2], posix=True) 
      step.whitespace += "\t" 
      step.whitespace_split = True 
      step = list(step) 
      lastCoordinateInCurrentStep = int(step[4]) 
      if lastCoordinateInCurrentStep < firstCoordinateInCurrentWindow: continue 
      elif lastCoordinateInCurrentStep <= lastCoordinateInCurrentWindow: 
       nSNPsInThisStep = int(step[6]) 
       nMismatchesInThisStep = int(step[7]) 
       nSNPsInThisWindow = nSNPsInThisWindow + nSNPsInThisStep 
       nMismatchesInThisWindow = nMismatchesInThisWindow + nMismatchesInThisStep 
      elif lastCoordinateInCurrentStep > lastCoordinateInCurrentWindow: break 
     if nSNPsInThisWindow ==0 : 
      S = 0 
      D = 0 
     else: 
      S = nSNPsInThisWindow/currentWindowSize 
      pi = nMismatchesInThisWindow/(currentWindowSize*numPairwiseComparisons) 
      print(nSNPsInThisWindow,nMismatchesInThisWindow,currentWindowSize,S,pi) 
      D = (pi-(S/a1))/math.sqrt(c1*S + c2*S*(S-1/currentWindowSize)) 
     outputFile.write(str(("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n")%(segmentNumber,segmentName,segmentStartPos,segmentStopPos,windowNumber,firstCoordinateInCurrentWindow,lastCoordinateInCurrentWindow,currentWindowSize,nSNPsInThisWindow,S,D))) 

Answer 1

빠른 검색하면 110 라인에 stepSize을 변경 않음을 보여줍니다 우리는 더 이상 디버깅 할 수 없습니다.