生信分析中基本Linux命令的使用

yjt2004us 2017-10-14

展开全文

内容写的特别的“简洁”，存在疑惑的部分，欢迎讨论。

Linux基本命令能做的事

学习了cat, head, tail, less, more,cut,sort,wc,uniq等基本命令后，如何使用这些命令对生物信息数据做简单的分析呢。大致可以完成以下任务：

了解数据内容
数据基本信息，例如文件大小，有多少行
数据提取，排序和去重

所以本文假定你掌握了基本的Linux命令，对于不知道的命令会用man或者help去了解这些命令的作用。

数据准备

这里采用的实验数据是拟南芥的参考基因组及其注释文件，可在TAIR中下载，命令如下：

wget http://www./download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.faswget http://www./download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

基本上从NCBI, EBI或其他数据库下载的数据都是以ASCII编码，可以用file命令检查。如果不是ASCII编码的，你需要使用hexdump或其他命令删除里面的特殊符号。

$ file TAIR10_GFF3_genes.gffTAIR10_GFF3_genes.gff: ASCII text

了解数据内容

在拿到一个纯文本文件后，第一步肯定是想看下这个文件的大致内容。但是如果在文件特别大的时候直接用cat，结果就是瞬间爆炸，啥都看不清，比较好的命令就是head,tail,less。

1.查看文件前几行:`head`

head -n 5 TAIR10_chr_all.fas

2.查看文件后几行:`tail`

tail -n 5 TAIR10_chr_all.fas

3.逐页显示文本: `less`

less TAIR10_chr_all.fas

在less显示的界面中，你可以移动光标和寻找关键字。

一些小技巧：

1.显示文件前后几行：

(head -n 2;tail -n 2) < tair10_chr_all.fas#="" 可以上述操作到.bashrc文件中作为函数function="" i()="" {=""  =""  (head="" -n="" 2;="" tail="" -n2="" )="">< '$1'="" |="" column="" -t}#　重新登录terminal或者source="">

2.去除前面的comment line：

tail -n + 2 xxxx.gff

3.调试管道命令(pipeline)：

command1 | command 2 | lesscommand1 | command 2 | head -n

对于管道命令的输出结果，可以及时使用less或者head查看，如果有错误可以及时用ctrl+c停止操作。

4.从头(de novo)管道创建：

根据第3个小技巧，我们也可以在创建多个管道的时候逐渐增加，每一步可以及时调试。

数据基本信息

查看文本数据大小

了解文本数据大小可以帮助我们简单判断处理结果，假设处理后的数据过大（好几十G）或过小（0 kb），与以往经验或期望不符，你就知道自己的处理方式存在问题了。使用ls就可以完成这个任务：

$ ls -lh TAIR10_chr_all.fa-rw-r--r-- 1 1030 users 116M Aug  8  2016 TAIR10_chr_all.fa-l 以长格式显示-h 以G,M,K为单位来显示数据大小

文件的行数

通过wc可以统计文件有多少行：

$wc -l TAIR10_chr_all.fa1514792 TAIR10_chr_all.fa

注意：实际上你可能不希望统计到commend line(以#开头的部分)以及无意义的空白行，所以你需要用grep -v排除那些无意义的行。

grep -v '#' target_file.txt | grep -v '^$' | wc -l

文件的列数

对于BED,VCF或者其他文件，你希望了解文件里包含有多少行，一个比较蠢的方法就是head -n 1然后一个一个数过去。一个比较好用的就是使用awk。这部分我将在下一篇中讲。

数据提取，排序和去重

可以使用cut提取某个特定的列，例如我只需要GFF文件的第1,3，4,5行也就是chr,feature,start,end。

$ cut -f 1,3,4,5 TAIR10_GFF3_genes.gff | head1       chromosome      1       304276711       gene    3631    58991       mRNA    3631    58991       protein 3760    56301       exon    3631    39131       five_prime_UTR  3631    37591       CDS     3760    39131       exon    3996    42761       CDS     3996    42761       exon    4486    4605# 可以保存为新的文件cut -f 1,3,4,5 TAIR10_GFF3_genes.gff | I() > part.txt

cut 可以用-d指定分隔符。

我们希望根据feature类型对part.txt文件进行进行排序：

$ sort -k2,2 part.txt | head -n31       CDS     1000112 10002311       CDS     1000112 10002311       CDS     10003966        10004523

或者是先按照chr逆序然后根据第二行排序：

$ sort -k1,1nr -k2,2 part.txt | head -n35       CDS     10001590        100017365       CDS     10004720        100048245       CDS     10004720        10004824

对于feature而言有许多相同部分，如果你想知道到底有哪几类的话，可以只提取feature，对其sort,然后统计每一个出现的次数：

$ cut -f 3 TAIR10_GFF3_genes.gff | sort | uniq -c 197160 CDS      7 chromosome 215909 exon  34621 five_prime_UTR  28775 gene    180 miRNA  35386 mRNA   3911 mRNA_TE_gene    480 ncRNA  35386 protein    924 pseudogene   1274 pseudogenic_exon    926 pseudogenic_transcript     15 rRNA     71 snoRNA     13 snRNA  30634 three_prime_UTR   3903 transposable_element_gene    689 tRNA

当然你还可以在这一步之后继续跟一个sort,找到出现最多的feature。

以上是实用一些Linux简单命令对纯文本格式数据的简单分析，之后会一篇Linux三大神器：grep,awk和sed在生物信息数据分析的应用。

敬请期待！