【问题标题】:compare text in two files and append text in a field比较两个文件中的文本并在字段中追加文本
【发布时间】:2012-12-13 13:48:03
【问题描述】:

我有两个文件: 文件 A 看起来像

ProbeID rsID    chr bp  strand  alleleA alleleB
SNP_A-1780270   rs987435    7   78599583    -   C   G
SNP_A-1780271   rs345783    15  33395779    -   C   G
SNP_A-1780272   rs955894    1   189807684   -   G   T
SNP_A-1780274   rs6088791   20  33907909    -   A   G
SNP_A-1780277   rs11180435  12  75664046    +   C   T
SNP_A-1780278   rs17571465  1   218890658   -   A   T
SNP_A-1780283   rs17011450  4   127630276   -   C   T
SNP_A-1780285   rs6919430   6   90919465    +   A   C
SNP_A-1780286   rs41528453  --- --- --- A   G
SNP_A-1780287   rs2342723   16  5748791 +   C   T

文件 B 看起来像

ProbeID call
SNP_A-1780270   2
SNP_A-1780271   0
SNP_A-1780272   2
SNP_A-1780274   1
SNP_A-1780277   0
SNP_A-1780278   2
SNP_A-1780283   2
SNP_A-1780285   2
SNP_A-1780286   0
SNP_A-1780287   0

我想要一个如下所示的输出:

ProbeID call    genotype
SNP_A-1780270   2   G G
SNP_A-1780271   0   C C
SNP_A-1780272   2   T T 
SNP_A-1780274   1   A G
SNP_A-1780277   0   C C
SNP_A-1780278   2   T T
SNP_A-1780283   2   T T 
SNP_A-1780285   2   C C
SNP_A-1780286   0   A A
SNP_A-1780287   0   C C

本质上,这与两个列表中的 ProbeID 匹配,并且在文件 B 中,检查调用列中相应的“调用”值。当 call = 0 时,在相邻列中打印 alleleA 的值两次。当 call = 1 时,打印 alleleA 和 alleleB 的值。当 call =2 时,两次打印 alleleB 的值。

【问题讨论】:

    标签: python r compare match


    【解决方案1】:

    使用pandas

    import pandas as pd
    import re
    
    A = pd.read_csv('FileA', delimiter = r'\s+')
    B = pd.read_csv('FileB', delimiter = r'\s+')
    A = A.set_index(['ProbeID'])
    B = B.set_index(['ProbeID'])
    C = pd.concat([A,B], axis = 1)
    
    idx = C['call'] == 0
    C['alleleB'][idx]  = C['alleleA'][idx]
    idx = C['call'] == 2
    C['alleleA'][idx]  = C['alleleB'][idx]
    print(C[['call', 'alleleA', 'alleleB']])
    

    产量

                   call alleleA alleleB
    ProbeID                            
    SNP_A-1780270     2       G       G
    SNP_A-1780271     0       C       C
    SNP_A-1780272     2       T       T
    SNP_A-1780274     1       A       G
    SNP_A-1780277     0       C       C
    SNP_A-1780278     2       T       T
    SNP_A-1780283     2       T       T
    SNP_A-1780285     2       C       C
    SNP_A-1780286     0       A       A
    SNP_A-1780287     0       C       C
    

    如果你有很多 Bfile,你可能会使用这样的东西:

    import pandas as pd
    import re
    
    A = pd.read_csv('FileA', delimiter = r'\s+')
    A = A.set_index(['ProbeID'])
    
    BFiles = ['FileB1', 'FileB2', 'FileB3']
    for i, bfile in enumerate(BFiles):
        B = pd.read_csv('FileB', delimiter = r'\s+')
        B = B.set_index(['ProbeID'])
        C = pd.concat([A,B], axis = 1)
    
        idx = C['call'] == 0
        C['alleleB'][idx]  = C['alleleA'][idx]
        idx = C['call'] == 2
        C['alleleA'][idx]  = C['alleleB'][idx]
        cfile = 'FileC{i}'.format(i = i)
        with open(cfile, 'w') as f:
            f.write(C[['call', 'alleleA', 'alleleB']])
    

    cfile 更改为合适的值。

    【讨论】:

    • 在某个时候,我应该学会使用这个神奇的pandas,我听说过很多...
    • @mgilson:我只是在自己学习。这很有趣,而且可能很有用。
    • @unutbu 当我将信息存储在文件中时,我尝试使用它,但无法获得任何成功。另外,如果有很多 FileB 我想要这样做并且每次都在不同的文件中输出,怎么办?
    • @jules:该文件使用可变数量的空格来分隔数据。幸运的是,pandas.read_csv 可以接受正则表达式模式作为分隔符。因此,您可以使用 pd.read_csv('FileA', delimiter = r'\s+') 将 FileA 读入 DataFrame。
    【解决方案2】:

    您可以使用嵌套字典轻松完成此操作:

    data = {}
    with open(fileA) as fA:
        header = next(fA).split()
        attributes = header[1:]
        for line in fA:
            lst = line.split()
            data[lst[0]] = dict(zip(attributes,l[1:])
    
    with open(fileB) as fB:
        header = next(fB).split()
        for line in fB:
            ID,call = line.split()
            data[ID]['call'] = int(call)
    

    现在您可以遍历您的数据并仅打印您需要的内容。

    或者,如果行完全对应,您可能一次使用itertools.izip 处理 1 行(或者如果使用 python 3,则只需简单的 zip):

    import itertools as it:
    
    with open(fileA) as fA,open(fileB) as fB:
        header_a = next(fA).split()
        header_b = next(fB).split()
        attrib_a = header_a[1:]
        attrib_b = header_b[1:]
        for line_a,line_b in it.izip(fA,fB):
            dat_a = line_a.split()
            dat_b = line_b.split()
            assert(dat_a[0] == dat_b[0])  #make sure they're the same ID
            dat = dict(zip(attrib_a,dat_a[1:]))
            dat.update(zip(attrib_b,dat_b[1:]))
            if (dat['call'] == '0'):
               print dat_a[0],dat['call'],dat['alleleA'],dat['alleleA']
    
            elif (dat['call'] == '1'):
               print dat_a[0],dat['call'],dat['alleleA'],dat['alleleB']
    
            elif (dat['call'] == '2'):
               print dat_a[0],dat['call'],dat['alleleB'],dat['alleleB']
    
            else:
                 raise AssertionError("Unknown call")
    

    【讨论】:

    • 当我有很多可以归类到fileB中的文件时,我尝试使用第二个。基本上我使用 import glob list = glob.glob('./*.txt') for file_name in list: 但是我得到一个错误 header_b = next(fB).split() StopIteration 我该如何解决这个问题?
    • @jules -- 据我所知,这意味着其中一个文件是完全空的。你可以把它放在try/except 子句中:docs.python.org/2/tutorial/errors.html#handling-exceptions——在这种情况下,你将使用except StopIteration: continue,它应该继续你的for 循环文件,而不是对文件做任何事情一个问题。
    【解决方案3】:

    这是一个 R 解决方案。

    my.data <- merge(df1, df2, by = "ProbeID")
    
    # select rows based on call
    zero <- my.data$call == 0
    one <- my.data$call == 1
    two <- my.data$call == 2
    
    # subset rows based on previous condition and calculate genotype
    my.data[zero, "genotype"] <- paste(my.data$alleleA[zero], my.data$alleleA[zero], sep = " ")
    my.data[one, "genotype"] <- paste(my.data$alleleA[one], my.data$alleleB[one], sep = " ")
    my.data[two, "genotype"] <- paste(my.data$alleleB[two], my.data$alleleB[two], sep = " ")
    
    my.data[, c("ProbeID", "call", "genotype")]
    
    
            ProbeID call genotype
    1  SNP_A-1780270    2      G G
    2  SNP_A-1780271    0      C C
    3  SNP_A-1780272    2      T T
    4  SNP_A-1780274    1      A G
    5  SNP_A-1780277    0      C C
    6  SNP_A-1780278    2      T T
    7  SNP_A-1780283    2      T T
    8  SNP_A-1780285    2      C C
    9  SNP_A-1780286    0      A A
    10 SNP_A-1780287    0      C C
    

    【讨论】:

      猜你喜欢
      • 2011-10-17
      • 1970-01-01
      • 1970-01-01
      • 2015-07-08
      • 1970-01-01
      • 1970-01-01
      • 2011-08-25
      • 2019-11-30
      • 2015-12-05
      相关资源
      最近更新 更多