1

実験データを 4 行のグループにまとめたテキスト ファイルがあります。まれなデータポイントを削除したいと思います。以下は最初の 8 行のファイル形式で、これが何千行も繰り返されます (行番号はファイルに存在しません)。

1 Info_1
2 Sequence_1
3 Info_1
4 Info_1
5 Info_2
6 Sequence_2
7 Info_2
8 Info_2
9 Info_3
10 Sequence_3
11 Info_3
12 Info_3

したがって、行 1 ～ 4 にはシーケンス 1 の情報が含まれ、行 5 ～ 8 にはシーケンス 2 の情報が含まれ、9 ～ 12 行にはシーケンス 3 の情報が含まれます。特定の状況では、完全に一意であるか、3 回未満しか検出されていないシーケンスを含む 4 行のグループを削除するのが一般的です。

私がやりたいのは、2行目を6、10、14、18行と比較することです...そしてそれが3回以上見つかった場合は何もしません。見つかった回数が 3 回以下の場合は、行 1 ～ 4 と、一致するシーケンスを含む 4 行の各グループを削除します。次に、ファイル内の 1 行おきに同じ比較を実行します。

したがって、上記のファイルのシーケンス 1 とシーケンス 3 が一致し、そのシーケンスが 3 回未満しか繰り返されていないため、4 行の各グループを削除すると、結果のファイルは次のようになります。

1 Info_2
2 Sequence_2
3 Info_2
4 Info_2

これが私が始めたものです：

awk 'FNR==NR {
    if (FNR%4==2) {
    a[$1]++
    if (a[$1]<3) b[int(FNR/4)]=1
    }
    next}
b[int(FNR/4)]==0' inputFile inputFile > outputFile

ただし、見つかった回数が 3 回未満のすべての行が削除されるわけではありません。助けていただければ幸いです。ありがとう。

要求された実際のテスト可能な例を次に示します。 入力:

Info1
AAGC
Info1
Info1
Info2
AACT
Info2
Info2
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4
Info5
AACT
Info5
Info5

AAGC回数が発生しますが、>= 3回数がAACT発生するため<3、出力は次のようになります。

Info1
AAGC
Info1
Info1
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4

うまくいけば、それが明確にするのに役立ちます。

4

3 に答える 3

2

この質問については、高水準のプログラミング言語と専門のバイオインフォマティクス ライブラリを使用することをお勧めします。私は python と biopython を使用します。

from Bio import SeqIO

input_file = open("input.fastq")

#to store sequence in dictionary using dna sequences as keys
#it consume RAM (if fastq is very very large, it could be a problem)
dict_seq = {}
for seq in SeqIO.parse(input_file, "fastq"):
  if not str(seq.seq) in dict_seq:
    dict_seq[str(seq.seq)] = []
  dict_seq[str(seq.seq)].append(seq)

#filter sequences
list_filter = [ dna for dna in dict_seq.keys() if len(dict_seq[dna]) >= 3]

#print output in fastq format
output_file = open("filter.fastq", "w")
for dna in list_filter:
  for seq in dict_seq[dna]:
    output_file.write(seq.format("fastq"))

output_file.close()

入力.fastq

@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGGG
+EAS54_6_R1_2_1_443_348
;;;;;;;;;;;9;7;;.7;393333
@EAS54_6_R1_2_1_413_789
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88
@EAS54_6_R1_2_1_413_329
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88

filter.fastq

@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88
@EAS54_6_R1_2_1_413_789
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88
@EAS54_6_R1_2_1_413_329
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;;88
于 2015-05-29T02:57:47.623 に答える
1 
$ cat tst.awk
{
    numRecs += (NR%4 == 1 ? 1 : 0)
    rec[numRecs,(NR-1)%4 + 1] = $0
    cnt[$0]++
}
END {
    for (recNr=1; recNr<=numRecs; recNr++) {
        if (cnt[rec[recNr,2]] > 2) {
            for (fldNr=1; fldNr<=4; fldNr++) {
                print rec[recNr,fldNr]
            }
        }
    }
}

$ awk -f tst.awk file
Info1
AAGC
Info1
Info1
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4
于 2015-05-29T04:27:53.093 に答える