多次复制文件,写入重复文件,对文件进行排序,排序后计算特定行的位置

多次复制文件,写入重复文件,对文件进行排序,排序后计算特定行的位置

在详细介绍之前,我想指出我已经问过这个问题的某些部分 ->你可以在这里找到它。我收到了一些很好的答案,但我需要做更多的事情,所以这次我将重复我的问题并添加更多详细信息:

所以我有一个具有独特内容的文件,如下所示(让我们称之为myUniqueFile):

chromosoom  start    end       phylop   GPS
chr1    28745756    28745756    7.905   5   
chr1    31227215    31227215    10.263  5
chr1    47562402    47562402    2.322   4
chr1    64859630    64859630    1.714   3
chr1    70805699    70805699    1.913   2
chr1    89760653    89760653    -0.1    0
chr1    95630169    95630169    -1.651  -1

如您所见,这些都是具有不同分数的不同位置。

我还有另一个文件,如下所示(我们称其为myDuplicationFile

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5

因此,首先我想从 myUniqueFileto 添加行(标题除外),但我希望以与从 添加的每个新行重复的myDuplicationFile方式添加它们 。因此保留他的标准内容 + 添加 1 个新行。它应该看起来像这样:myDublicationFilemyUniqueFilemyDublicationFilemyUniqueFile

myDublicatedFile1.txt:

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    0.905   1    <- first line from `myUniquefile`


myDublicatedFile2.txt:

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5
chr1    31227215    31227215    10.263  5    <- second line from `myUniquefile`

因此,对于添加的每个新行,都会创建一个新文件,等等myDublicatedFile3,4,5

在我有了myDublicatedFiles新添加的内容后,我想对这些文件的特定列从高到低进行排序(对于菲洛普列)我这样做, for f in myDublicatedFile* ; do sort -g -r -k 3 $f >> $f.method1.txt所以看起来像这样:

myDuplicatedFile1.method1.txt:

chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    7.905   5 <- count 2
chr3    33093371    33093371    3.753   4
chr3    30648039    30648039    2.214   3
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    31663820    31663820    0.713   3
chr3    15540407    15540407    -1.391  -1
chromosoom  start    end       phylop   GPS

因此,在对这些文件进行排序后,我想知道排序后添加的行的位置。对我来说,用“grep”做一些事情并使用“count”似乎是合乎逻辑的。

因此,myDublicatedFile1.method1.txt这个计数/等级为 2,因为添加的行myUniquefile最终位于文件中的第二位。

计算出计数/排名后phlop(方法一)我想对列进行排序GPS(方法2)列,然后再次计算添加的行的排名。 myDuplicatedFile1.method1.method2.txt 应如下所示:

chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    7.905   5 
chr3    33093371    33093371    3.753   4
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    15540407    15540407    -1.391  -1
chromosoom  start    end       phylop   GPS

如果计数/排名写入不同的文件中,那么我可以稍后使用它们进行统计,这很容易。所以最重要的文件是这些计数,因为我最终会使用这些:)

就像是:

countsForMethod1.txt:

29
3
5
6
50
etc.

countsForMethod2.txt:

7
3
21
45
etc..

答案1

split假设您有GNU的版本,并且有像, or 这样的coreutilsshell可用(对于此处使用的进程替换功能),您可以修改之前接受的答案来处理标题行和排序,例如bashkshzsh

tail -n +2 myUniqueFile | SHELL=$(command -v bash) split -l1 --filter='{ 
  head -n 1 myDuplicationFile &&
    sort -g -r -k4,4 <(tail -n +2 myDuplicationFile) -
  } > "$FILE"'

然后,您可以使用简单的一行来查找输出文件中条目awk的位置:myUniqueFile

awk 'FNR==NR && NR>1 {a[$0]++; next} ($0 in a) {print FILENAME, FNR}' myUniqueFile xa?
xaa 3
xab 2
xac 4
xad 5
xae 5
xaf 8
xag 9

冲洗并重复其他方法/排序顺序。

答案2

该脚本无需创建临时文件即可计算排名(几乎会创建一个文件 - sorted_file)。此外,它对myDuplicationFile每种方法排序一次,然后进一步使用它。

#!/bin/bash

rank_determination() {
    # Sorts the "myDuplicationFile" one time
    # The "sorted_file" will be used further.
    ###
    tail -n +2 myDuplicationFile | sort -g -r -k "$1","$1" > sorted_file

    # gawk iterates through "myUniqueFile" line by line (except the first line).
    gawk -v field_number="$1" '
    NR != 1 {
        # Stores the needed value for the each line
        ###
        search_value=$field_number
        cnt=1

        # then, it checks the specified column in the "sorted_file"
        # line by line for the value, which is less than 
        # the "search_value" from the "myUniqueFile".
        ###
        while((getline < "sorted_file") > 0) {
            if($field_number < search_value)
                break
            cnt++
        }

        print cnt
        # closing is needed for reading the file from the beginning
        # each time. Else, "getline" will read line by line consistently.
        ###
        close("sorted_file")
    }' myUniqueFile
}

# I create a function, which takes
# the number argument, which means the column number:
# "4" for "phylop" column, "5" for the "GPS" column.
#
# The function creates output, which you can redirect
# to the needed file.
# Call this function multiple times with different arguments
# for the each needed column.
rank_determination 4 > method1.txt
rank_determination 5 > method2.txt

输出

tail -n +1 -- method*
==> method1.txt <==
2
1
3
4
4
7
8

==> method2.txt <==
2
2
3
5
6
7
8

答案3

我同意@WeijunZhou 在他的评论中所说的,你不需要创建所有这些临时文件来执行此操作。

以下 perl 脚本将一次性遍历两个文件,计算方法 1 (phylops) 和方法 2 (GPS) 排序的计数。

它的工作原理是保留重复文件中 phylop 和 GPS 值的排序列表(数组),然后(对于唯一文件的每一行)计算 phylop 和 GPS 值将排序到各自排序数组中的位置。

#!/usr/bin/perl

use strict;

# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';

# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
my @phylop=();
my @GPS=();

open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
  chomp;
  next if (m/^chromosoom/);

  my($chr,$start,$end,$phylop,$GPS) = split;
  push @phylop, $phylop + 0; # add 0 to make sure we only ever store a number
  push @GPS, $GPS + 0;
};
close(DUPE);

# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$a <=> $b} @phylop;
@GPS = sort {$a <=> $b} @GPS;

print "Method1\tMethod2\n";

# Now find out where the phylop and GPS value from each line of uniqfile
# would have ended up if we had sorted it into dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
  next if (m/^chromosoom/);
  chomp;

  my $phylop_sort_line=1;
  my $GPS_sort_line=1;

  my($chr,$start,$end,$phylop,$GPS) = split;

  for my $i (0..@phylop-1) {
    $phylop_sort_line++ if ($phylop < $phylop[$i]);
    $GPS_sort_line++ if ($GPS < $GPS[$i]);
  };

  #printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
  printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;  
};
close(UNIQ);

当针对您上面提供的示例数据运行时,输出为:

$ ./counts-for-methods.pl
Method1 Method2
2       1
1       1
3       2
4       3
4       5
7       7
8       7

该脚本完全忽略两个文件中的标题行,因此如果您当前的算法对这些行号进行计数,则这些行号可能会少一。

它还假设唯一文件中的值始终紧邻重复文件中的相同值排序。如果这不是您想要的,请将循环<中的比较更改for my $i (0..@phylop)<=

如果您分别需要方法 1 和方法 2 的值,您可以使用 轻松提取它们awk。或者perl可以轻松修改脚本以打开两个输出文件,每个文件对应一种方法,并将相应的值打印到每个文件。


这是处理输入行中 151 个字段的版本。我没有这样的输入文件,所以我使用代码中已注释掉的“5 个字段版本”对其进行了测试。输出与上面的版本相同。

#!/usr/bin/perl

use strict;

# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';

my @phylop=();
my @GPS=();

# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
  chomp;
  next if (m/^chromosoom/);

  my @fields = split;

# 151 fields version:
  push @phylop, $fields[42]+0;
  push @GPS, $fields[150]+0;

# 5 fields version:
#  push @phylop, $fields[3]+0;
#  push @GPS, $fields[4]+0;

};
close(DUPE);

# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$b <=> $a} @phylop;
@GPS = sort {$b <=> $a} @GPS;

print "Method1\tMethod2\n";

# Now find out where the phylop and GPS from each line of uniqfile
# would have ended up if we had sorted it into the dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
  next if (m/^chromosoom/);
  chomp;

  my $phylop_sort_line=1;
  my $GPS_sort_line=1;

  my @fields = split;

  for my $i (0..@phylop-1) {

# 151 fields version:
    $phylop_sort_line++ if ($fields[42] < $phylop[$i]);
    $GPS_sort_line++ if ($fields[150] < $GPS[$i]);

# 5 fields version:
#    $phylop_sort_line++ if ($fields[3] < $phylop[$i]);
#    $GPS_sort_line++ if ($fields[4] < $GPS[$i]);
  };

  #printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
  printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;

};
close(UNIQ);

答案4

只是为了节省一些打字,我将调用myUniqueFile samplemyDuplicationFile standard

#!/bin/bash                                                                     

(
while read line; do
  echo $line|cat standard -|tail -n +2|sort -g -r -k 4|awk '/^chr1/{print FNR}' >> countsForMethod1.txt
  echo $line|cat standard -|tail -n +2|sort -g -r -k 5|awk '/^chr1/{print FNR}' >> countsForMethod2.txt
done
) <(tail -n +2 sample)

解释:整个 while 循环由一对括号括起来,以便将bash其视为单个命令。该命令将文件sample作为输入,并使用该命令删除标题行tail。然后它被read命令一次一行地消耗。这意味着循环内部$line是文件中的一行sample。该变量被回显并通过管道cat生成myDuplicated*文件,只是它们是“动态”生成的并且从未写入磁盘。在tail对文件进行排序之前,标题行会被删除。awk然后用于找出样本位于哪一行。

编辑:我认为split有其优点,而这个答案消除了对中间文件的需要。

相关内容