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

ファイル内の長さの異なる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

おすすめ記事