-
Notifications
You must be signed in to change notification settings - Fork 3
/
03.highlevel.Rmd
1323 lines (1029 loc) · 44.7 KB
/
03.highlevel.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Linux神器 {#LinuxGreatTools}
视频课见 <http://bioinfo.ke.qq.com>。
## 正则表达式替换文本随心所欲 {#regularExpr}
正则表达式 (regular expression)是用来做模糊匹配的,匹配符合特定模式的文本。最早来源于Unix系统中的`sed`和`grep`命令,在各个程序语言,如`perl`, `python`中也都有实现。不同程序语言中正则表达式语法大体通用,细节上又各自有自己的特色。
```{r, fig.cap="正则表达式基本语法"}
url = "http://www.ehbio.com/ehbio_resource/pyre.png"
if(!file.exists(pic_file <- "pyre.png")){
download.file(url, pic_file, mode="wb")
}
knitr::include_graphics(if(identical(knitr:::pandoc_to(),'html')) url else pic_file)
```
```bash
# 假如有这么一个测试文件
ct@ehbio: ~/$ cat <<END >url.list
http://www.ehbio.com/ImageGP
http://www.ehbio.com/Training
http://www.ehbio.com/Esx
www.ehbio.com
http://www.ehbio.com/ImageGP/index.php/Home/Index/Lineplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/GOenrichmentplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PHeatmap.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Volcanoplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Manhattanplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Histogram.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/VennDiagram.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/UpsetView.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Densityplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCoAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/CPCoAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Sankey.html
http://blog.genesino.com
https://blog.csdn.net/qazplm12_3/
https://blog.csdn.net/woodcorpse/
blog.csdn.net/woodcorpse/article/details/79313846
ImageGP is one of the Best online plot.
123456789
END
```
获取以`https`开头的行
```
ct@ehbio: ~/$ grep '^https' url.list
https://blog.csdn.net/qazplm12_3/
https://blog.csdn.net/woodcorpse/
```
获取包含数字的行
```
ct@ehbio:~/$ grep '[0-9]' url.list
https://blog.csdn.net/qazplm12_3/
blog.csdn.net/woodcorpse/article/details/79313846
123456789
```
获取空行
```
ct@ehbio:~/$ grep '^$' url.list
```
获取`html`结尾的行
```
ct@ehbio:~/$ grep 'html$' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Lineplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/GOenrichmentplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PHeatmap.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Volcanoplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Manhattanplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Histogram.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/VennDiagram.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/UpsetView.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Densityplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCoAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/CPCoAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Sankey.html
```
获取`Boxplot`和`Barplot`的地址
```
# 未能满足要求
ct@ehbio:~/$ grep 'B.*plot' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
ImageGP is one of the Best online plot.
# 一个办法:更长的匹配
ct@ehbio:~/$ grep 'B.*plot.html' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
# 限定中间不能有空格
ct@ehbio:~/$ grep 'B[^ ]*plot' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
# 限定中间只能有2个字符
ct@ehbio:~/$ grep 'B..plot' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
# 2个字符的另外一种写法
ct@ehbio:~/$ grep -P 'B.{2}plot' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
# 2个字符的再一种写法,只允许出现特定字符
ct@ehbio:~/$ grep -P 'B[arox]*plot' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/Boxplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/Barplot.html
```
获取PCA或PCoA相关的行
```
ct@ehbio:~/$ grep 'PCo*A' url.list
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/PCoAplot.html
http://www.ehbio.com/ImageGP/index.php/Home/Index/CPCoAplot.html
```
## awk-生信分析不可缺少 {#awk}
前面的学习过程中已经提到了`awk`和`sed`的使用,作为一个引子。现在则详细列举关于`awk`常用的操作和一些偏门的操作。
### awk基本参数解释 {#awk_explain}
awk擅长于对文件按行操作,每次读取一行,然后进行相应的操作。
awk读取单个文件时的基本语法格式是`awk 'BEGIN{OFS=FS="\t"}{print $0, $1;}' filename`。
读取多个文件时的语法是`awk 'BEGIN{OFS=FS="\t"}ARGIND==1{print $0, $1;}ARGIND==2{print $0;}' file1 file2`。
awk后面的命令部分是用引号括起来的,可以单引号,可以双引号,但注意不能与内部命令中用到的引号相同,否则会导致最相邻的引号视为一组,引发解释错误。**引号不可以嵌套**
* `OFS`: 文件输出时的列分隔符 (output field separtor)
* `FS`: 文件输入时的列分隔符 (field separtor)
* `BEGIN`: 设置初始参数,初始化变量
* `END`: 读完文件后做最终的处理
* 其它`{}`:循环读取文件的每一行
* `$0`表示一行内容;`$1`, `$2`, ... `$NF`表示第一列,第二列到最后一列。
* `NF (number of fields)`文件多少列;`NR (number of rows)` 文件读了多少行: `FNR` 当前文件读了多少行,常用于多文件操作时。
* `a[$1]=1`: 索引操作,类似于python中的字典,在`ID map`,`统计`中有很多应用。
### awk基本常见操作 {#awk_common_op}
* 针对特定列的计算,比如wig文件的标准化
```bash
# 注意除了第一行是空格,其它行都是tab键分割
ct@ehbio:~/sxbd$ cat <<END >ehbio.wig
variableStep chrom=chr2
300701 12
300702 10
300703 11
300704 13
300705 12.5
END
ct@localhost:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}\
{if(FNR>1) $2=$2*10^6/(2.5*10^6); print $0}' ehbio.wig
variableStep chrom=chr2
300701 4.8
300702 4
300703 4.4
300704 5.2
300705 5
```
* 计算某列内容出现的次数。
```bash
# 怎么获得count文件,应该不难吧
ct@ehbio:~/sxbd$ cat count
ID Type
Pou5f1 Pluripotency
Nanog Pluripotency
Sox2 Neuron
Tet1 Epigenetic
Tet3 Epigenetic
Myc Oncogene
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{if(FNR>1) a[$2]+=1;}END\
{print "Type\tCount"; for(i in a) print i,a[i];}' count
Type Count
Neuron 1
Epigenetic 2
Oncogene 1
Pluripotency 2
# 这个也可以用下面方式代替,但不直接
ct@ehbio:~/sxbd$ tail -n +2 count | cut -f 2 | sort | uniq -c | \
sed -e 's/^ *//' -e 's/ */\t/'
2 Epigenetic
1 Neuron
1 Oncogene
2 Pluripotency
```
* 之前也提到过的[列操作](http://mp.weixin.qq.com/s/rZ26i19hiS5ZOqIoqkL1Wg),从GTF文件中提取启动子区域
GRCh38.gtf可以从<ftp://ftp.ensembl.org/pub/release-91/gtf/homo_sapiens/Homo_sapiens.GRCh38.91.gtf.gz>下载,或使用提供的测试文件。
```bash
ct@ehbio:~/sxbd$ sed 's/"/\t/g' GRCh38.gtf | \
awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") {ensn=$10; symbol=$16; \
if($7=="+") {start=$4-1; up=start-1000; if(up<0) up=0; dw=start+500; \
print $1,up, dw, ensn, symbol, $7;} else \
if($7=="-") {start=$5-1; up=start+1000; dw=start-500; \
if(dw<0) dw=0; print $1,dw,up,ensn,symbol,$7}}}' | sort -k1,1 -k2,2n \
>GRCh38.promoter.bed
```
* 数据矩阵的格式化输出
```bash
ct@ehbio:~/sxbd$ cat numeric.matrix
ID A B C
a 1.002 1.234 1.999
b 2.333 4.232 0.889
ct@ehbio:~/sxbd$ awk '{if(FNR==1) print $0; \
else {printf "%s%s",$1,FS; for (i=2; i<=NF; i++) \
printf "%.1f %s", $i, (i==NF?RS:FS)}}' numeric.matrix
ID A B C
a 1.0 1.2 2.0
b 2.3 4.2 0.9
```
* 判断FASTQ文件中,输出质量值的长度是与序列长度不一致的序列ID
```bash
ct@ehbio:~/sxbd$ cat <<END | gzip -c >Test_2.fq.gz
>ehbio1
ACGTCGACGACGAGAGGAGAGGAGCCCTCTCGCCCGCCCTACTACCACCCACACACAACACAAGTGT
+
FFFFFFA$A#$$AFEEEEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
>ehbio2
ACGTCGACGACGAGAGGAGAGGAGCCCTCTCGCCCGCCCTACTACCACCCACACACAACACAAGTGT
+
FFFFFF$A#$$AFEEEEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
>ehbio3
ACGTCGACGACGAGAGGAGAGGAGCCTCTCGCCCGCCCTACTACCACCCACACACAACACAAGTGT
+
FFFFFFA$A#$$AFEEEEFFFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
END
```
```bash
ct@ehbio:~/sxbd$ zcat Test_2.fq.gz | \
awk '{if(FNR%4==1) ID=$0; else if(FNR%4==2) seq_len=length($0); \
else if(FNR%4==0) {quality_len=length($0); if(seq_len!=quality_len) print ID; }}'
```
* 筛选差异基因
```bash
# TAB键分割的文件
ct@ehbio:~/sxbd$ cat de_gene
ID log2fc padj
A 1 0.001
B -1 0.001
C 1 0.001
D 2 0.0001
E -0.51 0.051
F 0.1 0.1
G 1 0.1
ct@ehbio:~/sxbd$ awk '$3<0.05 || NR==1' de_gene
ID log2fc padj
A 1 0.001
B -1 0.001
C 1 0.001
D 2 0.0001
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) print $0; \
else {abs_log2fc=($2<0?$2*(-1):$2);if(abs_log2fc>=1 && $3<0.05) print $0;}}' de_gene
ID log2fc padj
A 1 0.001
B -1 0.001
C 1 0.001
D 2 0.0001
```
* 筛选差异基因存储到不同的文件
```bash
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"; up="up"; dw="dw";}\
{if(FNR==1) {print $0 >up; print $0 >dw;} else \
if ($3<0.05) {if ($2>=1) print $0 >up; else if($2<=-1) print $0 >dw;}}' de_gene
ct@ehbio:~/sxbd$ head up dw
==> up <==
ID log2fc padj
A 1 0.001
C 1 0.001
D 2 0.0001
==> dw <==
ID log2fc padj
B -1 0.001
```
* 筛选差异基因存储到不同的文件(自动版)
```bash
awk 'BEGIN{OFS=FS="\t"; up=ARGV[1]".up"; dw=ARGV[1]".dw";}\
{if(FNR==1) {print $0 >up; print $0 >dw;} \
else if ($3<=0.05) {if($2<=-1) print $0 >up; else if ($2>=1) print $0 >dw;}}' de_gene
```
* ID map,常用于转换序列的ID、提取信息、合并信息等
```bash
# TAB键分割的文件
ct@ehbio:~/sxbd$ cat id_map
ENSM Symbol Entrez
ENSG00000280516 TMEM42 693149
ENSG00000281886 TGM4 7047
ENSG00000280873 DGKD 8527
ENSG00000281244 ADAMTS13 11093
ENSG00000280701 RP11-272D20.2
ENSG00000280674 ZDHHC3 51304
ENSG00000281623 Y_RNA
ENSG00000280479 CACFD1 11094
ENSG00000281165 SLC2A6 11182
ENSG00000281879 ABO 28
ENSG00000282873 BCL7A 605
ENSG00000280651 AC156455.1 100506691
ct@ehbio:~/sxbd$ vim ensm
ct@ehbio:~/sxbd$ cat ensm
ENSG00000281244
ENSG00000281165
ENSG00000282873
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}ARGIND==1{if(FNR>1) ensm2entrez[$1]=$3;}\
ARGIND==2{print ensm2entrez[$1];}' id_map ensm
11093
11182
605
# 替代解决方案,注意 -w的使用,避免部分匹配。最稳妥的方式还是使用awk。
ct@ehbio:~/sxbd$ grep -w -f ensm id_map | cut -f 3
11093
11182
605
```
* 转换大小写, `toupper`, `tolower`
```bash
ct@ehbio:~/sxbd$ cat symbol
Tgm4
Dgkd
Abo
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}ARGIND==1{if(FNR>1) ensm2entrez[$2]=$3;}\
ARGIND==2{print ensm2entrez[toupper($1)];}' id_map symbol
7047
8527
28
```
* awk数值操作
```bash
ct@ehbio:~/sxbd$ cat <<END >file
2
4
3.1
4.5
5.4
7.6
8
END
# log2对数
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS="\t";FS="\t"}{print log($0)/log(2)}' file
# 取整,四舍五入
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS="\t";FS="\t"}{print int($1+0.5);}' file
```
* awk定义函数
```bash
ct@ehbio:~/sxbd$ cat <<END | sed 's/ */\t/g'>file
1 2 3 4
5 6 7 8
9 10 11 12
END
```
```bash
ct@ehbio:~/sxbd$ awk 'function abs(x){return ((x < 0.0) ? -x : x)}BEGIN{OFS="\t";FS="\t"}\
{pos[1]=$1;pos[2]=$2;pos[3]=$3;pos[4]=$4; len=asort(pos); \
for(i=len;i>1;i--) print abs(pos[i]-pos[i-1]);}' file
```
* 字符串匹配
```bash
# TAB键分割的文件
ct@ehbio:~/sxbd$ cat ens.bed
1 100 105
2 100 105
3 100 105
Mt 100 105
X 100 105
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{if($1~/^[0-9XY]/) $1="chr"$1; else \
if($1~/M.*/) gsub(/M.*/, "chrM", $1); print $0}' ens.bed
chr1 100 105
chr2 100 105
chr3 100 105
chrM 100 105
chrX 100 105
```
* 字符串分割
```bash
ct@ehbio:~/sxbd$ cat trinity_id
Trinity_C1_g1_i1
Trinity_C1_g1_i2
Trinity_C1_g1_i3
Trinity_C2_g1_i1
Trinity_C3_g1_i1
Trinity_C3_g3_i2
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{count=split($1, geneL, "_"); gene=geneL[1]; \
for(i=2;i<count;i++) gene=gene"_"geneL[i]; print gene,$1;}' trinity_id
Trinity_C1_g1 Trinity_C1_g1_i1
Trinity_C1_g1 Trinity_C1_g1_i2
Trinity_C1_g1 Trinity_C1_g1_i3
Trinity_C2_g1 Trinity_C2_g1_i1
Trinity_C3_g1 Trinity_C3_g1_i1
Trinity_C3_g3 Trinity_C3_g3_i2
```
* awk脚本
```bash
ct@ehbio:~/sxbd$ cat <<END >grade.awk
BEGIN{OFS=FS="\t"; up=ARGV[1]".up"; dw=ARGV[1]".dw";}
{if(FNR==1) {print $0 >up; print $0 >dw;}
else if ($3<=0.05) {
if($2<=-1) print $0 >up;
else if ($2>=1) print $0 >dw;}
}
END
ct@ehbio:~/sxbd$ awk -f grade.awk de_gene
```
* awk给每行增加行号,使其变为唯一
```bash
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS="\t";FS="\t"}NR!=1{$4=$4"_"NR;print $0}' file
```
### awk糅合操作 - 命令组合体现魅力 {#awk_combine}
* awk中执行系统命令 (注意引号的使用)
```bash
# input_mat
ct@ehbio:~/sxbd$ cat <<END | sed 's/ */\t/g' >input_mat
SRR1 root
SRR2 leaf
SRR3 stem
END
ct@ehbio:~/sxbd$ touch SRR1.fq SRR2.fq SRR3.fq
ct@ehbio:~/sxbd$ ls
SRR1.fq SRR2.fq SRR3.fq
# 系统命令组成字符串,交给system函数运行
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{system("mv "$1".fq "$2".fq");}' input_mat
#
ct@ehbio:~/sxbd$ ls
leaf.fq root.fq stem.fq
```
* awk 引用系统变量
```bash
ct@ehbio:~/sxbd$ echo 1 | awk -v ehbio="shengxinbaodian" \
-v ehbio2="sxbd" '{print ehbio, ehbio2;}'
shengxinbaodian sxbd
```
## SED命令 - 文本替换舍我其谁 {#sed}
### sed基本参数解释 {#sed_basic}
sed是`stream editor`的简称,擅长对文件进行各种正则操作、插入操作、替换操作和删除操作,可以全局,可以指定特定范围的行或者特定特征的行。
`s/pat/replace/`: 正则替换
前插行`i`, 后插行`a`, 替换行`c`, 删除行`d`, 输出行`p`
`N`: 读入下一行,同时存储;`n`:读入下一行,抛弃当前行
### 常见操作 {#sed_common}
* 替换特定的文本
```bash
# 空格是我们不太喜欢出现在文件中的一个符号,尤其是作为列名字时
# 列使用TAB键分割
ct@ehbio:~/sxbd$ cat <<END | sed 's/;/\t/g' mat
ID;2 cell;4 cell;8 cell;embryo
Pou5f1_1;2;3;4;5
Nanog_1;2;3.2;4.3;5
c-Myc;2;3;4;5
Tet1_3;2;3;4;5
END
# 一次替换
ct@ehbio:~/sxbd$ sed 's/ /_/' mat
ID 2_cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
# 全部替换
ct@ehbio:~/sxbd$ sed 's/ /_/g' mat
ID 2_cell 4_cell 8_cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
```
* 获得逗号分隔的一组数
```bash
ct@ehbio:~/sxbd$ echo `seq 1 10` | sed 's/ /,/g'
1,2,3,4,5,6,7,8,9,10
```
* 针对指定行替换
```bash
ct@ehbio:~/sxbd$ sed '2,$ s/_[0-9]//g' mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1 2 3 4 5
Nanog 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1 2 3 4 5
```
* 替换特定出现位置
```bash
# 替换第一个空格
ct@ehbio:~/sxbd$ sed 's/ /_/1' mat
ID 2_cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
# 替换第二个空格
ct@ehbio:~/sxbd$ sed 's/ /_/2' mat
ID 2 cell 4_cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
# 替换第二个及以后的空格
ct@ehbio:~/sxbd$ sed 's/ /_/2g' mat
ID 2 cell 4_cell 8_cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
```
* 给序列起名字
```bash
ct@ehbio:~/sxbd$ cat seq
ACDGTFGGCATGCDTGD
ACDGAGCDTAGCDGTA
CAGDTAGDCTADTG
ct@ehbio:~/sxbd$ sed = seq
1
ACDGTFGGCATGCDTGD
2
ACDGAGCDTAGCDGTA
3
CAGDTAGDCTADTG
# 同时缓冲两行,但只对第一行行首操作
ct@ehbio:~/sxbd$ sed = seq | sed 'N;s/^/>/;'
>1
ACDGTFGGCATGCDTGD
>2
ACDGAGCDTAGCDGTA
>3
CAGDTAGDCTADTG
```
* 给文件增加标题行
```bash
ct@ehbio:~/sxbd$ tail -n +2 mat | sort -k2,2n
c-Myc 2 3 4 5
Nanog_1 2 3.2 4.3 5
Pou5f1_1 2 3 4 5
Tet1_3 2 3 4 5
# 1 表示第一行
# i 表示插入,在指定行前面插入新行
ct@ehbio:~/sxbd$ tail -n +2 mat | sort -k2,2n | sed '1 i ID\t2_cell\t4_cell\t8_cell\tembryo'
ID 2_cell 4_cell 8_cell embryo
c-Myc 2 3 4 5
Nanog_1 2 3.2 4.3 5
Pou5f1_1 2 3 4 5
Tet1_3 2 3 4 5
```
* 提取特定或指定范围的行
```bash
# -n是必须的,阻止程序自动输出匹配行,不然会导致重复输出
ct@ehbio:~/sxbd$ sed -n '2,4p' mat
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
ct@ehbio:~/sxbd$ sed -n '4p' mat
c-Myc 2 3 4 5
```
* 提取符合特定模式的行
`/pattern/`支持普通字符串和正则表达式匹配
```bash
ct@ehbio:~/sxbd$ sed -n '/_/ p' mat
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
Tet1_3 2 3 4 5
ct@ehbio:~/sxbd$ sed -n '/-/ p' mat
c-Myc 2 3 4 5
```
* 去除文件中的空行
```bash
ct@ehbio:~/sxbd$ cat mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
# 空行就是只有行首和行尾的行
ct@ehbio:~/sxbd$ sed '/^$/d' mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
```
* 原位删除
```bash
ct@ehbio:~/sxbd$ cat mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
# -i 参数的使用
ct@ehbio:~/sxbd$ sed -i '/^$/d' mat
ct@ehbio:~/sxbd$ cat mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc 2 3 4 5
Tet1_3 2 3 4 5
```
* 删除指定范围的行
```bash
ct@ehbio:~/sxbd$ cat mat
ID 2 cell 4 cell 8 cell embryo
Pou5f1_1 2 3 4 5
Nanog_1 2 3.2 4.3 5
c-Myc_2 2 3 4 5
Tet1_3 2 3 4 5
ct@ehbio:~/sxbd$ sed '2,3d' mat
ID 2 cell 4 cell 8 cell embryo
c-Myc_2 2 3 4 5
Tet1_3 2 3 4 5
```
* 记忆匹配
`\(\)`启动记忆匹配;`\1`为第一个匹配项,`\2`为第二个匹配项;匹配项的计数根据左括号出现的位置来定,第一个`(`包括起来的为`\1`。
```bash
ct@ehbio:~/sxbd$ echo "hah ehbio hah" | sed 's/ \(.*\) /\t\1\t\1\t/'
hah ehbio ehbio hah
```
* 奇偶数行处理
```bash
ct@ehbio:~/sxbd$ echo -e "odd\neven\nodd\neven"
odd
even
odd
even
# 奇偶数行合并
ct@ehbio:~/sxbd$ echo -e "odd\neven\nodd\neven" | sed 'N;s/\n/\t/'
odd even
odd even
# 取出偶数行,比较简单
# 注意 n (小写)撇掉了奇数行
ct@ehbio:~/sxbd$ echo -e "odd\neven\nodd\neven" | sed -n 'n;p'
even
even
# 取出奇数行
# 先都读进去,然后替换偶数行为空值,再输出
ct@ehbio:~/sxbd$ echo -e "odd\neven\nodd\neven" | sed -n 'N;s/\n.*//;p'
odd
odd
```
* Windows/Linux换行符困境
Windows下的换行符是`\r\n`, Linux下换行符是`\n`, MAC下换行符是`\r`。所以Windows下的文件拷贝到Linux后,常会出现行尾多一个`^M`符号的情况,从而引起匹配或其它解析问题。
`^M`的输是 `ctrl+v+M` `ctrl+v;ctrl+m`,不是简单的输入`^`,再输入`M`。
```bash
ct@ehbio:~/sxbd$ cat -A windows.txt
ID^M$
A^M$
B^M$
C^M$
ct@ehbio:~/sxbd$ sed 's/^M//' windows.txt | cat -A
ID$
A$
B$
C$
```
* sed中使用bash变量
```bash
# 注意双引号的使用
ct@ehbio:~/sxbd$ bash_variable='ehbio'
ct@ehbio:~/sxbd$ echo "sheng xin bao dan " | sed "s/$/$bash_variable/"
sheng xin bao dan ehbio
```
## VIM的使用 {#vim}
`VIM`是一款功能强大的文本编辑工具,也是我在`Linux`,`Windows`下编辑程序和文本最常用的工具。
### 初识VIM {#vim_first}
VIM分多种状态模式,`写入`模式,`正常`模式,`可视化`模式。
* `正常`模式:打开或新建文件默认在`正常`模式,可以浏览,但不可以写入内容。这个模式也可以称作`命令行`模式,这个模式下可以使用VIM强大的命令行和快捷键功能。其它模式下按`ESC`就可以到`正常`模式。
* `写入`模式:在`正常`模式下按字母`i` (光标前插入), `o` (当前光标的下一行操作), `O` (当前光标的上一行操作),`a` (光标后插入)都可以进入`写入`模式,就可以输入内容了。
* `可视化`模式:通常用于选择特定的内容。
进入`写入`模式后,VIM使用起来可以跟`记事本`一样了。在写入文字时,可以利用组合键`CTRL+n`和`CTRL+p`完成写作单词的自动匹配补全,从而加快输入速度,保证输入的前后一致。
`正常`模式有更强大的快捷键编辑功能,把手从鼠标上解放出来。
* `dd`: 删除一行
* `3dd`: 删除一行
* `dw`: 删除一个单词
* `d3w`: 删除3个单词
* `yy`: 复制一行
* `3yy`: 复制三行
* `yw`: 复制一个单词
* `p`: (小写p)粘贴到下一行
* `P`: (大写P)粘贴到上一行
* `>>`: 当前行右缩进一个TAB
* `3>>`: 当前行及后2行都向右缩进一个TAB
* `<<`: 当前行左缩进一个TAB
* `3<<`: 当前行及后2行都向左缩进一个TAB
* `/word`: 查找特定单词
* `u`: 撤销上一次操作
* `.`: 重复上一次操作
* `CTRL+r`: 重做撤销的操作
* `y$`: 从当前复制到行尾
* `d$`: 从当前删除到行尾
跳转操作
* `gg`: 跳到文件开头
* `G`: 跳到文件结尾
* `zt`: 当前行作为可视屏幕的第一行
* `5G`: 跳到第5行
`正常`模式下输入`冒号`进入更强大的命令行定制功能。
* `:5d`: 删除第5行
* `:20,24y`:复制20到24行
* `:.,+3y`:复制当前行和下面3行
* `:2,11>`: 右缩进
* `:w`: 保存文件
* `:q`: 退出编辑器
* `:vsplit`: 分屏
键盘操作不容易被捕获,看右下角可以得到一点信息。动图请[点击查看](//blog.genesino.com/images/vim/vim_basic_operation.gif)。
```{r}
if (identical(knitr:::pandoc_to(), "html")){
url = "http://blog.genesino.com/images/vim/vim_basic_operation.gif"
if(!file.exists(vim_basic_operation <- "vim_basic_operation.gif")){
download.file(url, vim_basic_operation, mode='wb')
}
knitr::include_graphics(if(identical(knitr:::pandoc_to(), "html")) url else vim_basic_operation)
}
```
VIM还有不少**魔性**操作,具体可以看这两个帖子:
* [http://coolshell.cn/articles/5426.html](http://coolshell.cn/articles/5426.html)
* [http://coolshell.cn/articles/11312.html](http://coolshell.cn/articles/11312.html)
### VIM中使用正则表达式 {#vim_re}
这儿以提取生信宝典公众号中发过的原创文章的HTML代码为例子,获得原创文章的名字和链接,用以制作文章列表。
部分数据如下所示,利用正则表达式的第一步就是找规律。
* 这段文字是JSON格式,列表和字典的组合,使用`json`函数可以很容易解析。但我们这通过正则表达式解析。
* `title`后面跟随的文章的题目; `url`后面跟随的是文章的链接。
* `{"`和`"}`标记每篇文章的信息的开始和结束。
* `auth_apply_num`是目前不关注的信息。
```{r}
url = "http://blog.genesino.com/images/vim/wechatSXBD_source.png"
if(!file.exists(wechatSXBD_source <- "wechatSXBD_source.png")){
download.file(url, wechatSXBD_source, mode='wb')
}
knitr::include_graphics(if(identical(knitr:::pandoc_to(), "html")) url else wechatSXBD_source)
```
下面的[动画](http://blog.genesino.com/images/vim/vim_bregexpr.gif)展示了如何通过正则表达式,把这段文字只保留题目和链接,并转成`Markdown`的格式。
```{r}
if (identical(knitr:::pandoc_to(), "html")){
url = "http://blog.genesino.com/images/vim/vim_bregexpr.gif"
if(!file.exists(vim_bregexpr <- "vim_bregexpr.gif")){
download.file(url, vim_bregexpr, mode='wb')
}
knitr::include_graphics(if(identical(knitr:::pandoc_to(), "html")) url else vim_bregexpr)
}
```
* `:set wrap`: 折行显示
* `:s/"}, {"/\r/g`: `:`开启命令行模式;`s`: 是替换,之前讲Linux命令时也多次提及;`/`作为分割符,三个一起出现,前两个`/`中的内容为被替换内容,后两个`/`中的内容为替换成的内容;这里没有使用正则表达式,直接是原字符的替换,`\r`表示换行符。这样把每篇文章的信息单行显示,方便后续处理。
* `:%s/auth_apply.*"title":"/[/`:`%`表示对所有行进行操作;被替换的内容是`auth_apply`和`title":"`及其之间的内容(`.*`表示,`.`表示任意字符,`*`表示其前面的字符出现任意次)
* `:%s/".*"url":"/](/`:从题目到url之间的内容替换掉;第一次替换时忘记了第一行中开头还有引号,结果出现了误操作,后面又退回去,手动删除特殊部分,其它部分继续匹配。
* `:%s/$/)/`:表示在行尾(`$`)加上`)`, 就组成了Markdown中完整的链接形式`[context](link)`。
* `:%s/^/* /`:表示在行首(`^`)加上`* `变成Markdown格式的列表
至此就完成了生信宝典公众号文章到Markdown链接的转换,可以放到菜单栏`文章集锦`里面方便快速查询了。
一步步的处理也有些麻烦,有没有办法更简单些呢? 动画可查看[链接](http://blog.genesino.com/images/vim/vim_bregexpr_record.gif)。
```{r}
if (identical(knitr:::pandoc_to(), "html")){
url = "http://blog.genesino.com/images/vim/vim_bregexpr_record.gif"
if(!file.exists(vim_bregexpr_record <- "vim_bregexpr_record.gif")){
download.file(url, vim_bregexpr_record, mode='wb')
}
knitr::include_graphics(if(identical(knitr:::pandoc_to(), "html")) url else vim_bregexpr_record)
}
```
* 首先也是把每篇文章的信息处理为单行显示,一样的模式更容易操作,去掉第一行行首不一致的部分
* 使用上下箭头可以回溯之前的命令,类似于Linux终端下的操作
* `%s/.*title":"\([^"]*\).*url":"\(.*\)/* [\1](\2)/c`: 这个是记忆匹配,记录下匹配的内容用于替换,`\(`和`\)`表示记忆匹配的开始和结束,自身不匹配任何字符,只做标记使用;从左只右, 第一个`\(`中的内容记录为`\1`, 第二个`\(`中的内容记录为`\2`,以此类推。尤其在存在括号嵌套的情况下,注意匹配位置,左括号出现的顺序为准。在匹配文章题目时使用了`[^"]*`而不是`.*`,是考虑到正则表达式的匹配是贪婪的,会囊括更多的内容进来,就有可能出现非预期情况,所以做这么个限定,匹配所有非`"`内容。
正则表达式在数据分析中有很多灵活的应用,可以解决复杂的字符串抽提工作。常用的程序语言或命令如`pytho`, `R`, `grep`, `awk`, `sed`都支持正则表达式操作,语法也大体相似。进一步学习可参考一下链接:
* VIM正则表达式 [http://blog.csdn.net/u014015972/article/details/50688837](http://blog.csdn.net/u014015972/article/details/50688837)
* Pyton正则表达式 [https://www.cnblogs.com/huxi/archive/2010/07/04/1771073.html](https://www.cnblogs.com/huxi/archive/2010/07/04/1771073.html)
## 有了这些,文件批量重命名还需要求助其它工具吗?{#rename_all}
### 简单重命名 {#rename_simple}
Linux下文件重命名可以通过两个命令完成,`mv`和`rename`。
* `mv`: 直接运行可以进行单个文件的重命名,如 `mv old_name.txt new_name.txt`
* `rename`: 默认支持单个文件或有固定规律的一组文件的批量重命名,示例如下:
#### rename演示 {#rename_simple1}
使用`touch`新建文件,两个样品(分别是易生信a,易生信b),各自双端测序的FASTQ文件
```
ysx@ehbio:~/test$ touch YSX_a_1.fq.gz YSX_a_2.fq.gz YSX_b_2.fq.gz YSX_b_1.fq.gz
ysx@ehbio:~/test$ ls
YSX_a_1.fq.gz YSX_a_2.fq.gz YSX_b_1.fq.gz YSX_b_2.fq.gz
```
把文件名中的 易生信(`YSX`)改为易汉博 (`ehbio`)
```
# rename '被替换文字' '要替换成的文字' 操作对象
ysx@ehbio:~/test$ rename 'YSX' 'ehbio' *.gz
ysx@ehbio:~/test$ ls
ehbio_a_1.fq.gz ehbio_a_2.fq.gz ehbio_b_1.fq.gz ehbio_b_2.fq.gz
```
不同操作系统,`rename`的使用方法略有不同。印象中:
* 在CentOS都是上面的语法 `rename old new file_list`
* 在Ubuntu都是下面的语法 `rename s/old/new/ file_list`
```
# 在Centos下,该命令未起作用
ysx@ehbio:~/test$ rename 's/ehbio_//' *
ysx@ehbio:~/test$ ls
ehbio_a_1.fq.gz ehbio_a_2.fq.gz ehbio_b_1.fq.gz ehbio_b_2.fq.gz
# 如果写的rename命令没发挥作用,使用man rename查看写看其具体使用方法, 个人经验,无外乎上面提到的两种用法。
ysx@ehbio:~/test$ man rename
# NAME
# rename - rename files
#
# SYNOPSIS
# rename [options] expression replacement file...
```
替换后缀
```
# 替换后缀
ysx@ehbio:~/test$ rename 'fq' 'fastq' *.gz
ysx@ehbio:~/test$ ls