☰

DNA序列计数

tag-icon • tag-icon text-processing bioinformatics
DNA序列计数

我有一段由一些空格分隔的 DNA 序列。我需要删除空格并返回不带空格字符的序列计数。关于如何做到这一点有任何帮助吗？我正在使用带有 bash 的终端。

例如，该序列如下所示：

GTCGATTGCAAGGATCCGCATGGGATAAAGGAATCGCAGTTCGAACAGGCAATGCCGCAG
CTATGATAGGACATCTCTTGGAGACACCTATTAATGTTTCAGAAACGGATACCTTGGTTG
TCCAGTACGAAATTAAGTTGGACAATTCTTTGACGTGCGGC
CTATATTAAAATTGTGGGTACATCACTCTCTTACCTGAGAATTCCAACAGAGCAGGACGC
TAACCCAGTGTCTATACCAGTCTGTGGCTTTGAAAGATTAGACACATTTCTGGATGAATT
TTCAAATTCTAAATTGATCGTTCAGTCTACACTAAGACATTCGTACGTTAGTCTTGAGAA

我想删除空白并准确计算有多少个碱基。或者，我们可以计算存在多少个 A、C、G 或 T，然后将它们相加，而不计算空格。

答案1

使用 GNU awk 进行多字符 RS 和 RT：

$ awk -v RS='[^\n]' 'RT{cnt[RT]++} END{for (base in cnt) print base, cnt[base]}' file
A 101
C 68
T 98
G 74

我假设您的描述中的“基础”是示例中不是换行符的任何字符。

答案2

假设没有空行或尾随空格等，您可以使用fold生成单个字母流，然后sort结合 来uniq -c计算每个字母有多少个：

$ fold -w 1 file | sort | uniq -c
 101 A
  68 C
  74 G
  98 T

如果输入中存在垃圾空白字符，则使用初始tr步骤删除这些字符：

$ tr -d -c 'ACGT' <file | fold -w 1 | sort | uniq -c
 101 A
  68 C
  74 G
  98 T

此处的命令tr将删除输入中除A、C、G或 之外的任何字符T。

可以用sort | uniq -c单个命令替换管道末尾的位，awk该命令计算输入中每​​个字符的出现次数，然后报告这些内容：

$ tr -d -c 'ACGT' <file | fold -w 1 | awk '{ count[$0]++ } END { for (ch in count) printf "%4d %s\n", count[ch], ch }'
 101 A
  68 C
  74 G
  98 T

但如果我们要引入awk，那么我们不妨去掉fold：

$ tr -d -c 'ACGT' <file | awk '{ for (i = 1; i <= length; ++i) count[substr($0,i,1)]++ } END { for (ch in count) printf "%4d %s\n", count[ch], ch }'
 101 A
  68 C
  74 G
  98 T

...也可能是tr：

$ awk '{ gsub("[^ACGT]", ""); for (i = 1; i <= length; ++i) count[substr($0,i,1)]++ } END { for (ch in count) printf "%4d %s\n", count[ch], ch }' file
 101 A
  68 C
  74 G
  98 T

代码awk，打印得很漂亮：

{
    gsub("[^ACGT]", "")  # removes anything not A, C, G, or T
    for (i = 1; i <= length; ++i)
        count[substr($0, i, 1)]++
}
END {
    for (ch in count) {
        printf "%4d %s\n", count[ch], ch
    }
}

第一个块（解析每一行输入）可以重写为使用gsub()而不是substr()：

{
    count["A"] += gsub("A", "A")
    count["C"] += gsub("C", "C")
    count["G"] += gsub("G", "G")
    count["T"] += gsub("T", "T")
}
END {
    for (ch in count) {
        printf "%4d %s\n", count[ch], ch
    }
}

...但是除了获得稍微少一点的嵌套代码之外，这可能不会比以前的代码有太大改进（除非它有助于某些用户的可读性）。

答案3

使用 Perl one 衬垫：

perl -F'' -e '
    BEGIN{my %h}
    map { /\S/ and $h{$_}++ } @F;
    END{print map { "$_ $h{$_}\n" } keys %h}
' file

输出

C 68
A 101
G 74
T 98

答案4

使用乐（以前称为 Perl_6）

raku -e '.say for slurp.comb(/\S/).Bag.pairs;'

示例输出：

G => 74
T => 98
A => 101
C => 68

或制表符分隔的输出（更改.say为.put）：

~$ raku -e '.put for slurp.comb(/\S/).Bag.pairs;' file
G   74
A   101
T   98
C   68

如果需要对输出进行排序，请添加.sort到末尾：

~$ raku -e '.put for slurp.comb(/\S/).Bag.pairs.sort;' file
A   101
C   68
G   74
T   98

或按最高核苷酸计数排序：

~$ raku -e '.put for slurp.comb(/\S/).Bag.pairs.sort: -*.value;' file
A   101
T   98
G   74
C   68

或者只是数一下（无字体空白）：

~$ raku -e '.put for slurp.comb(/\S/).elems;' file
341

最后，如果您正在处理非常大的文件，您可能需要尝试使用,lines.join来代替slurp, 以获得更好的内存管理。

输入示例：

GTCGATTGCAAGGATCCGCATGGGATAAAGGAATCGCAGTTCGAACAGGCAATGCCGCAG
CTATGATAGGACATCTCTTGGAGACACCTATTAATGTTTCAGAAACGGATACCTTGGTTG
TCCAGTACGAAATTAAGTTGGACAATTCTTTGACGTGCGGC
CTATATTAAAATTGTGGGTACATCACTCTCTTACCTGAGAATTCCAACAGAGCAGGACGC
TAACCCAGTGTCTATACCAGTCTGTGGCTTTGAAAGATTAGACACATTTCTGGATGAATT
TTCAAATTCTAAATTGATCGTTCAGTCTACACTAAGACATTCGTACGTTAGTCTTGAGAA

https://raku.org

相关内容