2017-04-05 67 views
1

我有2個文件。我如何知道位置是否在使用python的區域中

文件A有3列:染色體,起始位置,結束位置。

CHR,START,END 
chr1,1203245,1203374 
chr1,1202020,1202213 
chr1,1201293,1201465 
chr1,1200844,1201128 
chr1,1200527,1200585 

文件B有2列:染色體,位置。

CHR,POS 
chr1,1579264 
chr1,1641372 
chr1,3020521 
chr2,2097836 
chr3,2374462  

這兩個文件都很大。

如何決定文件B的每個位置在文件A的任何區域中使用python?如果只有一個位置和一個區域,我可以編寫代碼,但我不知道區域列表。

Position is in a region意味着CHR應該是相同的,POS >= START and POS <= END,例如,文件B的chr1,1203250應在文件chr1,1203245,1203374一個

我拿oscarbranson的建議,編寫代碼:

for i,r in B.iterrows(): 
    B.loc[i, 'in_A'] = any((r.CHR == A.CHR) & (r.POS >= A.SATRT) & (r.POS <= A.END)) 

但兩者的2文件很大,代碼仍在運行。如果有什麼辦法可以讓這個更快?

+0

您能通過編輯您的問題向我們展示任何輸入和輸出的例子嗎? –

+0

我猜正則表達式會做。爲了給你舉個例子,我們必須知道POS列和START/END之間的關係。你願意再次更新你的問題,解釋什麼是比賽? –

回答

0
  1. 如果你的文件是大,你不能只加載都在數據幀,因爲它會填滿整個內存。

  2. 不是比較B的位置與A中的每個區間,因爲這非常低效(明白:非常非常慢)。通過CHROM

    1. 排序兩個文件,開始,結束:

    相反,做到這一點。

  3. 同時對兩個文件的行進行迭代。像這樣的事情(爲理念;它沒有考慮染色體進去了,可能有錯誤的話,再加上你需要定義保存位置的方式):

    a = open("A") 
    b = open("B") 
    # init 
    (chromA,start,end) = a.readline().split(",") 
    (chromB,pos) = b.readline().split(",") 
    while 1: 
        try: 
         if pos < start: 
          # Read more lines from B 
          (chromB,pos) = b.readline().split(",") 
         elif pos > end: 
          # Read more lines from A 
          (chromA,start,end) = a.readline().split(",") 
         else: 
          # Pos is in the [start, end] interval 
          # Save the result and go to the next pos 
          savePos() 
          (chromB,pos) = b.readline().split(",") 
        except StopIteration: 
         break 
    

我真正的建議儘管不會像過去每個生物信息學家那樣手動完成此操作,但要使用「bedtools」或現有的Python庫(如BioPython)。重新研發輪子是一個可怕的時間和資源的研究損失。

-1

我把你的問題解釋爲'在A的所有行中的任何START和END值之間的B的每一行中的POS值?如果這是正確的,這樣的事情會的工作,使用pandas數據分析庫:

import pandas as pd 

A = pd.read_csv('path/to/A.csv') 
B = pd.read_csv('path/to/B.csv') 

for i, r in B.iterrows(): 
    B.loc[i, 'in_A'] = any((r.POS > A.START) & (r.POS < A.END)) 

print(B) 

>>  CHR  POS in_A 
    0 chr1 1579264 False 
    1 chr1 1641372 False 
    2 chr1 3020521 False 
    3 chr2 2097836 False 
    4 chr3 2374462 False 

這遍歷B中的每一行,並生成一個布爾值數組,它是TRUE其中(POS> START)&(POS < END),然後在B數據幀中創建一個新列,如果POS中存在任何的A開始/結束限制,則爲真。

有意義嗎?

+0

非常感謝您的幫助。我使用這樣的熊貓:'B.loc [i,'in_A'] = any((r.CHR == A.CHR)&(r.POS> = A.SATRT)&(r.POS <= A. END))',但它太慢了。任何建議,以使這更快? – kelloggs

+0

這很慢,因爲此答案將每個位置與其他文件中的每個間隔進行比較,因此它隨長度(A)x長度(B)而增長。如果你之前對文件進行排序,你可以想出更聰明的東西。這將兩個文件完全加載到內存中; OP表示他們「很大」,如果他們真的這樣做是行不通的。 – JulienD

相關問題