內(nèi)容寫(xiě)的特別的“簡(jiǎn)潔”舶担，存在疑惑的部分坡疼，可以討論
Unix基本命令能做的事
學(xué)習(xí)了cat, head, tail, less, more,cut,sort,wc,uniq等基本命令后，如何使用這些命令對(duì)生物信息數(shù)據(jù)做簡(jiǎn)單的分析呢衣陶。大致可以完成以下任務(wù)：
- 了解數(shù)據(jù)內(nèi)容
- 數(shù)據(jù)基本信息柄瑰，例如文件大小闸氮，有多少行
- 數(shù)據(jù)提取，排序和去重
所以本文假定你掌握了基本的Unix命令教沾，對(duì)于不知道的命令會(huì)用man或者help去了解這些命令的作用蒲跨。
數(shù)據(jù)準(zhǔn)備
這里采用的實(shí)驗(yàn)數(shù)據(jù)是擬南芥的參考基因組及其注釋文件，可在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或其他數(shù)據(jù)庫(kù)下載的數(shù)據(jù)都是以ASCII編碼或悲，可以用file命令檢查。如果不是ASCII編碼的堪唐，你需要使用hexdump或其他命令刪除里面的特殊符號(hào)巡语。
$ file TAIR10_GFF3_genes.gff
TAIR10_GFF3_genes.gff: ASCII text
了解數(shù)據(jù)內(nèi)容
在拿到一個(gè)純文本文件后，第一步肯定是想看下這個(gè)文件的大致內(nèi)容淮菠。但是如果在文件特別大的時(shí)候直接用cat男公，結(jié)果就是瞬間爆炸，啥都看不清合陵，比較好的命令就是head,tail,less.
1.查看文件前幾行:head
head -n 5 TAIR10_chr_all.fas
2.查看文件后幾行:tail
tail -n 5 TAIR10_chr_all.fas
3.逐頁(yè)顯示文本: less
less TAIR10_chr_all.fas
在less顯示的界面中理澎，你可以移動(dòng)光標(biāo)和尋找關(guān)鍵字
一些小技巧：
1.顯示文件前后幾行
(head -n 2;tail -n 2) < TAIR10_chr_all.fas
# 可以上述操作到.bashrc文件中作為函數(shù)
function i() {
(head -n 2; tail -n2 ) < "$1" | column -t
}
# 重新登錄terminal或者source .bashrc就可以快捷使用了
ｉ TAIR10_chr_all.fas
2.去除前面的comment line
tail -n + 2 xxxx.gff
3.調(diào)試管道命令(pipeline)
command1 | command 2 | less
command1 | command 2 | head -n
對(duì)于管道命令的輸出結(jié)果，可以及時(shí)使用less或者h(yuǎn)ead查看曙寡，如果有錯(cuò)誤可以及時(shí)用ctrl+c停止操作
4.從頭(de novo)管道創(chuàng)建
command1 | less
command1 | command2 | less
command1 | command2 | command3 | less
根據(jù)第3個(gè)小技巧，我們也可以在創(chuàng)建多個(gè)管道的時(shí)候逐漸增加寇荧，每一步可以及時(shí)調(diào)試
數(shù)據(jù)基本信息
查看文本數(shù)據(jù)大小
了解文本數(shù)據(jù)大小可以幫助我們簡(jiǎn)單判斷處理結(jié)果举庶，假設(shè)處理后的數(shù)據(jù)過(guò)大（好幾十G）或過(guò)小（0 kb）揩抡，與以往經(jīng)驗(yàn)或期望不符户侥，你就知道自己的處理方式存在問(wèn)題了。使用ls就可以完成這個(gè)任務(wù)：
$ ls -lh TAIR10_chr_all.fa
-rw-r--r-- 1 1030 users 116M Aug 8 2016 TAIR10_chr_all.fa
-l 以長(zhǎng)格式顯示
-h 以G,M,K為單位來(lái)顯示數(shù)據(jù)大小
文件的行數(shù)
通過(guò)wc可以統(tǒng)計(jì)文件有多少行峦嗤。
$wc -l TAIR10_chr_all.fa
1514792 TAIR10_chr_all.fa
注意：實(shí)際上你可能不希望統(tǒng)計(jì)到commend line(以#開(kāi)頭的部分)以及無(wú)意義的空白行蕊唐，所以你需要用grep -v排除那些無(wú)意義的行。
grep -v "#" target_file.txt | grep -v '^$' | wc -l
文件的列數(shù)
對(duì)于BED,VCF或者其他文件烁设，你希望了解文件里包含有多少行替梨，一個(gè)比較蠢的方法就是head -n 1然后一個(gè)一個(gè)數(shù)過(guò)去。一個(gè)比較好用的就是使用awk
數(shù)據(jù)提取装黑，排序和去重
可以使用cut提取某個(gè)特定的列副瀑，例如我只需要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指定分隔符.
我們希望根據(jù)feature類型對(duì)part.txt文件進(jìn)行進(jìn)行排序
$ sort -k2,2 part.txt | head -n3
1 CDS 1000112 1000231
1 CDS 1000112 1000231
1 CDS 10003966 10004523
或者是先按照chr逆序然后根據(jù)第二行排序
$ sort -k1,1nr -k2,2 part.txt | head -n3
5 CDS 10001590 10001736
5 CDS 10004720 10004824
5 CDS 10004720 10004824
對(duì)于feature而言有許多相同部分糠睡，如果你想知道到底有哪幾類的話，可以只提取feature疚颊，對(duì)其sort,然后統(tǒng)計(jì)每一個(gè)出現(xiàn)的次數(shù)
$ 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
當(dāng)然你還可以在這一步之后繼續(xù)跟一個(gè)sort,找到出現(xiàn)最多的feature.
以上是實(shí)用一些Unix簡(jiǎn)單命令對(duì)純文本格式數(shù)據(jù)的簡(jiǎn)單分析狈孔，之后會(huì)一篇Unix三大神器：grep,awk和sed在生物信息數(shù)據(jù)分析的應(yīng)用信认。
