☰

DNA序列的bash shell命令

tag-icon • tag-icon text-processing command-line
DNA序列的bash shell命令

我有包含 DNA 序列的文件，如下例所示。我怎样才能找到每个中A、T、G、C的计数柱子（是的，不是行，而是列）使用 bash shell 命令。

该文件有 846975 行，如下所示：

AAGAAAGAAGAGGAACTTCTCTCCATCCAGCCTCATTCCACTGCACCAACTCTTCTGTGT
CGGGTTGTGCAGGAGAGAAAGGGAGCTTGGCAACTCTTTGCTGTGCTGAGTTGTGGTAGC
CCATCACTGGGTTGTAAAGTGCCTTGCCTCCTTTCCTCCCCTCCTTTTTTTTTGAGACAG
AGTCTCACTCTGTCGTCCAGGCTGAGGTGCAGTGGTGCGATCTCTGCTCACTGCAACCTC
AGCCTCCTGGGTTCAAGTGATTCTCCTGCCTCAGCCTCCCAGGAAGCTGGGACTACAGGC
ACATGCCACCACACCTGGCTAACTTTTTTTTATTTTTAGTAGAGAAAGGGTATCACCATG
TTGGCCAGGCTGGTCTTGAACTCCTGACTTCAGGTGATCCACCCACCTTGGCCCCCCAAA
GTGCTGGGGTTAAAGGCATGAGACACTGCGCCCGTCCACCTCCTCTTTTACTTGGGAGAA
ATGCACAGATTCTGGGTGCCATGTGCATTTGTTTTGGGAGTGATAATTGATCTAACTTAT
GGAAATAATACTAGATAGTTAGCGGATGGATTCTGTATCTGATGAGAGTTTTGGGCAAAA
CGAATTCCTAGTTTCTGAGTCTTATTTTTCCCCTGATTCAAGAAAACTGTGAATTATCCA
GCCAGTAAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGACACTGGATTTCTACCACTC
ATTTTCTTACTACTTTTCCTGTGCAAGGATCATGGCACAAGTTGCAGTTTCCACCCTGCC
CATTGAAGATGAGGAGTCTGTTGAAGATGAGGAGTCCTTGGAGAGCAGGATGGTGGTGAC
ATTCCTGTCAGCTCTCGCCTCCATGGTCAGACCTTCTGTTCTCACATTCTGTAGTTCGGT
AGGACTGGGCGGTAGATAAGGTTGATTTGTTTTCGTAGAACTTACAATTTTGTGATTTTT
AGTTCTAATGAGTAGACCTTTTTCGTGAATAGTAGTTACGATCAAACACCTCTGACCAAA

对于此示例，第一列的输出应为A=9,T=1,G=3,C=4。

答案1

在一般情况下，如果您awk确实支持空/空字段分隔符（这是gawk特定的），您可以执行如下操作（当您想要每行计数时）。

gawk -F '' '{for (x=1;x<=NF;x++) arr[NR$x]++;
    for (i in arr) printf("%s:%d\t", i, arr[i]);delete arr;print ""}' infile

不，如果您确实需要每列的计数（假设字段分隔符为空），首先需要将输入数据从列转换为行，然后执行与上面相同的操作。

gawk -F '' '{for (x=1;x<=NF;x++) arr[NR$x]++;
    for (i in arr) printf("%s:%d\t", i, arr[i]);delete arr;print ""}' \
<(gawk -F '' '{for (i=1;i<=NF;i++) ctr[i]=ctr[i]$i} 
    END{for (i in ctr) print ctr[i]}' infile)

通过下面的awk脚本，我们将列转置为行。

gawk -F '' '{for (i=1;i<=NF;i++) ctr[i]=ctr[i]$i} 
    END{for (i in ctr) print ctr[i]}' infile

您可以使用下面的命令代替上面的命令进行转置。

sed 's/./& /g' in|datamash -t ' ' transpose |sed 's/ //g'

答案2

我想你的意思是查找每个中 A、T、G、C 的计数排。

GNUawk解决方案：

awk -v FPAT="[A-Z]" '{ 
           printf "A: %s, T: %s, G: %s, C: %s\n",
           gsub("A",""),gsub("T",""),gsub("G",""),gsub("C","") 
       }' file

输出：

A: 17, T: 15, G: 9, C: 19
A: 11, T: 16, G: 24, C: 9
A: 8, T: 23, G: 10, C: 19
A: 9, T: 16, G: 16, C: 19
A: 11, T: 13, G: 16, C: 20
A: 17, T: 20, G: 10, C: 13
A: 11, T: 14, G: 13, C: 22
A: 12, T: 15, G: 17, C: 16
A: 14, T: 23, G: 15, C: 8
A: 20, T: 18, G: 17, C: 5
A: 16, T: 23, G: 8, C: 13
A: 17, T: 16, G: 12, C: 15
A: 12, T: 21, G: 10, C: 17
A: 16, T: 14, G: 24, C: 6
A: 8, T: 22, G: 11, C: 19
A: 14, T: 25, G: 16, C: 5
A: 20, T: 18, G: 10, C: 12

答案3

我假设用户想要计算每行中的碱基数。

使用 Perl：

$ perl -ne 'chomp; printf("%s\t%2d %2d %2d %2d\n", $_, tr/A/A/, tr/C/C/, tr/G/G/, tr/T/T/)' file
AAGAAAGAAGAGGAACTTCTCTCCATCCAGCCTCATTCCACTGCACCAACTCTTCTGTGT    17 19  9 15
CGGGTTGTGCAGGAGAGAAAGGGAGCTTGGCAACTCTTTGCTGTGCTGAGTTGTGGTAGC    11  9 24 16
CCATCACTGGGTTGTAAAGTGCCTTGCCTCCTTTCCTCCCCTCCTTTTTTTTTGAGACAG     8 19 10 23
AGTCTCACTCTGTCGTCCAGGCTGAGGTGCAGTGGTGCGATCTCTGCTCACTGCAACCTC     9 19 16 16
AGCCTCCTGGGTTCAAGTGATTCTCCTGCCTCAGCCTCCCAGGAAGCTGGGACTACAGGC    11 20 16 13
ACATGCCACCACACCTGGCTAACTTTTTTTTATTTTTAGTAGAGAAAGGGTATCACCATG    17 13 10 20
TTGGCCAGGCTGGTCTTGAACTCCTGACTTCAGGTGATCCACCCACCTTGGCCCCCCAAA    11 22 13 14
GTGCTGGGGTTAAAGGCATGAGACACTGCGCCCGTCCACCTCCTCTTTTACTTGGGAGAA    12 16 17 15
ATGCACAGATTCTGGGTGCCATGTGCATTTGTTTTGGGAGTGATAATTGATCTAACTTAT    14  8 15 23
GGAAATAATACTAGATAGTTAGCGGATGGATTCTGTATCTGATGAGAGTTTTGGGCAAAA    20  5 17 18
CGAATTCCTAGTTTCTGAGTCTTATTTTTCCCCTGATTCAAGAAAACTGTGAATTATCCA    16 13  8 23
GCCAGTAAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGACACTGGATTTCTACCACTC    17 15 12 16
ATTTTCTTACTACTTTTCCTGTGCAAGGATCATGGCACAAGTTGCAGTTTCCACCCTGCC    12 17 10 21
CATTGAAGATGAGGAGTCTGTTGAAGATGAGGAGTCCTTGGAGAGCAGGATGGTGGTGAC    16  6 24 14
ATTCCTGTCAGCTCTCGCCTCCATGGTCAGACCTTCTGTTCTCACATTCTGTAGTTCGGT     8 19 11 22
AGGACTGGGCGGTAGATAAGGTTGATTTGTTTTCGTAGAACTTACAATTTTGTGATTTTT    14  5 16 25
AGTTCTAATGAGTAGACCTTTTTCGTGAATAGTAGTTACGATCAAACACCTCTGACCAAA    20 12 10 18

或者

$ perl -ne 'BEGIN { print(" A  C  G  T\n") } printf("%2d %2d %2d %2d\n", tr/A/A/, tr/C/C/, tr/G/G/, tr/T/T/)' file
 A  C  G  T
17 19  9 15
11  9 24 16
 8 19 10 23
 9 19 16 16
11 20 16 13
17 13 10 20
11 22 13 14
12 16 17 15
14  8 15 23
20  5 17 18
16 13  8 23
17 15 12 16
12 17 10 21
16  6 24 14
 8 19 11 22
14  5 16 25
20 12 10 18

Perl 中的运算符进行音译，与 shell 中的实用程序tr非常相似。tr它返回音译的字符数。

相关内容