【发布时间】：2017-12-21 17:27:15
【问题描述】：
我有 FASTQ 格式的 DNA 序列数据，每条记录采用 4 行格式：
@sequence-header-信息
顺序
+
质量分数
序列行中的每个字符在质量得分行中都有一个对应的字符。所有的序列都包含一个 8-9 碱基的条形码序列，然后是一个引物区域，表示生物序列的开始。在条形码前面，可能有也可能没有一些垃圾碱基。我需要根据条形码将序列（及其相关信息）分成单独的文件，并删除引物之前的所有内容。对于序列本身，我可以使用 grep 和 sed 执行此操作（a 和 b 是通过循环带有样本编号的条形码序列文件设置的变量；正则表达式命中引物）：
a=AAAGCG
b=1
cat DNA.fastq | grep -A2 -B1 -E "${a}""CCTACGGG[ACGT]{1}GGC[AT]{1}GCAG" | grep -v "^--$" | sed -E 's/.+(CCTACGGG[ACGT]{1}GGC[AT]{1}GCAG)/\1/' > sample_${b}.fastq
这将很好地处理这样的输入数据，我想从第 2 行的开头删除 AAAGCG，从第 6 行的开头删除 ATAAAGCG）：
@M00816:90:000000000-AE7TD:1:2116:19022:11483 1:N:0:1
AAAGCGCCTACGGGTGGCAGCAGTGGGGAATTTTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGTGGGAAGAAGGCCTTCGGGTTGTAAACCACTTTTGTCAGGGAAGAAACGGTCTGAACTAATATTTCGGACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCAGCTTTGCAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGCACTGCA
+
CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGGGGCFFGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGCAGGGG?FGGGGGGGFGGGG7FG7CGGGGGGGGGGGGEGGGFFFG585C98DEEG?CFG*<<@>37;68CEFCFGF99EGGEEDGE54<9C<8CC>*854C<.
@M00816:90:000000000-AE7TD:1:2118:10209:10682 1:N:0:1
ATAAAGCGCCTACGGGGGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTATTCAAGTCAGAGGTGCAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAACCTGGATAGCTTGAATCAT
+
--CCCCCGGGGGGGGGGGGGGGGGGFGGGGDEGGGGGGFGGGFGGGGGGGGFGGGG8AEA9DF@FGGGGGGEFGFGFFCFGGG,@FFFG<,<FFGGGGGGDGGGGGGGGFGG7BEE>EDGGGFGFGGFFGGGGGGGGGGGCFGGGC:CGFGGG8,<CEEGGGGFGGGFCE5CEEGGGGGGGGGGGG88C5@CEEE5CCFGGGGGGFGCFGGCG:>G>@FGDCG=F?*1+?E>EE@F@FFG9<9EGGGGG477AFFFDG69F=04C89?FGG7C6CFGG@:EECC8/;;EGC?EC>FF:DA74)
问题是如何从质量得分的开头删除正确数量的字符。在上面的示例中，我需要从第 4 行 (CCCCCG) 的开头删除 6 个字符，从第 8 行 (--CCCCCG) 删除 8 个字符。这可以通过计算第 2 行和第 4 行中条形码和引物之前的碱基数并从分数中删除该字符数来完成。或者，可以在输出文件上使用第二个命令来计算每个序列的长度（sample_1.fastq 中的第 2 行和第 4 行）并将质量得分修剪为相同的长度。两者的问题在于它们需要从一行存储信息并在后续行的进程中使用它。
想法1：将序列长度存储在单独的文本文件中
cat sample_1.fastq | sed -n '2~4p' | awk 'BEGIN{print length($0)} > seq_lengths.txt'
问题是我想不出如何在文件的每一行与 sample_1.fastq 中的每条记录同时循环（与嵌套循环相反）。
想法2：
cat sample_1.fastq | awk '$0 == "+" {
a=length $NR-1
b=gensub(/^.{a}/, 1, $NR+1)
print $(NR-2), $(NR-1), $0, b
}' > sample1test.fastq
测试表明长度 $NR-1 给出的长度是 -1，而不是前一行的长度！ （我还是个 awk 新手。）
想法3：统计sample_1.fastq中每个序列行的长度，并在下一行打印出来。再次搜索文件，grep -B2 -A2 这个数字并存储为变量。使用 sed 将质量得分修剪为变量的值并删除多余的行。这需要在不中断管道流动的情况下存储变量，据我了解这是不可能的。
思路4：将2的倍数设为var1，将4的倍数设为var2。使用 awk 打印 sample_1.fastq 行号 var1 的长度并将其设置为 var3。修剪线 var2 到长度 var3。这需要能够按号码呼叫线路，我知道这也是不可能的。
任何想法或见解将不胜感激！非常感谢。
【问题讨论】：
-
您可以通过在重定向之前将文件名放在 awk 语句的末尾来避免使用 cat 管道。还需要将长度用作函数，即 a=length($NR-1)
-
不清楚你想做什么。试着举一个例子说明你需要计算什么以及你必须从行中删除什么。
-
这种格式足够先进，我强烈建议使用 Python、Ruby 或 Java 等更好的语言来完成这项工作。
-
我在示例中添加了更多细节，希望能更清楚。