开始我的魔幻的一天:
❝今天,我以为我发现了Linux的grep的bug,最后竟然发现是windows和linux系统的换行符不一样,知道真相的我留下了不学无术的眼泪。 ❞
事情是这个样子的:
今天,我像往常一样提取基因组的样本,我有一堆样本的ID,需要从所有的基因型的文件中提取出来。
以前我都是使用R语言,将基因型数据读进去,将所要提取的ID文件读进去,然后,我就有很多方法提取了 ,比如用match
匹配位置,然后提取写出。比如用merge
或者left_join
提取写出。比如用%in%
提取写出。
❝我有很多方法处理它,但是我今天想用
grep
函数,因为我知道grep -f file1 file2
可以根据file1的内容提取筛选file2. ❞
为什么我今天不用R语言处理了呢?
❝因为今天的基因型数据有点大,有90G,这个数据读到R中只为了筛选其中的几十行数据,不地道呀,太不地道了,虽然我们的服务器内存大,但是不是这样玩的,同事会投诉我滥用计算机资源的,我没有挖矿,为何用这么多资源?还不是编程水平差吗!想到这里,我再次流下了不学无术的眼泪。 ❞
于是,我就开始准备文件,需要提取的样本编号是这样的:
代码语言:javascript复制$ cat id1.txt
Sample_ID
202797280021_R04C03
202817020006_R03C01
202817020006_R10C03
202817020009_R02C02
202817020009_R20C01
202817170002_R07C02
202817170002_R11C01
202817170002_R11C03
202817170002_R13C03
原始的基因型数据,第一列是这样的,剩余的列都是进行数据,有1000多万位点。
说时迟那时快,我直接写下代码,是时候展示真正的实力了:
代码语言:javascript复制$ grep -f id1.txt total.txt >re_id1.txt
$ wc -l re_id1.txt
0
什么都没有!这不科学,我应该能提取出来的,应该都在文件中的,于是我用其中的一个基因型ID测试:
代码语言:javascript复制$ grep 202817020006_R10C03 total.txt
202817020006_R10C03
匹配出来了!
单个样本可以匹配出来,多个样本无法匹配出来,这是什么原因,我不仅陷入了沉思……
于是我开始了baidu,bing,google,查遍全网,也没有找到原因。
没有找到原因,我就模拟一个数据,自己测试一下吧,看看grep -f file1 file2
是不是如我理解的那样:
(base) [dfei@bogon ~]$ cat file1
2
4
e
(base) [dfei@bogon ~]$ cat file2
a1
b2
c3
d4
e5
如上所述,我模拟了两个文件,一个是另一个的子集,匹配结果如下:
代码语言:javascript复制(base) [dfei@bogon ~]$ grep -f file1 file2
b2
d4
e5
可以看到,例子是没问题的,grep -f
用起来是666的。为何我实际分析时会报错呢?我继续全网搜索。
我看了grep的参数,有一个-F
的参数,可以忽略正则表达式字符,直接用原始字符进行匹配,类似R中的fixed =T
,我好像发现了新大陆,迫不及待试了一下:
[dengfei@localhost test]$ grep -F -f id1.txt total.txt >re_id1.txt
[dengfei@localhost test]$ wc -l re_id1.txt
0
没有变化。
……漫长的分割线……
问题解决
原因是:windows和Linux换行符不一样所致。
代码语言:javascript复制cat -A name.txt
代码语言:javascript复制Sample_ID^M$
202797280021_R04C03^M$
202817020006_R03C01^M$
202817020006_R10C03^M$
202817020009_R02C02^M$
202817020009_R20C01^M$
202817170002_R07C02^M$
202817170002_R11C01^M$
202817170002_R11C03^M$
可以看到,我的文件每一行的分隔符是^M$
,这个是windows的换行符。而Linux是不支持它的,需要用dos2unix
才可以进行后续的分析。
「代码解决:」
代码语言:javascript复制dos2unix name.txt
然后再运行:
代码语言:javascript复制grep -f id1.txt total.txt >re_id1.txt
就匹配成功了。
那就顺便整理一下grep的用法吧
基本用法:
grep name list name: 为需要匹配的字符 list:为文件
1,直接查找
代码语言:javascript复制grep phoenix sample2
直接在sample2文件中,显示有phenoix的行
2,查找多个文件
代码语言:javascript复制grep phoenix sample1 sample2 sample3
在sample1,sample2,sample3三个文件中查找匹配到phoenix的行,并显示
3,查找所有文件(支持通配符)
代码语言:javascript复制grep phenoix sample*
grep phenoix *
4,忽略大小写 -i
代码语言:javascript复制grep -i phenoix *
5,递归查找,-r,查找当前文件夹的所有文件,包括所有子文件中的文件
代码语言:javascript复制grep -r phenoix *
6,反向显示 -v,显示不匹配的行
代码语言:javascript复制grep -v phenoix *
7,打印所有匹配的行,要全部匹配,而不是包含关系 -x
代码语言:javascript复制grep -x phenoix *
这里,只打印 phenoix的行,aphenoix是不打印的,因为不是完全匹配
8,仅显示匹配的文件名称,而不是所在的行 -l
代码语言:javascript复制grep -l phenoix *
9,显示匹配的个数 -c,类似uniq -c
代码语言:javascript复制grep -c phenoix *
10,显示匹配所在的行号,类似cat -n
代码语言:javascript复制grep -n phenoix *
11,匹配单词,而不是所有包含的行 -w
代码语言:javascript复制grep -w phenoix *
12,将匹配模式放到文件中 -f
代码语言:javascript复制grep -f file1 file2
会匹配file2中所有包括file1的行。注意:
- file1和file2要编码一致,可以用file查看
- file1和file2行结尾符号一致,可以用cat -A查看,如果不一致,用dos2unix命令编一下
- 匹配时,如果是单词,可以加-w,这样就不会匹配非单词
- 可以和-v一起使用,反向过滤
- -F,当有通配符是
.*
,-F将其作为固定字符,不会对其转义
❝欢迎关注我的公众号:
育种数据分析之放飞自我
。主要分享R语言,Python,育种数据分析,生物统计,数量遗传学,混合线性模型,GWAS和GS相关的知识。 ❞