☰

ファイル内の長さの異なる2つの連続した行を解析するスクリプト

tag-icon • tag-icon text-processing awk sed bioinformatics
ファイル内の長さの異なる2つの連続した行を解析するスクリプト

連続する2行の長さが同じ（テキストが完全に異なる）大容量ファイルを解析しようとしています。検索してみましたが、最初の記事がここにありますね。スクリプトを見つけて修正してみましたが、楽しかったです。 file はソート出力ファイルです。シーケンスと品質スコアを分析することで、ファイルは次のようになります。

CCTCGNAACCCAAAAACTTTGATTTCTNATAAGGTGCCAGCGGAGTCCTAAAAGCAACATCCGCTGATCCCTGGT
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
CCCCANCCAAACTCCCCACCTGACAATNTCCTCCGCCCGGATCGACCCGCCGAAGCGAGTCTTGGGTCTAAA
AAAAA#EEEEEEEEEEEAEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
ATCGTNTATGGTTGAGACTAGGACGGTNTCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAAAAC
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEAEEEEEEAEEEAEEEEEEEEEEEEEEEEEEEEEEE
CCCACNTGGAGCTCTCGATTCCGTGGGNTGGCTCAACAAAGCAGCCACCCCGTCCTACCTATTTAAAGTTTG
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE
GCATCNTTTATGGTTGAGACTAGGACGNTATCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAA
6AAAA#EEEEEAAAEEEEEEAEEAEEE#EEEEEEEAEAEEEEAEEAAA/EAEEEEAEEAEEAEEAEAAEEEEEE

質問：各シーケンスベースには、スコアのない破損した線のペアがあります。つまり、各ペアの2行の長さは同じでなければなりません。間違った行のペアをどのように解析しますか？ファイルには1億行があります。

私はparser.shというコードを試しました。

{ curr = $0 }
(NR%2)==0 {
    currLgth = length(curr)
    prevLgth = length(prev)
    maxLgth = (currLgth > prevLgth ? currLgth : prevLgth)
    if (prevLgth==currLgth) {
        print ""
        print prevLgth
        print currLgth
        for (i=1; i<=maxLgth; i++) {
        }
    }
}
{ prev = curr }

実行されますが、awk -f parser.sh filename 「等しくない」（'=='）を使用しても、すべての行の長さが印刷されます。

75
75

72
72

75
75

72
72

私はコーダーではないので、事前に謝罪し、助けが必要です。通常、コードを見つけて修正して動作させることは可能ですが、この場合はそうではありません。 -血

Fastqファイルは一度に4行を読み取ります。 Read#1 e,g には次の 4 行が含まれます。

@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE

最初の行はサンプル名、2行目は実際のシーケンス、3行目は「+」記号、4行目はシーケンスの各塩基のASCII「スコア」セットです。ベースごとに1つのスコアしかないため、2行目の長さは4行目の長さと同じにする必要があります。私は2行と4行を分析し、長さの異なるペアを見つけました。代わりに、ペアリングが欠落しているように見える結果が表示されます。

以下は、疑問符が欠落しているか解決されていない品質スコアを示すFASTQファイルの例です。

@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

これは私の（ライン2 + 4）解析されたファイルの外観です：

CGGCATCGTTTATGGTTGAGACTAGGACG
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

間に品質スコアラインがない2つの連続シーケンス行があります。

ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT

あなたが私に与えたコードを使って：

awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}{last=length($0)}' Fastq-seq-qual-parsed.txt
Bad pair at lines 5 and 6

または: ./new-try.awk

答え1

私が提案する

awk '
    { first = $0; getline; second = $0 }
    length(first) != length(second) {
        print "Error at line", NR-1
        print first
        print second
    }
' file

通常のbashを使用することもできますが、速度がはるかに遅くなります。

nr=1
while IFS= read -r first; IFS= read -r second; do 
    if (( ${#first} != ${#second} )); then 
        printf "%s\n" "problem at line $nr" "$first" "$second"
    fi
    ((nr+=2))
done < file

答え2

努力する：

awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file

はい

これをテストファイルとして試してみましょう。

$ cat file
good123
good345
bad12
bad123
good_again
good_also1

コマンドを使用すると、一致しないペアが正しく識別されます。

$ awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file
Bad pair at lines 3 and 4

どのように動作しますか？

  • NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}

    偶数行にある場合は、NR%2==0その行の長さが前の行と同じであることを確認してください。そうでない場合は、length($0)!=lastメッセージを印刷します。

  • last=length($0)

    これにより、現在の行の長さが変数に保存されますlast。

複数行バージョン

コードを複数行にわたって分散したい場合：

awk '
    NR%2==0 && length($0)!=last {
        print "Bad pair at lines",NR-1,"and",NR
    }

    {
        last=length($0)
    }' file

ファイルの特定の行を印刷する方法

たとえば、ファイルの3行を印刷するには、次のようにします。

$ awk 'NR==3' file
bad12

範囲を印刷するには、3から6までのすべての行を印刷するには、次のようにします。

$ awk 'NR>=3 && NR<=6' file
bad12
bad123
good_again
good_also1

あるいは、sedを使用して同様の結果を得ることができます。

$ sed -n '3p' file
bad12
$ sed -n '3,6p' file
bad12
bad123
good_again
good_also1

フィルタリングされていない入力データの使用

次の入力ファイルを検討してください。

$ cat File
@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@sample7
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE

?次のように、誤ったサンプル、つまり長さが異なる行があるか、2行目で始まるサンプルを検出できます。

$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(/^\?/ || length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample4 (line 11) is bad:
ATTTCGGGGGGGGGGGGGG
??????????????????????????????????
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE

または、2行目（「quality」）が次に始まるサンプルを無視したい場合?：

$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(!/^\?/ && length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE

答え3

まず、何かを見つけるために、5行と6行の長さが等しくないテストファイルを作成します（」cccccc「次のように）：

printf '%s\n' aaa aaa bbbb bbbb cccc ccc ddd ddd > foo

抽象的な金持ち2つの仮想ファイルに分割して使用するbash プロセスの交換そしてsed、各文字は次のように置き換えられます.。

  • これ最初仮想ファイルは物理ファイルを抽象化します。
  • これ2位仮想ファイルは抽象的です。奇妙な1行をコピーしてコピーします。2位連続変換奇妙なそしてでも線の長さは同じです。

...diff次のファイルは次のとおりです。

diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo)

出力には、6行が一致しないとマークされています。

6c6
< ...
---
> ....

上記の出力があまりにも冗長である場合は、diff同様のプログラムに多くのオプションがあるか、必要に応じてフィルタリングすることができます。行番号のみを表示：

diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo) | 
sed -n 's/c.*//p'

出力：

6

またはもう少し詳しく言うと、つまり数値が一致しない元のファイル行：

f=foo
diff <(sed 's/././g' $f) <(sed -n '1~2{s/././g;p;p}' $f) |  
sed -n 's/^\(.*\)c.*/\1/p' | grep -B 1 -wf - <(cat -n $f)

出力：

     5  cccc
     6  ccc

関連情報