如何比较一个文件中的多行并输出组合条目答案

【问题标题】：How to compare multiple lines in one file and output a combined entry如何比较一个文件中的多行并输出组合条目
【发布时间】：2014-02-21 17:20:52
【问题描述】：

我有一个显示四列的文件：

chr 开始结束脚本

像这样：

chrI    128980  129130  F53G12.5b  
chrI    132280  132430  F53G12.5c.2  
chrI    132280  132430  F53G12.5a  
chrI    132280  132430  F53G12.5b  
chrI    132280  132430  F53G12.5c.1  
chrI    133600  133750  F53G12.5c.2  
chrI    133600  133750  F53G12.5a  
chrI    133600  133750  F53G12.5b  
chrI    133600  133750  F53G12.5c.1  
chrI    136240  136390  F53G12.4  
chrI    139100  139250  F53G12.3  
chrI    163220  163370  F56C11.2a  
chrI    163220  163370  F56C11.2b  
chrI    173900  174050  F56C11.6a  
chrI    173900  174050  F56C11.6b  
chrI    173900  174050  F56C11.6c  
chrI    182240  182390  F56C11.3  
chrI    184080  184230  Y48G1BL.2a  
chrI    190720  190870  Y48G1BL.2a

并且许多区域（由 chr start end 描述）被重复，因为它们映射到超过 1 个转录本

例如：

chrI    133600  133750  F53G12.5c.2  
chrI    133600  133750  F53G12.5a  
chrI    133600  133750  F53G12.5b  
chrI    133600  133750  F53G12.5c.1

我想要的是一个代码，它采用具有相同列 1、2、3 的行并从中获取第 4 列的最短公共部分（在本例中为 F53G12.5）并输出一个精简条目，即：

chrI    133600  133750  F53G12.5

或者例如：

chrI    83280   83430   Y48G1C.10a  
chrI    90420   90570   Y48G1C.10b  
chrI    90420   90570   Y48G1C.10c  
chrI    90420   90570   Y48G1C.10a

应该给

 chrI    83280   83430   Y48G1C.10a  
 chrI    90420   90570   Y48G1C.10

您对此有什么建议吗？非常感谢

【问题讨论】：

对于我们非生物信息学家（！），您的意思是找到第 1,2 和第 3 列都相同的行，然后打印第 4 列中所有它们共有的最短部分？
是的！就是这样。感谢观看
有没有这样的日期chrI 83280 83430 Y48G1C.10a和chrI 90420 90570 Y48G1C.9a，你需要得到什么结果？

标签： python perl unix awk bioinformatics

【解决方案1】：

我怀疑这可以用 Pandas 完成，比这要好得多，但我对 Pandas 还不是很熟悉，所以……不调试就提交了。

def longest_identical_substring(words):
    result = words[0]
    for idx in range(len(words[0]), 0, -1):
        substrings = [w[:idx] for w in words]
        if max(substrings) == min(substrings): 
            result = substrings[0]
        else:
            return result

transcripts = defaultdict(list)
with open('myfile.csv') as infile:
    reader = csv.reader(infile)
    for row in reader:
        transcripts[row[:3]].append(row[3])
for ((chr, start, end), ts) in transcripts.items():
    print(chr, start, end, longest_identical_substring(ts))

【讨论】：

【解决方案2】：

awk 的一种方式。如果需要，您可以通过管道将其发送至 sort。

`script.awk`的内容

(a[$1" "$2" "$3]) {
    t=0; word=""; delete w1; delete w2;
    split($4,w1,""); 
    split(a[$1" "$2" "$3],w2,"");
    t=(length($4)<length(a[$1" "$2" "$3]))?length($4):length(a[$1" "$2" "$3])
    for (x=1;x<=t;x++) { 
        if (w1[x]==w2[x]) { 
            word=word""w1[x] 
        }
    a[$1" "$2" "$3]=word
    }
    next
} 

{
    a[$1" "$2" "$3]=$4
}

END {
        for (x in a)  print x,a[x]
}

您的文件：

$ cat file
chrI    128980  129130  F53G12.5b
chrI    132280  132430  F53G12.5c.2
chrI    132280  132430  F53G12.5a
chrI    132280  132430  F53G12.5b
chrI    132280  132430  F53G12.5c.1
chrI    133600  133750  F53G12.5c.2
chrI    133600  133750  F53G12.5a
chrI    133600  133750  F53G12.5b
chrI    133600  133750  F53G12.5c.1
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2a
chrI    163220  163370  F56C11.2b
chrI    173900  174050  F56C11.6a
chrI    173900  174050  F56C11.6b
chrI    173900  174050  F56C11.6c
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2a
chrI    190720  190870  Y48G1BL.2a

输出：

$ awk -f script.awk file
chrI 173900 174050 F56C11.6
chrI 128980 129130 F53G12.5b
chrI 182240 182390 F56C11.3
chrI 139100 139250 F53G12.3
chrI 136240 136390 F53G12.4
chrI 132280 132430 F53G12.5
chrI 163220 163370 F56C11.2
chrI 184080 184230 Y48G1BL.2a
chrI 190720 190870 Y48G1BL.2a
chrI 133600 133750 F53G12.5

【讨论】：

【解决方案3】：

这里没有所有的调试是一个简单的 awk 语句：

awk -F"." '{ trimmed=substr($2,RSTART,1);print $1"."trimmed;}' test.txt 
chrI    128980  129130  F53G12.5
chrI    132280  132430  F53G12.5
chrI    132280  132430  F53G12.5
chrI    132280  132430  F53G12.5
chrI    132280  132430  F53G12.5
chrI    133600  133750  F53G12.5
chrI    133600  133750  F53G12.5
chrI    133600  133750  F53G12.5
chrI    133600  133750  F53G12.5
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2
chrI    163220  163370  F56C11.2
chrI    173900  174050  F56C11.6
chrI    173900  174050  F56C11.6
chrI    173900  174050  F56C11.6
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2
chrI    190720  190870  Y48G1BL.2


awk -F"." '{ trimmed=substr($2,RSTART,1);print $1"."trimmed;}' test.txt |sort|uniq
chrI    128980  129130  F53G12.5
chrI    132280  132430  F53G12.5
chrI    133600  133750  F53G12.5
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2
chrI    173900  174050  F56C11.6
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2
chrI    190720  190870  Y48G1BL.2

按列排序 sort -nrk numeric reverse k for column id 在这种情况下我传递了 2 和 3

awk -F"." '{ trimmed=substr($2,RSTART,1);print $1"."trimmed;}' test.txt |sort -nrk2,3|uniq
chrI    190720  190870  Y48G1BL.2
chrI    184080  184230  Y48G1BL.2
chrI    182240  182390  F56C11.3
chrI    173900  174050  F56C11.6
chrI    163220  163370  F56C11.2
chrI    139100  139250  F53G12.3
chrI    136240  136390  F53G12.4
chrI    133600  133750  F53G12.5
chrI    132280  132430  F53G12.5
chrI    128980  129130  F53G12.5

根据列更新：

awk  '{ if( match($4, /[0-9a-zA-Z]+\.[0-9a-zA-Z]/)) {  trimmed=substr($4,RSTART,RLENGTH); } print $1"\t"$2"\t"$3"\t"trimmed;}' test.txt |sort|uniq
chrI    128980  129130  F53G12.5
chrI    132280  132430  F53G12.5
chrI    133600  133750  F53G12.5
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2
chrI    173900  174050  F56C11.6
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2
chrI    190720  190870  Y48G1BL.2

【讨论】：

抱歉，您似乎是通过修剪而不是通过查找最长的公共子字符串来完成的。因此，当我尝试在第 4 列的条目（如 Y48G1C.10b、Y48G1C.10c、Y48G1C.10a）上运行它时，它给了我 Y48G1C.1 而不是 Y48G1C.10
解决方案是根据提供的输出给出的。即 awk -F"。"是它基于hangon的地方，我将根据真实专栏对其进行更新
尝试最新的更新版本，它可能仍然会中断，因为正如我所说的那样，它是基于给出的内容，您的文件可能会完全不同

【解决方案4】：

我认为 Python 的 itertools.groupby 和 itertools.takewhile 可以处理您的问题的两个部分，按前三列中的值对行进行分组，并将第四列修剪为其公共前缀。

import itertools
from operator import itemgetter

def combine(data):
    for group, group_lines in itertools.groupby(data, itemgetter(0,1,2)):
        names = [line[3] for line in group_lines]
        prefix = "".join(t[0] for t in itertools.takewhile(lambda x:len(set(x))==1,
                                                           zip(*names)))
        yield group + (prefix,)

用类似的东西运行它：

with open(filename) as f:
    for item in combine(line.split() for line in f):
        print("{:8}{:8}{:8}{}".format(*item))

运行示例：

>>> data = """chrI    128980  129130  F53G12.5b
chrI    132280  132430  F53G12.5c.2
chrI    132280  132430  F53G12.5a
chrI    132280  132430  F53G12.5b
chrI    132280  132430  F53G12.5c.1
chrI    133600  133750  F53G12.5c.2
chrI    133600  133750  F53G12.5a
chrI    133600  133750  F53G12.5b
chrI    133600  133750  F53G12.5c.1
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2a
chrI    163220  163370  F56C11.2b
chrI    173900  174050  F56C11.6a
chrI    173900  174050  F56C11.6b
chrI    173900  174050  F56C11.6c
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2a
chrI    190720  190870  Y48G1BL.2a""".splitlines()
>>> for item in combine(line.split() for line in data):
        print("{:8}{:8}{:8}{}".format(*item))


chrI    128980  129130  F53G12.5b
chrI    132280  132430  F53G12.5
chrI    133600  133750  F53G12.5
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2
chrI    173900  174050  F56C11.6
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2a
chrI    190720  190870  Y48G1BL.2a

【讨论】：

【解决方案5】：

使用 awk

awk  '{sub(/[^0-9]+/,"",$2);NF=2} !a[$1]++' FS=. OFS=. file

chrI    128980  129130  F53G12.5
chrI    132280  132430  F53G12.5
chrI    133600  133750  F53G12.5
chrI    136240  136390  F53G12.4
chrI    139100  139250  F53G12.3
chrI    163220  163370  F56C11.2
chrI    173900  174050  F56C11.6
chrI    182240  182390  F56C11.3
chrI    184080  184230  Y48G1BL.2
chrI    190720  190870  Y48G1BL.2

【讨论】：

【解决方案6】：

这是昨晚的努力完成了:-)

#!/bin/bash
sort file                                           |\
   awk '
      NR==1 {f123=$1" "$2" "$3;trans=$4;next}       # NR=1, i.e. first line

      {                                             # NR!=1, i.e. subsequent lines
         if(f123!=$1" "$2" "$3){                    # Fields 1-3 have changed
            printf "%s %s\n",f123,trans
            f123=$1" "$2" "$3;trans=$4
         }else{                                     # Fields 1-3 unchanged, do common transcript
            newtrans=$4
            x=length(newtrans)                      # Get shorter of two transcripts
            if(length(trans)<x) x=length(trans)     # Copy common part
            common=""
            for(c=1;c<=x;c++){
               if(substr(trans,c,1)==substr(newtrans,c,1))common=common""substr(trans,c,1)
            }
            trans=common
         }
      }
      END {if(common)printf "%s %s\n",f123,common}
   '

一些注释...基本上输入文件已排序，因此具有相似 char/start/end 值的记录彼此相邻。然后将其通过管道传输到 awk。读取第一行时，字段（列）1 到 3 聚集在一起并保存为变量“f123”。在读取后续行时，将前 3 列与最后看到的 3 列进行比较。如果前三列的任何部分发生了变化，则看到的最后一行将连同其副本一起输出。如果前三列没有改变，那么我们有一个新的成绩单要处理。最后一个副本和当前副本共有的最短前缀然后通过复制字母来计算，直到一个不同，并保存新的副本以供下次第 1 到 3 列更改时输出。当我们达到最后一条记录时，我们可能正在积累一个新的通用成绩单，如果是，我们就输出它。

【讨论】：

【解决方案7】：

你们太棒了！感谢您提交的所有答案。这是使用 Perl 的另一种解决方案：

#!/usr/bin/env perl
#use Data::Dumper qw(Dumper);
use strict;
use warnings;


my $filename = $ARGV[0];
my @matrix;
my @transcripts;
my @transcript;
my %referenceTable;

my $count=0;
my $oldkey="";
my $key="";
my @keys;
my @key;
my %hash;
    open FILE,"< $filename" or die "can not open file\n";
    while (my $line=<FILE>) {
      my ($chromosome, $start, $stop, $transcript) = split("\t", $line);
        $key = $chromosome . "SPACE" . $start . "SPACE" . $stop;
        if ($oldkey ne $key) {

            $count = 0;
            $oldkey = $key;
        }
        push @{$referenceTable{$key}}, $transcript;

        $count++;
     }
my $output;
my ($k, $v); #Not @v -- $v will contain string that will be a reference to an array
while (($k, $v) = each(%referenceTable)){
 my ($chromosome, $start, $stop) = split(/SPACE/, $k);
 print "chromosome start stop \: $chromosome\t $start\t $stop \t";
 print "Common prefix \: \t ";
 $output = getleastcommonprefix(@{$v});
  print $output . "\n";
}

#print Dumper \%referenceTable;


sub getleastcommonprefix {
    my @searcharray = @_;
    my $common      = $searcharray[0];
    foreach my $index (1 .. $#searcharray) {
        $_ = $searcharray[0] . reverse $searcharray[$index];
        m/(.*)(.*)(??{quotemeta reverse $1})/s;
        if (length $1 < length $common) {
            $common = $1;
        }
    } ## end foreach my $index (1 .. $#searcharray)
    return $common;
} ## end sub getleastcommonprefix

#print 'Common prefix for file $filename [' . getleastcommonprefix(@array_of_test_names) . "]\n";

【讨论】：

script.awk的内容

您的文件：

输出：

`script.awk`的内容