转录组分析实战第一天就踩的坑——sed与换行符的恩怨
题目是这样的:
统计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
相关文章
- Python文本情感分析_Python数据分析实战
- 在线客服系统源码开发实战总结:需求分析及前端代码基本技术方案
- Http实战之Wireshark抓包分析
- Dubbo的spi机制分析和实战案例
- Linux PCI驱动程序之MSI-X实现分析
- 【深入浅出Java原理及实战】「源码分析系列」深入分析JDK动态代理的分析原理机制
- R语言中的时间序列分析模型:ARIMA-ARCH / GARCH模型分析股票价格
- 【经验科普】实战分析C工程代码可能遇到的编译问题及其解决思路
- Vector底层结构和源码分析
- python实战:分析网站的m3u8文件下载ts文件并解密
- 技术组件优化分析:原理、方法与实战分享
- 【SpringBoot技术专题】「实战指南」从实战开发角度去分析操作RestTemplate的应用及使用技巧
- Redis实现消息的发布订阅原理分析
- Hadoop项目实战-用户行为分析之分析与设计详解大数据
- GNU/Linux程序崩溃分析框架漏洞导致内核提权风险
- Oracle数据库管理及其原理实战分析(数据库原理 oracle)
- Linux 传输文件的rz与sz命令实战分析(linux rz sz)
- MS SQL中表名与架构的必要性分析(mssql 表名 架构)
- 实战示范Oracle事务示例分析(oracle事务的例子)
- 分析Redis项目实战源码剖析(redis 项目实战源码)
- 节点Redis集群最佳实践节点数量分析(redis集群多少个)
- window系统的Rsync同步实战分析