☰

解析文件中长度不等的连续两行的脚本

tag-icon • tag-icon text-processing awk sed bioinformatics
解析文件中长度不等的连续两行的脚本

我正在尝试解析一个大文件，其中每两连续行具有相同的长度（文本完全不同）。我已经搜索过了，我的第一篇文章在这里。我找到了一个脚本并尝试修改它，但没有任何乐趣。文件是排序输出文件。我已经解析出序列和质量分数，因此文件如下所示：

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

问题：某处有一对损坏的线，使得每个序列碱基没有相应的分数，即每对两条线的长度应该相等，我如何解析出不正确的线对？文件有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

我不是编码员，所以提前道歉，但需要帮助。通常可以找到代码并修改它以使其工作，但在这种情况下不行。 -p

Fastq 文件每次读取有四行。 Read#1 e,g 将包含以下 4 行：

@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE

第一行是样本名称，第二行是实际序列，第三行是“+”符号，第四行是序列中每个碱基的一组 ASCII“分数”。每个碱基只有一个分数，因此第 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

有连续的两条序列线，它们之间没有质量得分线：

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

我们可以检测坏样本，即行长度不等或以 开头的第二行的样本?，如下所示：

$ 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

或者，如果我们想忽略第二行（“质量”）以 开头的样本?，则：

$ 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

抽象的富分成两个虚拟文件使用bash 流程替代和sed，其中每个字符都替换为.：

  • 这第一的虚拟文件抽象了实际文件，
  • 这第二名虚拟文件仅抽象奇怪的行，然后将其复制——这样在第二名将每个连续的奇怪的和甚至线具有相同的长度。

...然后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

或者更详细一点，IE编号不匹配的原始文件行：

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

相关内容