生信数据分析中基本Unix命令的运用

内容写的特别的“简洁”,存在疑惑的部分,可以讨论

Unix基本命令能做的事

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

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

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

数据准备

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

wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

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

$ file TAIR10_GFF3_genes.gff
TAIR10_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显示的界面中,你可以移动光标和寻找关键字

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 .bashrc就可以快捷使用了
i TAIR10_chr_all.fas

2.去除前面的comment line

tail -n + 2 xxxx.gff

3.调试管道命令(pipeline)

command1 | command 2 | less
command1 | command 2 | head -n

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

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

command1 | less
command1 | command2 | less
command1 | command2 | command3 | less

根据第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.fa
1514792 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 | head
1       chromosome      1       30427671
1       gene    3631    5899
1       mRNA    3631    5899
1       protein 3760    5630
1       exon    3631    3913
1       five_prime_UTR  3631    3759
1       CDS     3760    3913
1       exon    3996    4276
1       CDS     3996    4276
1       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 -n3
1       CDS     1000112 1000231
1       CDS     1000112 1000231
1       CDS     10003966        10004523

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

$ sort -k1,1nr -k2,2 part.txt | head -n3
5       CDS     10001590        10001736
5       CDS     10004720        10004824
5       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.

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,847评论 4 362
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,208评论 1 292
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,587评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,942评论 0 205
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,332评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,587评论 1 218
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,853评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,568评论 0 198
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,273评论 1 242
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,542评论 2 246
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,033评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,373评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,031评论 3 236
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,073评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,830评论 0 195
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,628评论 2 274
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,537评论 2 269

推荐阅读更多精彩内容