如何将脏的fastq文件排序为交错的fastq

编程问答 2022-06-04

问题描述

我有一个大小约为80GB的fastq文件（file.fastq），其中有一个标题行和三个后续的信息行。我需要在匹配标题信息的标题行中匹配/ 1和/ 2，并将它们放在一个按/ 1和/ 2（交错）排序的文件中。
例如，我想为@HS2000-1015_160:5:2114:14081:69256/2找到@HS2000-1015_160:5:2114:14081:69256/1，如果找到了，我想将它们堆叠在一起（也对其他标头重复此操作），并得到如下所示的结果。

我当时在想这样的事情，但不知道该怎么做 cat file.fastq | grep -A3 --no-group-separator ...

file.fastq 

@HS2000-1015_160:5:1114:14809:86636/2
GTTTGGAAAACAGAAGGGATTCGTCTGGAAAGGGCGGGAAACAAAGGACTGGAAGGGAGTGATGGAAAAAAACATGATGAAACTGATCTGGAAAGAGCCA
+
@B@FFFFFHGHHHIIJJJHIJJJJIIJJJIJJJIJJIHFDDEDDBddddDBdddddddd2?CddddDCCddddddddCCD@CDDDCDDEDDCddddC9<A
@HS2000-1015_160:5:2114:14081:69256/1
CACCCTGATCCAAGGTGGACCTCATATGCTATCCACCAATTCTTTCATTAATCTTTATCTGAGTTCCTATTACATATTTCAAAGAAACTGTCAAACAACT
+
CCCFFFFFHHHHHJJCGIIHIIJJJIJJJIIJJJJIJJ;FHIJIJJIGIGDIIJJJJJIJJIIFGIJIEIJJIJJIIIJJJHHFGHFFFFDFDEDEEDDD
@HS2000-1015_160:7:1108:13370:100570/1
GTGGGTCTTTGATCTTACAGGACAAATCCTTCAGAAAGATCAGGAAATCAAGGACCTAAAACAGAAGATAGCCGAAGTCATGGCCGTCATGCCCAGCATA
+
CCCFFFFFHHHFHGGIJIJDHEIIIGEHGGII>CFII<BG<DGCADDHBD>>BC?(8CGCCHIEC;@9@D47?BAD4,.(.6.;(3=8<A>@:<@AB9AA
@HS2000-1015_160:7:1108:13370:100570/2
CTTGACTGCCAGAGACGCTCCTTTGCAATGCCTTCCGGTAACCAAATTTTTGGGCACAACACACAGCTGGCCTTCATTTCTTCAGGGGCTGGTAAACAGA
+
@@@ADFFFHHHFD=EF@:GHIIFHH<ECHGF@DDBB:6@D60?F=888)8='--(=5@EAE5?'(..((.;?@>>A>3;@####################
@HS2000-1015_160:6:1306:14563:31825/1
TTTTAAACTTATGAACATCCAACAGCAACTAATGGTAAGTCAGAAAGCTTTTATTTTTAAGTATATGCATGTGTATGTGTGTGTTTTGTATGTGTGTGTC
+
CCCFFFFFHHGHHJEDHJIJJJJIIJJJJJJGGIICFHIDGHH@HGIIIGGHFHHG>GGDHEBFIJJJHGI@FCDEEEHHGHHEHHFFFEFFFEEEEEDD
@HS2000-1015_160:7:2308:13476:86835/2
ATGGTCCTAAGTAACGTGAGGCCCAGAGGGAGGAGACCTTCAGAGGAGGCCAGGATGGGGGTGATGGGGAACTTTTCCTGAAGCTGGAGTAGGGAAAGAT
+
@@CFDFFFHFHDHHGIGJGIJIJJEG@FHGCFGBF=DFGFHGG>GGGGHIHCHFCEHHFFD'8<ACdddd38CCddddCDDACDCDC?B:>@C?>A?A>>
@HS2000-1015_160:5:2306:10070:71746/1
AGAAGACGAGCACAGTGAGGAAGGAGATGAGCAGGCAGGGGATGATGAGGTTGATGGTGTAGAACAAGGGCAGGCGCCGGATGTACAGCGAGTATGTGAT
+
CCCFFFFFHHHHHIICGGHHHAFGEFGIJJJIIJIIGIIGG;FGHEHHGI=CHCEHGHE?DFFDFECEddddddddDBDDBD@ACDAACDBD5<CDC>@>
@HS2000-1015_160:5:2306:10070:71746/2
GAACCTCAAGGACTATTGGGAGAGCGGCGAGTGGGCCATCATCAAAGCCCCAGGCTACAAACACGACATCAAGTACAACTGCTGCGAGGAGATCTACCCC
+
@CCFFFFDHGHFHIJJJJJJGGGIIJJIGHI@FHGIIGHHEFGHHFFFFFBCddddDCdddddddd;@BDCCDACDD@>ACCddddBDB<BA?C@CC@BD
@HS2000-1015_160:6:2116:4077:79041/1
TTGGGAGGGACATATCTAGGGTGGACCAAGCGCTGGCCAGGACTCACCTGGGCTGAAGAGGTGAGCCCAGAGGAGCTGGTGGTCTTGAAGCAGCTCCGCC
+
CCCFFFFFHHHHHIJIIIJFGCHIIIJJJIJJIJIJJJJJJJJJJJJJJJGHHHFFDDFC=ACEDDDB=?BD@ABDBDDCCBDCCCCCCCCCDDCdddd?
@HS2000-1015_160:6:2116:4077:79041/2
GGTCCCCGCCTACGCCCACTGGGTTGGTGCACCTGGTGGTGGTGGCCGCCAAGAAGCTGGTGAACCGCCTCCAAGTGGCTCCCAAGACGCAGCTGGATGA
+
CCCFFFFFHHFHHJHJJJJJJJJGGHJHIGAGIIJIFHJ;@F;CHHFHFddddDCDDCDD9CCCDDBDDBBDDCDACDD8@BD3>?BCDBDDDACCDC@>
@HS2000-1015_160:5:1110:16577:55703/2
AGAAGCTGAAGTACGGTGTAGGAATGTTTCTACCACGCATTTTCCCTGCACCCTCTCACATCATTCATGGGCCATTTGTAGACAGATGGGAAAGAGATGT
+
CC@FFFFFHGDDHGIG:C<FBHEFHIGGGIIJIIIJIFJJJJIIHIJIHHIIDHIIIJJJJJJJIHHGEHHGEFFFCEDECDDEE>@CDDCCDBCDCCDC
@HS2000-1015_160:5:2113:11446:94436/1
AATTGCTTGAACCCGGGAGGCAGAGGTTGCAGTGAGCCAAGATCGCACCACTGCAGTCCAGCCTGGGAAACCAGAGCAAAACTTGGTCTCAAAAAAAAAA
+
CCCFFFFFHHHHGHGFGFHjibF1DC<@;)9D*099?F;8=(BFF5(;'.=CEHFBDC7.>;?C@BBB35>(<?<<A>AC?B>4+(+>A:(:@(2)5@B#
@HS2000-1015_160:5:2113:11446:94436/2
CGTCAGGGCCAACCCCGCCCCACCCTGACCCTACCTGGCACCCCTCACCTGTGGCCTGCCAGCACAGCCTCGCCCCTGCTGGCCAATGTGTCCCCCGTCA
+
?@@DA@DDFHH?DHI)<@@FHDBGGCHCBDH;DFA<)6.=7D;@CBCHD)).7@=>;?==AABC95<(5(5309@D########################
@HS2000-1015_160:5:1216:14609:76938/2
TCAAAAGGCCTCTGGAGAGATTTCTGGGTGTGGGAGGTGCTATTTTGGGGCTTGGCCAGTCATATGGAGATAAGCCTACAAGGTTGGGGACCTGGCAGAT
+
CCCFFFFFFHHGHJJHHCCAFHIGHHJJ:E?FHHHDG?FGGIIJJIJIIJIJJIJHGCEECC@CBEFECEC@C@@CCddddDD>CCBDDD@??CBDC72:
@HS2000-1015_160:6:2209:18284:44195/1
GCAGCCAGCGGACGTCCAGGAACCGGGATGCCTCCAGCAGTGAGGCGGTCAGCCTGCAGCATGGGATGGCTGTGGATCTTTGGGGCAGCCCTGGGGCAGT
+
CCCFFFFFHHHHHIJHIIJJJJJJIJJEHIJJGIJJJJGHCHEHGFddddDCddddDDCDDCDDDCCBDCDDBCddddddddDBBDBddddDBddddBDC
@HS2000-1015_160:6:2209:18284:44195/2
TAAAATGTCACAAAGCTGGAAACTCTTCCCTATCACAAACCAAAACTTAAAAGGACGTTACCTGGCTGGGTCTAAACTCCACATAACTCGCTTGCAGTTG
+
CCCFFFFEHHHGHJIIIJJIJJHIIJEHJJHIJJJIIJJIJIJJIJIIHJJIJGGHGHGIIHHIIIIHFH@DFFFDEEEECDDDCddddBDDBBDCDACC

我想要的结果：

@HS2000-1015_160:7:1108:13370:100570/1
GTGGGTCTTTGATCTTACAGGACAAATCCTTCAGAAAGATCAGGAAATCAAGGACCTAAAACAGAAGATAGCCGAAGTCATGGCCGTCATGCCCAGCATA
+
CCCFFFFFHHHFHGGIJIJDHEIIIGEHGGII>CFII<BG<DGCADDHBD>>BC?(8CGCCHIEC;@9@D47?BAD4,.(.6.;(3=8<A>@:<@AB9AA
@HS2000-1015_160:7:1108:13370:100570/2
CTTGACTGCCAGAGACGCTCCTTTGCAATGCCTTCCGGTAACCAAATTTTTGGGCACAACACACAGCTGGCCTTCATTTCTTCAGGGGCTGGTAAACAGA
+
@@@ADFFFHHHFD=EF@:GHIIFHH<ECHGF@DDBB:6@D60?F=888)8='--(=5@EAE5?'(..((.;?@>>A>3;@####################
+
@HS2000-1015_160:5:2306:10070:71746/1
AGAAGACGAGCACAGTGAGGAAGGAGATGAGCAGGCAGGGGATGATGAGGTTGATGGTGTAGAACAAGGGCAGGCGCCGGATGTACAGCGAGTATGTGAT
+
CCCFFFFFHHHHHIICGGHHHAFGEFGIJJJIIJIIGIIGG;FGHEHHGI=CHCEHGHE?DFFDFECEddddddddDBDDBD@ACDAACDBD5<CDC>@>
@HS2000-1015_160:5:2306:10070:71746/2
GAACCTCAAGGACTATTGGGAGAGCGGCGAGTGGGCCATCATCAAAGCCCCAGGCTACAAACACGACATCAAGTACAACTGCTGCGAGGAGATCTACCCC
+
@CCFFFFDHGHFHIJJJJJJGGGIIJJIGHI@FHGIIGHHEFGHHFFFFFBCddddDCdddddddd;@BDCCDACDD@>ACCddddBDB<BA?C@CC@BD
@HS2000-1015_160:6:2116:4077:79041/1
TTGGGAGGGACATATCTAGGGTGGACCAAGCGCTGGCCAGGACTCACCTGGGCTGAAGAGGTGAGCCCAGAGGAGCTGGTGGTCTTGAAGCAGCTCCGCC
+
CCCFFFFFHHHHHIJIIIJFGCHIIIJJJIJJIJIJJJJJJJJJJJJJJJGHHHFFDDFC=ACEDDDB=?BD@ABDBDDCCBDCCCCCCCCCDDCdddd?
@HS2000-1015_160:6:2116:4077:79041/2
GGTCCCCGCCTACGCCCACTGGGTTGGTGCACCTGGTGGTGGTGGCCGCCAAGAAGCTGGTGAACCGCCTCCAAGTGGCTCCCAAGACGCAGCTGGATGA
+
CCCFFFFFHHFHHJHJJJJJJJJGGHJHIGAGIIJIFHJ;@F;CHHFHFddddDCDDCDD9CCCDDBDDBBDDCDACDD8@BD3>?BCDBDDDACCDC@>
+
@HS2000-1015_160:5:2113:11446:94436/1
AATTGCTTGAACCCGGGAGGCAGAGGTTGCAGTGAGCCAAGATCGCACCACTGCAGTCCAGCCTGGGAAACCAGAGCAAAACTTGGTCTCAAAAAAAAAA
+
CCCFFFFFHHHHGHGFGFHjibF1DC<@;)9D*099?F;8=(BFF5(;'.=CEHFBDC7.>;?C@BBB35>(<?<<A>AC?B>4+(+>A:(:@(2)5@B#
@HS2000-1015_160:5:2113:11446:94436/2
CGTCAGGGCCAACCCCGCCCCACCCTGACCCTACCTGGCACCCCTCACCTGTGGCCTGCCAGCACAGCCTCGCCCCTGCTGGCCAATGTGTCCCCCGTCA
+
?@@DA@DDFHH?DHI)<@@FHDBGGCHCBDH;DFA<)6.=7D;@CBCHD)).7@=>;?==AABC95<(5(5309@D########################
+
@HS2000-1015_160:6:2209:18284:44195/1
GCAGCCAGCGGACGTCCAGGAACCGGGATGCCTCCAGCAGTGAGGCGGTCAGCCTGCAGCATGGGATGGCTGTGGATCTTTGGGGCAGCCCTGGGGCAGT
+
CCCFFFFFHHHHHIJHIIJJJJJJIJJEHIJJGIJJJJGHCHEHGFddddDCddddDDCDDCDDDCCBDCDDBCddddddddDBBDBddddDBddddBDC
@HS2000-1015_160:6:2209:18284:44195/2
TAAAATGTCACAAAGCTGGAAACTCTTCCCTATCACAAACCAAAACTTAAAAGGACGTTACCTGGCTGGGTCTAAACTCCACATAACTCGCTTGCAGTTG
+
CCCFFFFEHHHGHJIIIJJIJJHIIJEHJJHIJJJIIJJIJIJJIJIIHJJIJGGHGHGIIHHIIIIHFH@DFFFDEEEECDDDCddddBDDBBDCDACC

解决方法

假设/断言：

  • 所有组都有4行（即没有3行或5行组）
  • 文件大小（80 GB）不允许将整个内容加载到内存（例如awk数组）中以进行分组/排序
  • 文件中不存在~字符
  • 每个唯一的标头最多匹配文件中的2行，并且所说的匹配具有后缀/1和/2

一个awk/sort/awk想法：

  • 使用~字符将每组4行附加在一起（有效地将\n替换为~）
  • 使用通用sort对附加行进行排序
  • 使用awk仅删除具有匹配的/1和/2条目的行
  • 将~字符替换为换行符（\n）

注意：fastq.dat包含问题中提供的“之前”数据的副本

awk '{ (FNR % 4 == 0) ? sfx="\n" : sfx="~"; printf "%s%s",$0,sfx }' fastq.dat | 
sort | 
awk -F '/' 'prevhdr != $1 { prevhdr=$1 ; prevline=$0 ; next } { gsub(/~/,"\n",prevline) ; print prevline ; gsub(/~/,$0) ; print }' > fastq.final

注释：

  • 显然将需要足够的空间来容纳另外80 GB的文件（也许在不同的文件系统上，以便通过2个在同一基础设备上读写80 GB的进程来最大程度地减少FS抖动吗？）
  • 假设此时操作系统可以处理机上内存需求，否则在写入单个输出文件（???）时可能需要分别执行每个步骤

OP已评论了内存使用情况和输出的缓慢生成。可能会感兴趣的是，看看将进程分解为单独的命令是否可能有所帮助，尽管是以需要更多磁盘空间为代价的，例如：

awk '{ (FNR % 4 == 0) ? sfx="\n" : sfx="~"; printf "%s%s",sfx }' fastq.dat > fastq.appended

sort fastq.appended > fastq.sorted 

awk -F '/' 'prevhdr != $1 { prevhdr=$1 ; prevline=$0 ; next } { gsub(/~/,$0) ; print }' fastq.sorted > fastq.final

'rm' -rf fastq.appended fastq.sorted 2>/dev/null

假设sort操作可能是大问题（内存和速度较慢），那么此answer about speeding up sort可能会很有趣。总体思路是增加sort使用的内存并允许使用并行线程，例如：

sort -S 50% --parallel=4

显示已维护4行组：

$ head -8 fastq.final
@HS2000-1015_160:5:2113:11446:94436/1
AATTGCTTGAACCCGGGAGGCAGAGGTTGCAGTGAGCCAAGATCGCACCACTGCAGTCCAGCCTGGGAAACCAGAGCAAAACTTGGTCTCAAAAAAAAAA
+
CCCFFFFFHHHHGHGFGFHJIBF1DC<@;)9D*099?F;8=(BFF5(;'.=CEHFBDC7.>;?C@BBB35>(<?<<A>AC?B>4+(+>A:(:@(2)5@B#
@HS2000-1015_160:5:2113:11446:94436/2
CGTCAGGGCCAACCCCGCCCCACCCTGACCCTACCTGGCACCCCTCACCTGTGGCCTGCCAGCACAGCCTCGCCCCTGCTGGCCAATGTGTCCCCCGTCA
+
?@@DA@DDFHH?DHI)<@@FHDBGGCHCBDH;DFA<)6.=7D;@CBCHD)).7@=>;?==AABC95<(5(5309@D########################

仅拉@HS行以显示数据已根据需要进行了排序：

$ grep '@HS' fastq.final
@HS2000-1015_160:5:2113:11446:94436/1
@HS2000-1015_160:5:2113:11446:94436/2

@HS2000-1015_160:5:2306:10070:71746/1
@HS2000-1015_160:5:2306:10070:71746/2

@HS2000-1015_160:6:2116:4077:79041/1
@HS2000-1015_160:6:2116:4077:79041/2

@HS2000-1015_160:6:2209:18284:44195/1
@HS2000-1015_160:6:2209:18284:44195/2

@HS2000-1015_160:7:1108:13370:100570/1
@HS2000-1015_160:7:1108:13370:100570/2

注意：仅出于显示目的而手动添加的空行

,

发表评论的时间有点长...想知道以下两种方法中的哪一种可能更适合于处理这些较大的数据集。

注意：我不使用fastq数据，以下内容基于大约15分钟的谷歌搜索，因此请稍候...

1-以允许更快处理的方式生成fastq文件

fastq-dump的文档中提到使用--split-3标志将把fastq数据转储到3x文件中……一个“左”文件，一个“右”文件和一个“单个”文件

“左”和“右”文件分别是指您的/1和/2后缀吗？如果为“是”，那么我想知道是否可以生成3x文件（然后丢弃“单个”文件）？

也将感兴趣的是查看“左”和“右”文件中的数据的标题是否具有相同的顺序（在这种情况下，2x文件之间的某种“合并”联接可能会比手动进行预排序然后匹配）。

假设可以生成3x文件，这引出了第二个想法...

2-是否有用于匹配左/右（/1 & /2）数据集的工具？

fastq-pair的简要介绍似乎表明，多个文件（例如fastq-dump生成的“左”文件和“右”文件-??上方）可以比许多蛮力更快地匹配。解决方案（例如，引用的awk解决方案与我的第一个答案有点相似）。

awkbashfastqlinuxlinuxshell

相关问答

Selenium Web驱动程序和Java元素在(x，y)点处不可单击其他元素将获得点击?
Selenium Web驱动程序和Java。元素在(x，y)点处不可单击。其...
Python-如何使用点“” 访问字典成员？
Python-如何使用点“。” 访问字典成员？
Java 字符串是不可变的到底是什么意思？
Java 字符串是不可变的。到底是什么意思？
Java中的“ final”关键字如何工作？我仍然可以修改对象
Java中的“ final”关键字如何工作？（我仍然可以修改对象。...
“loop:”在Java代码中这是什么，为什么要编译？
“loop:”在Java代码中。这是什么，为什么要编译？
java.lang.ClassNotFoundException：sun.jdbc.odbc.JdbcOdbcDriver发生异常为什么？
java.lang.ClassNotFoundException：sun.jdbc.odbc.JdbcOdbc...