DNA配列に対するbashシェルコマンド

DNA配列に対するbashシェルコマンド

以下の例に示すように、DNA配列を含むファイルがあります。それぞれのA、T、G、Cの数をどのように見つけることができますか?(はい、行ではなく列)bashシェルコマンドを使用します。

ファイルには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

通常、awknull/null フィールド区切り文字 (特定gawk) をサポートしている場合は、次のことができます (1 行あたりの数が必要な場合)。

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の数を求めてください。

牛に似た一種の栄養アッ解決策:

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 -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

trPerlの演算子はシェルのユーティリティと同様に翻訳されます。tr翻訳された文字数を返します。

関連情報