提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

转录组分析实战第一天就踩的坑——sed与换行符的恩怨

2023-02-26 09:48:38 时间

题目是这样的：

统计SRR1039510_1.fastq.gz碱基总数 1575000

由于上一题已经列出了该文件中所有的序列

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | head
TGGGAGGCTGAGGCAGGAGAATCACTTAAACCTGGGAGGCAGAGGTTACAGTGAGCCGAGATT
AAAGAAGGCGACAGTGAGAAGGAGTCCGAGAAGAGTGATGGAGACCCAATAGTCGATCCTGAG
CTGCTGGGCCCCAAGGTCCTCCTGGTCCCAGTGGTGAAGAAGGAAAGAGAGGCCCTAATGGGG
CTTGGCTGCAGCCATCCCGCTTAGCCTGCCTCACCCACACCCGTGTGGTACCTTCAGCCCTGG
TGAGACAGGTAATTCAGTATAGTAGATTAATATTTTTAATATATATTTTCCCTTAAGATTTCC
ATTTCTCAGTGTAGAAATCATGTCTTCTTAATTGCTGAACCTTACTGCAAAAACTTGTGATGT
ATCAAGAATACCAAAACAGTTTCCTAATATACAGTATTTGAAAGTGCTTGCCATATTGGCTCT
CTCATTTTCATCTTCACCATCAACAGAGAGAGCAGCATACTTGCTTGCAGAACTGAACTTAGA
TCCAACCGCAGCTTGGCATCTTCGGTGGCCTGCAGCTCGTCCTCCAGCTCTTCCAGCTGCGTC
CGGCCTCCCAAAGTGCTGGGATTACAGGCATGAGCCACCGCGCTCTGCGAGGTACTTTTTCTA

于是我想，这还不简单，直接wc -c

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | wc -c
1600000

但是，这居然和答案不对！居然多统计了25000.。。。等等，25000不就是第一题统计多少读长的答案吗

$ zcat SRR1039510_1.fastq.gz | grep '@SRR'| wc -l
25000

插播一点题外话，昨天“小二黑”直接统计行数

$ zcat SRR1039510_1.fastq.gz | wc -l
100000

老师说，再计算一下就好：每4行一个read，除以4就好了。于是我联想起之前“萌哥”讲过bc这个命令，还自己搞了一个“花样”出来：

$ echo $(zcat SRR1039510_1.fastq.gz | wc -l)
100000
$ echo $(zcat SRR1039510_1.fastq.gz | wc -l)/4
100000/4
$ echo $(zcat SRR1039510_1.fastq.gz | wc -l)/4 | bc
25000

这姑且不说。

多统计出来25000个碱基，看起来似乎是每一行多统计了一个字符，那么多出来的这个字符是什么呢？联想一下之前在某处似乎听过，每行末尾的换行符是会计算进去的。于是我cat -A了一下。

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | cat -A | head
TGGGAGGCTGAGGCAGGAGAATCACTTAAACCTGGGAGGCAGAGGTTACAGTGAGCCGAGATT$
AAAGAAGGCGACAGTGAGAAGGAGTCCGAGAAGAGTGATGGAGACCCAATAGTCGATCCTGAG$
CTGCTGGGCCCCAAGGTCCTCCTGGTCCCAGTGGTGAAGAAGGAAAGAGAGGCCCTAATGGGG$
CTTGGCTGCAGCCATCCCGCTTAGCCTGCCTCACCCACACCCGTGTGGTACCTTCAGCCCTGG$
TGAGACAGGTAATTCAGTATAGTAGATTAATATTTTTAATATATATTTTCCCTTAAGATTTCC$
ATTTCTCAGTGTAGAAATCATGTCTTCTTAATTGCTGAACCTTACTGCAAAAACTTGTGATGT$
ATCAAGAATACCAAAACAGTTTCCTAATATACAGTATTTGAAAGTGCTTGCCATATTGGCTCT$
CTCATTTTCATCTTCACCATCAACAGAGAGAGCAGCATACTTGCTTGCAGAACTGAACTTAGA$
TCCAACCGCAGCTTGGCATCTTCGGTGGCCTGCAGCTCGTCCTCCAGCTCTTCCAGCTGCGTC$
CGGCCTCCCAAAGTGCTGGGATTACAGGCATGAGCCACCGCGCTCTGCGAGGTACTTTTTCTA$

然后突然想起

$ zcat SRR1039510_1.fastq.gz |head
@SRR1039510.1 HWI-ST177:290:C0TECACXX:1:1101:1373:2104 length=63
TGGGAGGCTGAGGCAGGAGAATCACTTAAACCTGGGAGGCAGAGGTTACAGTGAGCCGAGATT
+SRR1039510.1 HWI-ST177:290:C0TECACXX:1:1101:1373:2104 length=63
HJJJIJJJJJJJJIJJJGHHIJIIIIIIJJEHGGIJGIJIJJIJHHHGGFFDFFFDEDDDBDC
@SRR1039510.2 HWI-ST177:290:C0TECACXX:1:1101:1340:2124 length=63
AAAGAAGGCGACAGTGAGAAGGAGTCCGAGAAGAGTGATGGAGACCCAATAGTCGATCCTGAG
+SRR1039510.2 HWI-ST177:290:C0TECACXX:1:1101:1340:2124 length=63
HJJJJJJJJJJJIJIIGIJJJJGJHJJJHHDFFFE@CEEEDDDDDDDDDDDDDDDBDDDDDDD
@SRR1039510.3 HWI-ST177:290:C0TECACXX:1:1101:1273:2183 length=63
CTGCTGGGCCCCAAGGTCCTCCTGGTCCCAGTGGTGAAGAAGGAAAGAGAGGCCCTAATGGGG

人家这里都告诉你每个read的长度是63，于是计算一下

$ echo 63*25000 |bc
1575000

果然就是正确答案。

说明：1，的确每行多统计了一个字符；2，且多统计的是行末的换行符。

然后我想怎么把准确的数字统计出来，那就需要把行末的换行符去掉。于是：

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | sed 's/\n//g' | wc -c
1600000

发现并没有删除掉换行符。于是陷入了各种尝试，验证，怀疑……

最后在网上搜到了一个答案：

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | sed -e ':a;N;s/\n//g;ta' | wc -c
1575001

这个。。。似乎是去掉换行符了，但是。。怎么还多一个。。。

于是又开始猜想，可能是最后一个换行符没有删掉，为什么没有删掉呢，可能这个命令不适合我，需要再调整。于是又学习了这里面的冒号：t a N都是些什么意思，怎么用的。。。。

然后发现，确实没有错，就是最后一个换行符没有删掉。

zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | sed -e ':a;N;s/\n//g;ta' | cat -A 
刷屏中……GATGGGGATCTGGCTCATGTCCAAAGGTGCAATATCAAGGAAGGGCAGGCGTGATGGCTTATTTGTTTTGTATTCAATGAATTTACACAGTCTCAGATTATCTCCTAACTCCCCACAAGACACTTATTTATTCCAAAGAAAGAACCAAAACTCCTGATGATGAAGGCCAGGGCCTCACGGGGCACCTCTCGGTTCAGGAAGAACTTCCATCGCTCGGCTGTTGGACCGGAACCAGGATGCAACTGAGGACACTGACGTGCAGAACATGAGGTGGCTTCTTGCAGTTCTCCCTGTAGACCCACGATACCAGCTGTCGGTTTTGTCAATGAAGTCCAGGCTCAGGCCCCCGGGCCCGATCATGGCCCCAGTCCGCACAGAGCGCCCGGCCCTGGCCGGTGTGGACAGAACTGGGGGCAAAAACAAAAAAGAAGGAAGGAAGGAAAGAAAGAAATGGGTAT$

N在这里把所有的序列合成1行了，然后再输出的，最后带一个换行符。

到这里彻底放弃sed了，脑雾之后想起tr还可以用，于是又尝试一下。

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | tr -d '\n' | wc -c 
1575000

终于得到了正确答案！

但是为什么sed替换掉换行符，最后一个还留着呢？为什么不用标签来写script就一个也替换不掉呢？

因为（参考：http://yysfire.github.io/linux/sed-usage-summary.html 或者https://blog.csdn.net/u011729865/article/details/71773840）

sed 维护着两个数据缓存区：模式空间和保留空间。两者均被初始化为空。 sed 对输入的每一行运行一次如下所述的执行周期：首先，sed 从输入流中读入一行，并删除行末的换行符，将此行的内容放入模式空间。然后，脚本里的命令被执行；可以对每一个命令指定地址（地址相当于一种条件，只有条件被满足，才会执行紧跟其后的命令。当到达脚本的结尾，模式空间的内容（如果之前行末的换行符被删除，此时会被加回来）被写入到输出流（除非使用了选项'-n'）。然后，对下一行开始下一个执行周期。

标签的使用参考上述地址或者：http://blog.chinaunix.net/uid-20943515-id-403.html

man一下sed，发现有个-z选项，可以用空字符作为每行的分割

-z, --null-data
​
              separate lines by NUL characters
              
              
$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | sed -z 's/\n//g'| wc -c 
1575000

其实不纠结于sed，可以用awk来实现：

$ zcat SRR1039510_1.fastq.gz | sed -n '2~4p' | awk '{print $0}' ORS='' | wc -c
1575000

相关文章